text:
# X-ray states and radio emission in the black hole candidate XTE J1550-564

## 1 Introduction

The Soft X-ray Transient (SXT) XTE J1550$``$564 was first detected by the All Sky Monitor (ASM) on-board the Rossi X-ray Timing Explorer (RXTE) on 1998, September 7 (Smith 1998). The outburst was characterized by a “double peaked” profile and a strong (6.8 Crab) and brief (1 day) flare. XTE J1550$``$564 went through all canonical black hole states (Sobczak et al. 1999, 2000b; Homan et al. 2001) before its return to quiescence in May 1999.

Low frequency (0.08-18 Hz) X-ray quasi-periodic oscillations (QPOs) as well as a high frequency variable (100-285 Hz) QPO were detected during some of the RXTE/PCA observations. This is the fourth BHC to display such a high frequency QPO. Based on its strong aperiodic variability, QPOs, and X-ray spectrum, XTE J1550$``$564 is believed to harbor a black hole, but a mass function has still to be measured (Cui et al. 1999; Remillard et al. 1999; Wijnands, Homan, & van der Klis 1999; Sobczak et al. 2000a and Homan et al. 2001).

Soon after its discovery in X-rays, an optical counterpart was reported by Orosz, Bailyn and Jain (1998). A brightening of 4 magnitudes in V band over the quiescence level was noted by Jain et al. (1999). Based on the interstellar absorption lines, a distance of 2.5 kpc and an optical extinction of 2.2 $`\pm `$ 0.3 magnitudes have been deduced by Sánchez-Fernández et al. (1999). A radio counterpart to XTE J1550$``$564 has been detected at 843 MHz with a flux density of $``$ 10 mJy by the Molonglo Observatory Synthesis Telescope (MOST) (Campbell-Wilson et al. 1998). Subsequent radio observations indicated that the strong X-ray flare was accompanied by a large radio flare with ejection of relativistic plasma, possibly at superluminal velocities (Hannikainen et al. 2001).

XTE J1550$``$564 became active in soft X-rays again on 2000 April 2 (MJD 51636; Smith et al. 2000). The Burst and Transient Source Experiment (BATSE) on-board the Compton Gamma Ray Observatory (CGRO) detected XTE J1550$``$564 up to 300 keV (McCollough et al. 2000). Simultaneously, a new optical outburst from XTE J1550$``$564 was reported (Masetti & Soria 2000; Jain & Bailyn 2000).

In this paper we report ATCA observations of the radio counterpart of XTE J1550$``$564, together with daily X-ray observations performed during the 2000 outburst with RXTE/ASM and CGRO/BATSE. These observations have been performed in different X-ray states and we discuss the relation of radio emission within these states.

## 2 Observations and results

### 2.1 Soft and Hard X-ray observations

The ASM on-board RXTE (Levine et al. 1996) monitors the activity in soft X-rays of any source in the sky in three energy bands (1.5–3, 3–5, 5–12 keV). The ASM light-curve in the energy band 1.5–12 keV is presented in the top panel of Figure 1 and in the lower panel, we plot the hardness ratio (5–12 keV/3–5 keV) using one-day averages. The 20–100 keV light-curve, obtained with BATSE on-board CGRO using an Earth occultation technique (Harmon et al. 1992), is presented in the middle panel of Figure 1. A power law model with a variable spectral index has been folded through the BATSE response matrix for daily flux calculation. Due to precession of the CGRO orbit, XTE J1550$``$564 stopped being occulted on MJD 51682, a few days before the end of the CGRO mission.

Figure 1 illustrates the soft and hard X-ray behavior of XTE J1550$``$564 during the 2000 X-ray outburst. During the initial rise, the soft X-ray flux increased slowly, whereas the hard X-ray flux was already strong when the orbital configuration of CGRO allowed its detection by BATSE. On April 26 (MJD 51660), we note a rapid increase in soft X-rays as well as a decrease in hard X-rays, which then settled down in a plateau until the termination of CGRO. The decay of this outburst was a gradual decline in soft X-rays.

Analysis of the RXTE/ASM hardness ratio curve and RXTE/PCA data (Miller et al. in prep) allowed us to divide the outburst into three different parts. The outburst started and finished with a hard spectrum, whereas the central part was soft, but with detectable hard X-ray emission. This can be understood as follows: an initial Low Hard state followed by a transition to the Intermediate state (or Very High state, as this state is now recognized as a higher count rate version of the Intermediate state, see Homan et al. 2001) on April 26 and then a return to the Low Hard state after May 13 (MJD 51677). Preliminary analysis of RXTE/PCA observations confirm a change in power-law photon index from 2.30 to 1.74 between May 16 and May 24 (Tomsick et al., in prep.). We note that SXTs are rarely observed in the Low Hard state in the rising phase of their outbursts.

It is interesting to note that this new X-ray outburst occurred less than a year after the discovery outburst which took place from 1998 September to 1999 May. Before that, XTE J1550$``$564 was believed to be in quiescence. This renewed activity is unusual for SXTs with long recurrence time (XTE J1550$``$564 probably belongs to this class, as it has never been detected in the past years despite sensitive surveys of the Galactic plane). A parallel can be drawn with GRO J1655$``$40 which also had a new outburst a few months after the end of series of outbursts and a period of quiescence (e.g. Hameury et al. 1997 and references therein). We should also note that both systems have relatively long orbital periods: 2.6 day for GRO J1655$``$40 (Bailyn et al. 1995) and 1.5 day for XTE J1550$``$564 (Jain et al. 2001). The spectral type (F3–F6) of its stellar companion (Orosz and Bailyn 1997) and its peculiar position in the Hertzsprung gap (Kolb et al. 1997) may be responsible for the recurrent outbursts of GRO J1655$``$40 (but see Regös, Tout, & Wickramasinghe 1998 for alternative explanations). The spectral type of the secondary in XTE J1550$``$564 has not been identified yet but photometric data are consistent with a low mass K0–K5 star (Sánchez-Fernández et al. 1999). Its precise determination would be of great interest to understand these SXTs with repeated outbursts. If XTE J1550$``$564 does have a K-type optical companion, it must be evolved to be large enough to fill its Roche lobe. An evolved companion could conceivably be the similarity that is responsible for the repeated outbursts from these systems.

### 2.2 Radio observations

The radio continuum observations of XTE J1550$``$564 were carried out using the Australia Telescope Compact Array (ATCA), located in Narrabri, NSW, Australia. The ATCA is a synthesis telescope consisting of 6 antennae, each of 22 m diameter, aligned on a 6 km east-west array. The ATCA radio observations have been conducted on three different dates in 2000: April 30, May 6 and June 1. The first two observations were performed in the 0.750D compact configuration, while the last observation was carried out in the high spatial resolution 6B array. The first and third observations were made at the central frequencies of 1384, 2496, 4800, and 8640 MHz, with a total bandwidth of 128 MHz; observations were made at only 4800 and 8640 MHz during the second observation.

The amplitude and bandpass calibrator was PKS B1934$``$638 and the antennae gain and phase calibration was derived from regular observations (every 20 minutes) of the point source calibrators B1554$``$64 (at 4800 and 8640 MHz) and B1549$``$790 (at 1384 and 2496 MHz). The editing, calibration, Fourier transformation, deconvolution and image analysis were performed using the MIRIAD package (Sault, Teuben, & Wright 1995). Due to the proximity of a very strong ($``$ 145 Jy at 1 GHz) supernova remnant (G326.3$``$1.8) close to the field of view, the shortest baselines have not been used in the analysis of the data at 1384 and 2496 MHz. The ATCA flux densities of XTE J1550$``$564 are tabulated in Table 1.

The first ATCA observation (on MJD 51665) allows a detection of the radio counterpart at a weak level of 7.45 $`\pm `$ 0.12 mJy and 5.68 $`\pm `$ 0.06 mJy at 4800 and 8640 MHz respectively. For a flux density Sν $``$ $`\nu ^\alpha `$, a spectral index of $`\alpha `$ = –0.46 $`\pm `$ 0.05 is derived. Before this first ATCA observation, the MOST performed several observations while XTE J1550$``$564 was in its initial Low Hard state and detected it at a level of 8–15 mJy at 843 MHz. The second ATCA observations performed during the Intermediate/Very High state (on MJD 51670) did not reveal significant radio emission from XTE J1550$``$564 with a strong 3 $`\sigma `$ upper limit of 0.15 mJy at 8640 MHz.

After a transition back to the Low Hard state, we performed our third ATCA observation (on MJD 51697) at nearly the end of the X-ray outburst. Radio emission is again detected at a level of $``$ 1 mJy at all four frequencies. The spectrum is slightly inverted with a spectral index of 0.37 $`\pm `$ 0.10 (Figure 2a). The best position of the radio counterpart is: $`\alpha `$(J2000) = 15h50m58.7s and $`\delta `$(J2000) = –56°28′35.″2 with a total uncertainty of 0.″3; it is consistent with the position of the optical counterpart (Jain et al. 1999).

## 3 Discussion

### 3.1 A compact jet in the Low Hard State

The detection (on June 1st) of a radio source with an inverted spectrum (spectral index of 0.37 $`\pm `$ 0.10), while XTE J1550$``$564 was in the Low Hard state, is reminiscent of the behavior of the few other BHCs (transient or persistent) which have been observed at radio frequencies in the Low Hard state. Indeed, the persistent BHCs from our Galaxy Cyg X$``$1 (Brocksopp et al. 1999), GX 339$``$4 (Corbel et al. 2000) and GRS 1758$``$258 (Lin et al. 2000) are found most of the time in the Low Hard hard state and are detected in radio with a flat (spectral index $``$ 0.0) or inverted spectrum (and also probably 1E 1740.7$``$2942). Among the SXTs, we should note that GS 2023$`+`$338 (Han & Hjellming 1992), GRO J0422$`+`$32 (Shrader et al. 1994), GS 1354$``$64 (Brocksopp et al. 2001) and XTE J1118$`+`$480 (Hynes et al. 2000, Fender et al. 2001a) have been detected with similar properties (see discussion in Fender 2001b). XTE J1550$``$564 should now be added to this list.

A flat or inverted radio spectrum can be interpreted as optically thick synchrotron emission from a compact and conical jet following the model of Hjellming & Johnston (1988). In this model, the higher frequencies come from the base of the compact jet where the optical depth is the highest. Low frequency radio emission arises from the region further from the base. Summing all contributions over the spatial extent of the jet results in a flat or inverted spectrum. Recent VLBA imaging of a compact jet in the radio core of Cyg X$``$1 (Stirling et al. 2001) confirms this interpretation. Also the stronger orbital modulation of the higher frequencies favors this scenario (Pooley, Fender, & Brocksopp 1999). The physical size of these compact jets is believed to be of the order of $``$ 10 a.u at 8.6 GHz. (corresponding to the m.a.s angular scale). More VLBA observations of BHC in the Low Hard state should increase the number of compact jet detections.

It is expected that the flat spectrum must cut off at high frequency. The spectrum of Cyg X$``$1 is flat up to 220 GHz (Fender et al. 2000a). But in the case of GX 339$``$4 (Corbel & Fender in prep.) and XTE J1118$`+`$480 (Hynes et al. 2000, Fender et al. 2001a), the inverted radio spectrum probably extends up to the near-infrared range. In Figure 2b, we show the simultaneous radio and optical observations of XTE J1550$``$564 taken on 2000 June 1 (see Jain et al. 1999, 2001 for a discussion of the YALO data). XTE J1550$``$564 is detected at V = 18.9 $`\pm `$ 0.1 mag. and I = 16.7 $`\pm `$ 0.1 mag. An optical extinction of AV = 2.2 $`\pm `$ 0.3 mag. (Sánchez-Fernández et al. 1999) has been used to deredden these data. Considering the uncertainties in the optical extinction, the spectral index between the two optical bands (–2.6 $`\pm `$ 1.0) appears to be non-thermal and is not compatible with the thermal spectrum of an optically thick accretion disk (spectral index $`>`$ +0.3 in the optical range). It is therefore possible that a significant fraction of the optical-near IR emission in XTE J1550$``$564 during the Low Hard state is (optically thin) synchrotron emission from the compact jet. If this interpretation is correct, the cut-off frequency lies in the near infrared range (as in GX 339$``$4 and XTE J1118$`+`$480).

Using the June 1 observations, the radio luminosity (1–10 GHz band) of the compact jet is $``$ 7$`\times `$10$`{}_{}{}^{28}\mathrm{erg}\mathrm{s}^1`$ for a distance of 2.5 kpc. The luminosity of the compact jet in the radio band during the first part of the outburst may be an order of magnitude higher, as XTE J1550$``$564 was detected at a level of 8–15 mJy at 843 MHz with the MOST. Assuming that the inverted spectrum extends to the near infrared ($`10^{14}`$ Hz), as discussed above, we find that the total radiative luminosity of the compact jet would be of order $``$ 2$`\times `$$`10^{34}\mathrm{erg}\mathrm{s}^1`$ for the June 1 observation, several per cent of the 2-20 keV band X-ray luminosity of 5$`\times `$$`10^{35}\mathrm{erg}\mathrm{s}^1`$ measured on the same day (Tomsick et al., in prep.). Furthermore, taking account of internal energy of particles within the jet and relativistic bulk motion is likely to increase the total power required to energize the jet (a precise estimate would require the knowledge of the inclination angle of the jet). Thus, the compact jet is likely very powerful with a total power close to the bolometric luminosity of the system (accretion disk and corona).

### 3.2 Discrete ejection events at a state transition

The first ATCA observation was performed five days after the transition from the Low Hard state to the Intermediate/Very High state. XTE J1550$``$564 is detected around $``$ 6 mJy but contrary to the Low Hard state, the radio spectrum has a spectral index of –0.46 $`\pm `$ 0.05, indicative of optically thin synchrotron emission. It is likely that this optically thin emission is due to a discrete ejection of relativistic plasma at the time of the state transition when the accretion disk is unstable. We would then be detecting the decaying radio emission from this event and the emission from near XTE J1550$``$564 would be low, consistent with our non-detection of radio emission later in the Intermediate/Very High state. During the second part of the 1998-99 X-ray outburst, Homan et al. (2001) also reported an optically thin synchrotron event after a state transition. Similar discrete ejection events have been inferred during state transitions in GX 339$``$4 (Corbel et al. 2000). It is possible that almost all X-ray binaries (including neutron star systems) produce similar weak radio events during state transitions (see also Hjellming & Han 1995 and recent work by Fender & Kuulkers 2001c).

We note that it was hypothesized that the transition to the Very High state was accompanied by massive superluminal ejections (Fender 2000b) and bright flaring radio events. These new ATCA observations rule out such an interpretation as no bright ejection event was observed. The formation of these superluminal ejection events might possibly require a transition from a quiescent state to a Very High state, i.e. a large and rapid increase of the mass accretion rate and not a transition from a state where significant accretion is already occurring (such as the Low Hard state). Indeed, the discovery outburst of XTE J1550$``$564 was a transition from quiescence to a Very High state. A bright radio event was detected after the transition and possibly resolved in two components moving at apparent superluminal velocities (Hannikainen et al. 2001).

### 3.3 Quenched radio emission in the Intermediate/Very High State

The second ATCA observation was performed when XTE J1550$``$564 was in the Intermediate/Very High state. At that time, XTE J1550$``$564 was not detected with ATCA at 4800 and 8640 MHz, with a 3 $`\sigma `$ upper limit of 0.15 mJy at 8640 MHz, i.e. a reduction of radio emission by a factor greater than 50 if we take into account the MOST detection in the initial Low Hard state. This is the first time radio observations of a BHC have been performed during an Intermediate (or Very High) state. Radio observations in this state are particularly important because usually the radio emission is dominated by the decaying optically thin synchrotron emission associated with the quiescent/active state transition (and therefore decoupled from the BH/disk system); here, we view the intrinsic radio properties of the system in this state.

This non-detection is similar to the quenching of the compact jet of GX 339$``$4 during the 1998 High Soft state (Fender et al. 1999). Following the discussion in Fender et al. (1999), it indicates that the compact jet is probably physically suppressed in the Intermediate state. These new constraints lead us to conclude that from a radio point of view the Intermediate (or Very High) and High Soft states are identical. We should note that Belloni et al. (1999) classified the 1996 outburst of Cyg X$``$1 as an Intermediate state. At that time, radio observations (only performed at the end of the Intermediate state when Cyg X$``$1 was returning to its standard Low Hard state) revealed an increase of radio emission (Zhang et al. 1997), and probably indicating that the radio emission was indeed lower in the Intermediate state, as observed here in XTE J1550$``$564.

In the Low Hard state, there is a strong coupling between the compact jet and the luminosity of the corona (see Brocksopp et al. 1999 for Cyg X$``$1 and Corbel et al. 2000 for GX 339$``$4). Despite the fact that the Intermediate state has significant hard X-ray emission, it is not a sufficient condition to sustain the emission from the compact jet (the power law index is different between the Low Hard state and the Intermediate/Very High state, but the hard X-ray emission is currently believed to originate from the corona in both states). Therefore, the Low Hard state does produces radio emission via a compact jet, whereas the states with a stronger soft component (High Soft and Intermediate/Very High) lead to the quenching of these compact jets for a reason which is not understood. Evaporation processes in the standard geometrically thin disk, below a certain accretion rate, may possibly lead to coronal outflows (e.g. Meyer, Liu, & Meyer-Hofmeister 2000). It is also possible that high magnetic energy dissipation might be responsible for the formation of outflows (and therefore a compact jet above the corona) in the Low Hard state, whereas the soft states would be dominated by viscous energy dissipation in the accretion disk (e.g. Beloborodov 1999; Poutanen 1999; Di Matteo, Celotti, & Fabian 1999).

## 4 Conclusions

Radio emission has been detected from XTE J1550$``$564 on two occasions during its 2000 outburst. Like other BHCs, observed (too rarely) in the Low Hard state, the radio spectrum of XTE J1550$``$564 is inverted. This likely implies that the radio emission is synchrotron emission from a compact jet. The total power in the compact jet may be a significant fraction of the total accretion luminosity. The Intermediate/Very High state (like the High Soft state of GX 339$``$4) is characterized by a quenching of this compact jet. It appears that Low Hard states in BHCs produce compact outflows, whereas the states with a strong soft component from the accretion disk leads to suppression of a radio emitting outflow. State transitions are accompanied with discrete ejection of relativistic plasma. More multiwavelength observations of BHCs in the Low Hard state are needed in order to understand the physics of these compact jets.

We thank Drs. Ron Ekers, Dave McConnell and the ATCA TAC for allocating us target-of-opportunity observing time, Steven Tingay for conducting one of the ATCA observations, and Dr. George Nicolson for performing HartRAO radio observations of XTE J1550$``$564 (which resulted in only upper limits due to the proximity of the SNR). S.C. would like to thank Jean-Marie Hameury and Richard Hunstead for useful discussions and informations. The Australia Telescope is funded by the Commonwealth of Australia for operation as a National Facility managed by the CSIRO. RXTE ASM results are kindly provided by the ASM/RXTE teams at MIT and at the RXTE SOF and GOF at NASA’s GSFC. SC and PK acknowledge support from NASA grant NAG5-7405.


format: text

text:
# Kronecker’s Double Series and Exact Asymptotic Expansion for Free Models of Statistical Mechanics on Torus

## I Introduction

It is well known that the singularities in thermodynamic functions associated with a critical point occur only in the thermodynamic limit when dimension $`L`$ of the system under consideration tends to infinity. In such a limit, the critical fluctuations are correlated over a distance of the order of correlation length $`\xi _{\mathrm{bulk}}`$ that may be defined as the length scale governing the exponential decay of correlation functions. Besides these two fundamental lengths, $`L`$ and $`\xi _{\mathrm{bulk}}`$, there is also the microscopic length of interactions $`a`$. Thermodynamic quantities thus may in principle depend on the dimensionless ratios $`\xi _{\mathrm{bulk}}/L`$ and $`a/L`$. The Finite-Size Scaling (FSS) hypothesis  assumes that in the scaling interval, for temperatures so close to the critical point that $`a\xi _{\mathrm{bulk}}L`$, the microscopic length drops out and the behaviour of any thermodynamic quantity can be described in terms of the universal scaling function of the scaling variable $`t=\xi _{\mathrm{bulk}}/L`$. However, non-universal corrections to FSS do exist. These sometimes can be viewed as asymptotic series in powers of $`a/L`$.

Two-dimensional models of statistical mechanics have long served as a proving ground in attempts to understand critical behavior and to test the general ideas of FSS. Very few of them have been solved exactly . Notably, the so-called free models that can be treated in terms of non-interacting quasi-particles on the lattice. The Ising model  (which is equivalent to the model of free lattice fermions) and the Gaussian model (which can be viewed as the model of free lattice bosons) are the most prominent examples.

For all free models, exact asymptotic expansion of the free energy on an infinite cylinder of circumference $`L`$ can easily be obtained by direct application of the Euler-Maclaurin summation formula . However, derivation of such an expansion on a torus of area $`S`$ and aspect ratio $`\rho `$ is much more difficult problem. For the Ising model on torus, such an expansion has first been studied by Ferdinand and Fisher . Starting with the explicit expression for the partition function , they have calculated two leading terms, $`f_{\mathrm{bulk}}`$ and $`f_0(\rho )`$, of the expansion

$$F_{T=Tc}(\rho ,S)=f_{\mathrm{bulk}}S+f_0(\rho )+\underset{p=1}{\overset{\mathrm{}}{}}f_p(\rho )S^p$$
(1)
Their calculations have recently been pushed forward to get two next sub-leading terms .

It is the purpose of this paper to derive all terms of the exact asymptotic expansion of the logarithm of the partition function on torus for a class of free exactly solvable models of statistical mechanics. Our approach is based on an intimate relation between the terms of the asymptotic expansion and the so-called Kronecker’s double series. Besides the aesthetic appeal of the exact expansion, there is also physical motivation to study non-universal corrections to FSS. The problem is that in numerical simulations of lattice models one usually studies relatively small lattices. Therefore, to compare the results of high precision numerical simulations to the theoretical predictions one cannot neglect sub-leading corrections to FSS . Non-universal terms in the asymptotic expansion also provide important information about the structure of irrelevant operators in conformal field theory .

## II Free Models of Statistical Mechanics

In this section, we formulate three basic models of statistical mechanics: Ising model, dimer model and Gaussian model. These models are often refereed to as ”free” since they were shown to be equivalent to free fermions or free bosons on the lattice. Partition functions of all these models can be written in terms of the only object — the partition function with twisted boundary conditions $`Z_{\alpha ,\beta }(\mu )`$. Exact asymptotic expansion of the later will be our main objective in the subsequent section.

### A Ising model

The Ising model is usually formulated as follows. Consider a planar square lattice of size $`M\times N`$ with periodic boundary conditions, i.e. torus. To each site $`(m,n)`$ of the torus a spin variable is ascribed, $`s_{mn}`$, with two possible values: $`+1`$ or $`1`$. Two nearest neighbor spins, say $`s_{mn}`$ and $`s_{mn+1}`$ contribute a term $`Js_{mn}s_{mn+1}`$ to the Hamiltonian, where $`J`$ is some fixed energy. Therefore, the Hamiltonian is simply the sum of all such terms, one for each edge of the lattice

$$H(s)=J\underset{n=0}{\overset{N1}{}}\underset{m=0}{\overset{M1}{}}\left(s_{mn}s_{m+1n}+s_{mn}s_{mn+1}\right)$$
(2)
The partition function of the Ising model is given by the sum over all spin configurations on the lattice

$$Z_{\mathrm{Ising}}(J)=\underset{\{s\}}{}e^{H(s)}$$

It is convenient to set up another parameterizations of the interaction constant $`J`$ in terms of the mass variable $`\mu =\mathrm{ln}\sqrt{\mathrm{sh}\mathrm{\hspace{0.17em}2}J}`$. Critical point corresponds to the massless case $`\mu =0`$.

An explicit expression for the partition function of the Ising model on $`M\times N`$ torus, which was given originally by Kaufmann , can be written as

$$Z_{\mathrm{Ising}}(\mu )=\frac{1}{2}\left(\sqrt{2}e^\mu \right)^{MN}\left\{Z_{\frac{1}{2},\frac{1}{2}}(\mu )+Z_{0,\frac{1}{2}}(\mu )+Z_{\frac{1}{2},0}(\mu )+Z_{0,0}(\mu )\right\}$$
(3)
where we have introduced the partition function with twisted boundary conditions

$$Z_{\alpha ,\beta }^2(\mu )=\underset{n=0}{\overset{N1}{}}\underset{m=0}{\overset{M1}{}}4\left[\mathrm{sin}^2\left(\frac{\pi (n+\alpha )}{N}\right)+\mathrm{sin}^2\left(\frac{\pi (m+\beta )}{M}\right)+2\mathrm{sh}^2\mu \right]$$

Here $`\alpha =0`$ corresponds to the periodic boundary conditions for the underlying free fermion in the $`N`$-direction while $`\alpha =\frac{1}{2}`$ stands for anti-periodic boundary conditions. Similarly $`\beta `$ controls boundary conditions in $`M`$-direction. With the help of the identity

$$4\left|\mathrm{sh}\left(M\omega +i\pi \beta \right)\right|^2=4\left[\mathrm{sh}^2M\omega +\mathrm{sin}^2\pi \beta \right]=\underset{m=0}{\overset{M1}{}}4\left[\mathrm{sh}^2\omega +\mathrm{sin}^2\left(\frac{\pi (m+\beta )}{M}\right)\right]$$

the partition function with twisted boundary conditions $`Z_{\alpha ,\beta }`$ can be transformed into simpler form

$$Z_{\alpha ,\beta }(\mu )=\underset{n=0}{\overset{N1}{}}2\left|\mathrm{sh}\left[M\omega _\mu \left(\frac{\pi (n+\alpha )}{N}\right)+i\pi \beta \right]\right|$$
(4)
where lattice dispersion relation has appeared

$$\omega _\mu (k)=\mathrm{arcsinh}\sqrt{\mathrm{sin}^2k+2\mathrm{sh}^2\mu }$$
(5)
This is nothing but the functional relation between energy $`\omega _\mu `$ and momentum $`k`$ of a free quasi-particle on the planar square lattice.

### B Dimer model

A ”dimer” is a two-atom molecule. When drawn on a lattice it covers two adjacent sites of the lattice and the bond that joins them. The ”dimer problem” is to determine the number of ways of covering of a given lattice with dimers, so that all sites are occupied and no two dimers overlap. If we consider planar square lattice of size $`2M\times 2N`$ wrapped on a torus then the number of dimers must be $`2MN`$ and the number of distinct dense coverings of the lattice (the partition function) has been calculated by Kasteleyn and Fisher . This can be expressed in terms of the same partition function with twisted boundary conditions Eq.(4) as

$$Z_{\mathrm{Dimer}}=\frac{1}{2}\left\{Z_{\frac{1}{2},\frac{1}{2}}^2(0)+Z_{0,\frac{1}{2}}^2(0)+Z_{\frac{1}{2},0}^2(0)Z_{0,0}^2(0)\right\}$$
(6)
First two leading terms of the asymptotic expansion of this partition function has been obtained by Ferdinand . Let us also mention that dimers are always in the critical point and have no phase transition.

### C Gaussian model

Let us now turn to a boson analog of Ising model, which is often referred to as Gaussian model. Again, we consider planar square lattice of size $`N\times M`$ wrapped on a torus. To each site $`(m,n)`$ of the lattice we assign a continuous variable $`x_{mn}`$. The Hamiltonian of the model is

$$H(x)=J\underset{n=0}{\overset{N1}{}}\underset{m=0}{\overset{M1}{}}\left(x_{mn}x_{m+1n}+x_{mn}x_{mn+1}2x_{mn}^2\right)$$
(7)
The partition function of the model can be written as

$$Z(J)=_{𝐑^{MN}}e^{H(x)}d\sigma (x)$$

If the measure $`\mathrm{d}\sigma (x)`$ in the phase space $`𝐑^{MN}`$ is chosen to be Gaussian

$$\mathrm{d}\sigma _{\mathrm{Gauss}}(x)=\pi ^{MN/2}\underset{n=0}{\overset{N1}{}}\underset{m=0}{\overset{M1}{}}e^{x_{mn}^2}\mathrm{d}x_{mn}$$

the integration can be done explicitly and the partition function of the free boson model can be written in terms of the partition function with twisted boundary conditions Eq.(4) and parameterization $`J^1=4\mathrm{ch}^2\mu `$ as

$$Z_{\mathrm{Gauss}}(\mu )=\left(\sqrt{2}\mathrm{ch}\mu \right)^{MN}\left[Z_{0,0}(\mu )\right]^1$$
(8)
This model exhibit phase transition at the point $`\mu _c=0`$ where the partition function is divergent. This is due to the presence of so-called zero mode, i.e. due to the symmetry transformation $`x_{mn}x_{mn}+\mathrm{const}`$, which leave the Hamiltonian (7) invariant. Correlation functions of disorder operator in this model have been studied by Sato, Miwa and Jimbo .

The reason why this model is often considered as boson analog of the Ising model is that, one can choose another measure in the phase space, which makes this model equivalent to the Ising model considered above

$$\mathrm{d}\sigma _{\mathrm{Ising}}(x)=2^{MN}\underset{n=0}{\overset{N1}{}}\underset{m=0}{\overset{M1}{}}\left[\delta (x_{mn}1)+\delta (x_{mn}+1)\right]\mathrm{d}x_{mn}$$

where $`\delta `$’s are Dirac $`\delta `$-functions. With such a definition the variables $`x_{mn}`$ can actually take only two values: $`+1`$ or $`1`$, so that $`x_{mn}^2=1`$. In this case integration can be replaced by summation over discrete values of $`x_{mn}=\pm 1`$ and the Hamiltonian (7) coincides with the Hamiltonian of the Ising model (2) up to a constant.

## III Asymptotic Expansion in Terms of Kronecker’s Double Series

In the previous section it was shown that partition functions of three basic free models of statistical mechanics can be expressed in terms of the only object: the partition function with twisted boundary condition $`Z_{\alpha ,\beta }(\mu )`$. In this section we shall obtain exact asymptotic $`1/S`$-expansion of the partition function near the critical point. For reader’s convenience, all the technical details of our calculations and the definitions of the special functions are summarized in the appendices attached to the paper.

First of all, let us mention the symmetry properties of the partition function $`Z_{\alpha ,\beta }(\mu )`$. From its definition (4) one can easily verify that it is even and periodic with respect to its arguments $`\alpha `$ and $`\beta `$

$`Z_{\alpha ,\beta }(\mu )`$ $`=`$ $`Z_{\alpha ,\beta }(\mu )=Z_{\alpha ,\beta }(\mu )`$ (9)
$`Z_{\alpha ,\beta }(\mu )`$ $`=`$ $`Z_{1+\alpha ,\beta }(\mu )=Z_{\alpha ,1+\beta }(\mu )`$ (10)
These imply that twist angles $`\alpha `$ and $`\beta `$ can be taken from the interval $`[0,1]`$. Then, one can note that for all twists $`(\alpha ,\beta )(0,0)`$ the partition function $`Z_{\alpha ,\beta }(\mu )`$ is even with respect to its mass argument $`\mu `$. Hence, near the critical point ($`\mu =0`$) we have

$$Z_{\alpha ,\beta }(\mu )=Z_{\alpha ,\beta }(0)+\frac{\mu ^2}{2!}Z_{\alpha ,\beta }^{\prime \prime }(0)+\mathrm{}(\alpha ,\beta )(0,0)$$
(11)
The only exception is the point where both $`\alpha `$ and $`\beta `$ are equal to zero. This case has to be treated separately since at this point the partition function turns to zero. As a result, we have

$$Z_{0,0}(\mu )=\mu Z_{0,0}^{}(0)+\frac{\mu ^3}{3!}Z_{0,0}^{\prime \prime \prime }(0)+\mathrm{}(\alpha ,\beta )=(0,0)$$
(12)
In what follows notation $`Z_{\alpha ,\beta }(\mu )`$ will imply $`(\alpha ,\beta )(0,0)`$.

### A Asymptotic Expansion of $`Z_{\alpha ,\beta }(0)`$

Considering the logarithm of the partition function with twisted boundary conditions, Eq.(4), we note, that it can be transformed as

$$\mathrm{ln}Z_{\alpha ,\beta }(0)=M\underset{n=0}{\overset{N1}{}}\omega _0\left(\frac{\pi (n+\alpha )}{N}\right)+\underset{n=0}{\overset{N1}{}}\mathrm{ln}\left|\mathrm{\hspace{0.17em}1}e^{2\left[M\omega _0\left(\frac{\pi (n+\alpha )}{N}\right)+i\pi \beta \right]}\right|$$
(13)
The second sum here vanishes in the formal limit $`M\mathrm{}`$ when the torus turns into infinitely long cylinder of circumference $`N`$. Therefore, the first sum gives the logarithm of the partition function with twisted angle $`\alpha `$ on that cylinder. Its asymptotic expansion can be found with the help of the Euler-Maclaurin summation formula (Appendix A)

$$M\underset{n=0}{\overset{N1}{}}\omega _0\left(\frac{\pi (n+\alpha )}{N}\right)=\frac{S}{\pi }_0^\pi \omega _0(x)dx\pi \lambda \rho \mathrm{B}_2^\alpha 2\pi \rho \underset{p=1}{\overset{\mathrm{}}{}}\left(\frac{\pi ^2\rho }{S}\right)^p\frac{\lambda _{2p}}{(2p)!}\frac{\mathrm{B}_{2p+2}^\alpha }{2p+2}$$
(14)
Where $`_0^\pi \omega _0(x)dx=2\gamma `$, $`\gamma =0.915965\mathrm{}`$ is Catalan’s constant and $`\mathrm{B}_p^\alpha `$ are so-called Bernoulli polynomials. We have also used the symmetry property, $`\omega _0(k)=\omega _0(\pi k)`$, of the lattice dispersion relation (5) and its Taylor expansion

$$\omega _0(k)=k\left(\lambda +\underset{p=1}{\overset{\mathrm{}}{}}\frac{\lambda _{2p}}{(2p)!}k^{2p}\right)$$
(15)
where $`\lambda =1`$, $`\lambda _2=2/3`$, $`\lambda _4=4`$, etc. In what follows, we shall not use the special values of these coefficients assuming the possibility for generalizations.

The second term in Eq.(13) we may transform as

$`{\displaystyle \underset{n=0}{\overset{N1}{}}}\mathrm{ln}\left|\mathrm{\hspace{0.17em}1}e^{2\left[M\omega _0\left(\frac{\pi (n+\alpha )}{N}\right)+i\pi \beta \right]}\right|=`$ (16)
$`=\mathrm{𝚁𝚎}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{m}}\{{\displaystyle \underset{n=0}{\overset{[N/2]1}{}}}e^{2m\left[M\omega _0\left(\frac{\pi (n+\alpha )}{N}\right)+i\pi \beta \right]}`$ (17)
$`+{\displaystyle \underset{n=0}{\overset{N[N/2]1}{}}}e^{2m\left[M\omega _0\left(\frac{\pi (n+1\alpha )}{N}\right)+i\pi \beta \right]}\}`$ (18)
The argument of the first exponent can be expanded in powers of $`1/S`$ if we replace the lattice dispersion relation $`\omega _0(x)`$ with its Taylor expansion (15)

$$\mathrm{exp}\left\{2\pi m\left[\lambda \rho (n+\alpha )+i\beta \right]2\pi m\rho \underset{p=1}{\overset{\mathrm{}}{}}\frac{\lambda _{2p}}{(2p)!}\left(\frac{\pi ^2\rho }{S}\right)^p(n+\alpha )^{2p+1}\right\}$$

Taking into account the relation between moments and cumulants (Appendix B), we obtain asymptotic expansion of the first exponent itself in powers of $`1/S`$

$`e^{2m\left[M\omega _0\left(\frac{\pi (n+\alpha )}{N}\right)+i\pi \beta \right]}=e^{2\pi m\left[\lambda \rho (n+\alpha )+i\beta \right]}`$ (19)
$`2\pi m\rho {\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\pi ^2\rho }{S}}\right)^p{\displaystyle \frac{\mathrm{\Lambda }_{2p}}{(2p)!}}(n+\alpha )^{2p+1}e^{2\pi m\left[\lambda \rho (n+\alpha )+i\beta \right]}`$ (20)
The differential operators $`\mathrm{\Lambda }_{2p}`$ that have appeared here can be expressed via coefficients $`\lambda _{2p}`$ of the expansion of the lattice dispersion relation as (15)

$`\mathrm{\Lambda }_2`$ $`=`$ $`\lambda _2`$ (21)
$`\mathrm{\Lambda }_4`$ $`=`$ $`\lambda _4+3\lambda _2^2{\displaystyle \frac{}{\lambda }}`$ (22)
$`\mathrm{\Lambda }_6`$ $`=`$ $`\lambda _6+15\lambda _4\lambda _2{\displaystyle \frac{}{\lambda }}+15\lambda _2^3{\displaystyle \frac{^2}{\lambda ^2}}`$ (23)
$`\mathrm{}`$ (24)
$`\mathrm{\Lambda }_p`$ $`=`$ $`{\displaystyle \underset{r=1}{\overset{p}{}}}{\displaystyle \left(\frac{\lambda _{p_1}}{p_1!}\right)^{k_1}\mathrm{}\left(\frac{\lambda _{p_r}}{p_r!}\right)^{k_r}\frac{p!}{k_1!\mathrm{}k_r!}\frac{^k}{\lambda ^k}}`$ (25)
Here summation is over all positive numbers $`\{k_1\mathrm{}k_r\}`$ and different positive numbers $`\{p_1,\mathrm{},p_r\}`$ such that $`p_1k_1+\mathrm{}+p_rk_r=p`$ and $`k=k_1+\mathrm{}+k_r1`$.

The expansion for the second exponent in Eq.(18) can be obtained along the same lines by substitution: $`\alpha 1\alpha `$. Plugging the expansion of both of the exponents back into Eq.(18) we obtain

$`{\displaystyle \underset{n=0}{\overset{N1}{}}}\mathrm{ln}\left|\mathrm{\hspace{0.17em}1}e^{2\left[M\omega _0\left(\frac{\pi (n+\alpha )}{N}\right)+i\pi \beta \right]}\right|=`$ (26)
$`=\mathrm{𝚁𝚎}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{m}}\left\{{\displaystyle \underset{n=0}{\overset{[N/2]1}{}}}e^{2\pi m\left[\lambda \rho \left(n+\alpha \right)+i\beta \right]}+{\displaystyle \underset{n=0}{\overset{N[N/2]1}{}}}e^{2\pi m\left[\lambda \rho \left(n+1\alpha \right)+i\beta \right]}\right\}`$ (27)
$`+2\pi \rho {\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\pi ^2\rho }{S}}\right)^p{\displaystyle \frac{\mathrm{\Lambda }_{2p}}{(2p)!}}\mathrm{𝚁𝚎}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\{{\displaystyle \underset{n=0}{\overset{[N/2]1}{}}}(n+\alpha )^{2p+1}e^{2\pi m\left[\lambda \rho (n+\alpha )+i\beta \right]}`$ (28)
$`+{\displaystyle \underset{n=0}{\overset{N[N/2]1}{}}}(n+1\alpha )^{2p+1}e^{2\pi m\left[\lambda \rho (n+1\alpha )+i\beta \right]}\}`$ (29)
In all these series, summation over $`n`$ can be extended to infinity. The resulting errors are exponentially small and do not affect our asymptotic expansion in any finite power of $`1/S`$.

The key point of our analysis is the observation that all the series that have appeared in such an expansion can be obtained by resummation of either elliptic theta function, $`\theta _{\alpha ,\beta }(\tau )`$ (Appendix C), or Kronecker’s double series, $`\mathrm{K}_p^{\alpha ,\beta }(\tau )`$ (Appendix D). Namely, with the help of the identities (C10) and (D3) we obtain

$`{\displaystyle \underset{n=0}{\overset{N1}{}}}\mathrm{ln}\left|\mathrm{\hspace{0.17em}1}e^{2\left[M\omega _0\left(\frac{\pi (n+\alpha )}{N}\right)+i\pi \beta \right]}\right|=\mathrm{ln}\left|{\displaystyle \frac{\theta _{\alpha ,\beta }(i\lambda \rho )}{\eta (i\lambda \rho )}}\right|+\pi \lambda \rho \mathrm{B}_2^\alpha `$ (30)
$`2\pi \rho {\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\pi ^2\rho }{S}}\right)^p{\displaystyle \frac{\mathrm{\Lambda }_{2p}}{(2p)!}}{\displaystyle \frac{\mathrm{𝚁𝚎}\mathrm{K}_{2p+2}^{\alpha ,\beta }(i\lambda \rho )\mathrm{B}_{2p+2}^\alpha }{2p+2}}`$ (31)
Substituting Eqs.(14) and (31) into Eq.(13) we finally obtain exact asymptotic expansion of the logarithm of the partition function with twisted boundary conditions in terms of the Kronecker’s double series

$`\mathrm{ln}Z_{\alpha ,\beta }(0)`$ $`=`$ $`{\displaystyle \frac{S}{\pi }}{\displaystyle _0^\pi }\omega _0(x)dx+\mathrm{ln}\left|{\displaystyle \frac{\theta _{\alpha ,\beta }(i\lambda \rho )}{\eta (i\lambda \rho )}}\right|`$ (33)
$`2\pi \rho {\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\pi ^2\rho }{S}}\right)^p{\displaystyle \frac{\mathrm{\Lambda }_{2p}}{(2p)!}}{\displaystyle \frac{\mathrm{𝚁𝚎}\mathrm{K}_{2p+2}^{\alpha ,\beta }(i\lambda \rho )}{2p+2}}`$
Note, that Bernoulli polynomials $`\mathrm{B}_p^\alpha `$ have finally dropped out from the asymptotic expansion on torus. This actually means that Kronecker’s double series can be considered as elliptic generalizations of Bernoulli polynomials.

### B Asymptotic Expansion of $`Z_{0,0}^{}(0)`$

As it has already been mentioned, we have to treat the case $`(\alpha ,\beta )=(0,0)`$ separately. Taking the derivative of Eq.(4) with respect to mass variable $`\mu `$ and then considering limit $`\mu 0`$ we obtain

$$Z_{0,0}^{}(0)=2\sqrt{2}M\underset{n=1}{\overset{N1}{}}2\left|\mathrm{sh}\left[M\omega _0\left(\frac{\pi n}{N}\right)\right]\right|$$

Asymptotic expansion of this expression can be found along the same lines as above. In terms of the Kronecker’s double series, the expansion can be written as

$`\mathrm{ln}Z_{0,0}^{}(0)`$ $`=`$ $`{\displaystyle \frac{S}{\pi }}{\displaystyle _0^\pi }\omega _0(x)dx+{\displaystyle \frac{1}{2}}\mathrm{ln}8\rho S+2\mathrm{ln}\left|\eta (i\lambda \rho )\right|`$ (35)
$`2\pi \rho {\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\pi ^2\rho }{S}}\right)^p{\displaystyle \frac{\mathrm{\Lambda }_{2p}}{(2p)!}}{\displaystyle \frac{\mathrm{𝚁𝚎}\mathrm{K}_{2p+2}^{0,0}(i\lambda \rho )}{2p+2}}`$
### C Asymptotic Expansion of $`Z_{\alpha ,\beta }^{\prime \prime }(0)`$

The analysis of the $`Z_{\alpha ,\beta }^{\prime \prime }(0)`$ is a little more involved. Taking the second derivative of Eq.(4) with respect to mass variable $`\mu `$ and then considering limit $`\mu 0`$ we obtain

$`{\displaystyle \frac{Z_{\alpha ,\beta }^{\prime \prime }(0)}{Z_{\alpha ,\beta }(0)}}`$ $`=`$ $`\mathrm{𝚁𝚎}M{\displaystyle \underset{n=0}{\overset{N1}{}}}\omega _0^{\prime \prime }\left(\frac{\pi (n+\alpha )}{N}\right)\mathrm{cth}\left[M\omega _0\left(\frac{\pi (n+\alpha )}{N}\right)+i\pi \beta \right]`$ (36)
$`=`$ $`M{\displaystyle \underset{n=0}{\overset{N1}{}}}\omega _0^{\prime \prime }\left(\frac{\pi (n+\alpha )}{N}\right)+2\mathrm{𝚁𝚎}M{\displaystyle \underset{n=0}{\overset{N1}{}}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\omega _0^{\prime \prime }\left(\frac{\pi (n+\alpha )}{N}\right)e^{2m\left[M\omega _0\left(\frac{\pi (n+\alpha )}{N}\right)+i\pi \beta \right]}`$ (37)
where $`\omega _0^{\prime \prime }(x)`$ is the second derivative of $`\omega _\mu (x)`$ with respect to $`\mu `$ at criticality

$`\omega _0^{\prime \prime }(x)={\displaystyle \frac{2}{\mathrm{sin}x\sqrt{1+\mathrm{sin}^2x}}}`$ (38)
Using Taylor’s theorem, the asymptotic expansion of the $`\omega _0^{\prime \prime }(x)`$ can be written in the following form

$`\omega _0^{\prime \prime }(x)={\displaystyle \frac{\kappa }{x}}\left[1+{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\kappa _{2p}}{(2p)!}}x^{2p}\right]`$ (39)
where $`\kappa =2`$, $`\kappa _2=2/3`$, $`\kappa _4=172/15`$, etc. Again, in what follows, we shall not use the special values of these coefficients assuming the possibility for generalizations.

The first sum in Eq. (36) we may transform as

$`M{\displaystyle \underset{n=0}{\overset{N1}{}}}\omega _0^{\prime \prime }\left(\frac{\pi (n+\alpha )}{N}\right)=M{\displaystyle \underset{n=0}{\overset{N1}{}}}f\left(\frac{\pi (n+\alpha )}{N}\right)+{\displaystyle \frac{\kappa S}{\pi }}{\displaystyle \underset{n=0}{\overset{N1}{}}}\left({\displaystyle \frac{1}{n+\alpha }}+{\displaystyle \frac{1}{n+1\alpha }}\right)`$ (40)
where we have introduce the function $`f(x)=\omega _0^{\prime \prime }(x)\kappa /x\kappa /(\pi x)`$. This function and all its derivatives are integrable over the interval $`(0,\pi )`$. Thus, for the first term in Eq. (40) we may use again the Euler-Maclaurin summation formula (Appendix A), and after a little algebra we obtain

$`M{\displaystyle \underset{n=0}{\overset{N1}{}}}f\left(\frac{\pi (n+\alpha )}{N}\right)`$ $`=`$ $`{\displaystyle \frac{S}{\pi }}{\displaystyle _0^\pi }f(x)dx\pi \rho \kappa {\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\pi ^2\rho }{S}}\right)^{p1}{\displaystyle \frac{\kappa _{2p}\mathrm{B}_{2p}^\alpha }{p(2p)!}}`$ (41)
$`+`$ $`{\displaystyle \frac{\kappa S}{\pi }}{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{B}_{2p}^\alpha }{p}}{\displaystyle \frac{1}{N^{2p}}}`$ (42)
where $`_0^\pi f(x)dx=2\mathrm{ln}24\mathrm{ln}\pi `$. The second sum in Eq. (40) can be written in terms of the digamma function $`\psi (x)`$.

$`{\displaystyle \underset{n=0}{\overset{N1}{}}}\left({\displaystyle \frac{1}{n+\alpha }}+{\displaystyle \frac{1}{n+1\alpha }}\right)=\left[\psi (N+\alpha )+\psi (N+1\alpha )\psi (\alpha )\psi (1\alpha )\right]`$ (43)
The asymptotic expansion of the digamma function $`\psi (x)`$ is given by (see Appendix E)

$$\psi (N+\alpha )=\mathrm{ln}N\underset{p=1}{\overset{\mathrm{}}{}}(1)^p\frac{\mathrm{B}_p^\alpha }{p}\frac{1}{N^p}$$
(44)
Using the symmetry properties of the Bernoulli polynomials $`\mathrm{B}_p^\alpha `$, namely, $`\mathrm{B}_{2p}^{1\alpha }=\mathrm{B}_{2p}^\alpha `$ and $`\mathrm{B}_{2p+1}^{1\alpha }=\mathrm{B}_{2p+1}^\alpha `$, the Eq. (43) can be rewritten as

$`{\displaystyle \underset{n=0}{\overset{N1}{}}}\left({\displaystyle \frac{1}{n+\alpha }}+{\displaystyle \frac{1}{n+1\alpha }}\right)=2\mathrm{ln}N{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{B}_{2p}^\alpha }{p}}{\displaystyle \frac{1}{N^{2p}}}\psi (\alpha )\psi (1\alpha )`$ (45)
Plugging Eqs. (42) and (45) back in Eq. (40) we have finally obtain

$`M{\displaystyle \underset{n=0}{\overset{N1}{}}}\omega _0^{\prime \prime }\left(\frac{\pi (n+\alpha )}{N}\right)`$ $`=`$ $`{\displaystyle \frac{2\kappa S}{\pi }}\left[{\displaystyle \frac{1}{2\kappa }}{\displaystyle _0^\pi }f(x)dx+\mathrm{ln}N{\displaystyle \frac{\psi (\alpha )+\psi (1\alpha )}{2}}\right]`$ (46)
$``$ $`\pi \rho \kappa {\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\pi ^2\rho }{S}}\right)^{p1}{\displaystyle \frac{\kappa _{2p}\mathrm{B}_{2p}^\alpha }{p(2p)!}}`$ (47)
Let us now consider the second sum in Eq. (36). Note that function $`\omega _0^{\prime \prime }(x)`$ can be represent as

$$\omega _0^{\prime \prime }(x)=\frac{\kappa }{x}\mathrm{exp}\left\{\underset{p=1}{\overset{\mathrm{}}{}}\frac{\epsilon _{2p}}{(2p)!}x^{2p}\right\}$$
(48)
where coefficients $`\epsilon _{2p}`$ and $`\kappa _{2p}`$ are related to each other through relation between moments and cumulants (Appendix B). Following along the same lines as in the section (3.1), the second sum in Eq. (36) can be written as

$`2\mathrm{𝚁𝚎}M{\displaystyle \underset{n=0}{\overset{N1}{}}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\omega _0^{\prime \prime }\left(\frac{\pi (n+\alpha )}{N}\right)e^{2m\left[M\omega _0\left(\frac{\pi (n+\alpha )}{N}\right)+i\pi \beta \right]}=`$ (49)
$`={\displaystyle \frac{2\kappa S}{\pi }}\mathrm{𝚁𝚎}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\left\{{\displaystyle \frac{1}{n+\alpha }}e^{2\pi m\left[\rho \left(n+\alpha \right)+i\beta \right]}+{\displaystyle \frac{1}{n+1\alpha }}e^{2\pi m\left[\rho \left(n+1\alpha \right)+i\beta \right]}\right\}`$ (50)
$`+\kappa \pi \rho \mathrm{\Omega }_2\mathrm{𝚁𝚎}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\left\{(n+\alpha )e^{2\pi m\left[\lambda \rho \left(n+\alpha \right)+i\beta \right]}+(n+1\alpha )e^{2\pi m\left[\lambda \rho \left(n+1\alpha \right)+i\beta \right]}\right\}`$ (51)
$`\kappa \pi \rho {\displaystyle \underset{p=2}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{\Omega }_{2p}}{p(2p)!}}\left({\displaystyle \frac{\pi ^2\rho }{S}}\right)^{p1}\mathrm{𝚁𝚎}\mathrm{K}_{2p}^{\alpha ,\beta }(i\lambda \rho )+\kappa \pi \rho {\displaystyle \underset{p=2}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\kappa _{2p}\mathrm{B}_{2p}^\alpha }{p(2p)!}}\left({\displaystyle \frac{\pi ^2\rho }{S}}\right)^{p1}`$ (52)
The differential operators $`\mathrm{\Omega }_{2p}`$ that have appeared here can be expressed via coefficients $`\omega _{2p}=\epsilon _{2p}+\lambda _{2p}\frac{}{\lambda }`$ as

$`\mathrm{\Omega }_2`$ $`=`$ $`\omega _2`$ (53)
$`\mathrm{\Omega }_4`$ $`=`$ $`\omega _4+3\omega _2^2`$ (54)
$`\mathrm{}`$ (55)
Let us introduce the function $`R_{\alpha ,\beta }(\rho )`$

$`R_{\alpha ,\beta }(\rho )`$ $`=`$ $`{\displaystyle \frac{\psi (\alpha )+\psi (1\alpha )}{2}}`$ (56)
$`+`$ $`\mathrm{𝚁𝚎}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\left\{{\displaystyle \frac{1}{n+\alpha }}e^{2\pi m\left[\rho \left(n+\alpha \right)+i\beta \right]}+{\displaystyle \frac{1}{n+1\alpha }}e^{2\pi m\left[\rho \left(n+1\alpha \right)+i\beta \right]}\right\}`$ (57)
which first derivative with respect to $`\rho `$ is given by (see Appendix C, Eq. (C14))

$`{\displaystyle \frac{}{\rho }}R_{\alpha ,\beta }(\rho )=2{\displaystyle \frac{}{\rho }}\mathrm{ln}\left|\theta _{\alpha ,\beta }(i\rho )\right|+{\displaystyle \frac{1}{2\pi }}\left({\displaystyle \frac{}{z}}\mathrm{ln}\left|\theta _{\alpha ,\beta }(i\rho )\right|\right)^2`$ (58)
For the cases $`(\alpha ,\beta )=(0,1/2),(1/2,0),(1/2,1/2)`$ the second term in above equation is equal to zero, and for $`R_{\alpha ,\beta }(\rho )`$ we obtain

$`R_{\alpha ,\beta }(\rho )=2\mathrm{ln}\left|\theta _{\alpha ,\beta }(i\rho )\right|+C_E+2\mathrm{ln}2`$ (59)
where $`C_E`$ is the Euler constant.

With the help of the identity (C12) the Eq. (52) can be rewritten as

$`2\mathrm{𝚁𝚎}M{\displaystyle \underset{n=0}{\overset{N1}{}}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\omega _0^{\prime \prime }\left(\frac{\pi (n+\alpha )}{N}\right)e^{2m\left[M\omega _0\left(\frac{\pi (n+\alpha )}{N}\right)+i\pi \beta \right]}=`$ (60)
$`={\displaystyle \frac{2\kappa S}{\pi }}\left[R_{\alpha ,\beta }(\rho )+{\displaystyle \frac{\psi (\alpha )+\psi (1\alpha )}{2}}\right]+{\displaystyle \frac{\kappa }{2}}\left(\kappa _2\rho {\displaystyle \frac{}{\rho }}+\lambda _2\rho ^2{\displaystyle \frac{^2}{\rho ^2}}\right)\mathrm{ln}\left|{\displaystyle \frac{\theta _{\alpha ,\beta }(i\rho )}{\eta (i\rho )}}\right|`$ (61)
$`\kappa \pi \rho {\displaystyle \underset{p=2}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{\Omega }_{2p}}{p(2p)!}}\left({\displaystyle \frac{\pi ^2\rho }{S}}\right)^{p1}\mathrm{𝚁𝚎}\mathrm{K}_{2p}^{\alpha ,\beta }(i\lambda \rho )+\kappa \pi \rho {\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\kappa _{2p}\mathrm{B}_{2p}^\alpha }{p(2p)!}}\left({\displaystyle \frac{\pi ^2\rho }{S}}\right)^{p1}`$ (62)
Substituting Eqs. (47) and (62) into Eq. (36) we finally obtain exact asymptotic expansion of the $`Z_{\alpha ,\beta }^{\prime \prime }(0)`$

$`{\displaystyle \frac{Z_{\alpha ,\beta }^{\prime \prime }(0)}{Z_{\alpha ,\beta }(0)}}`$ $`=`$ $`{\displaystyle \frac{2\kappa S}{\pi }}\left[{\displaystyle \frac{1}{2\kappa }}{\displaystyle _0^\pi }f(x)dx+\mathrm{ln}N+R_{\alpha ,\beta }(\rho )\right]`$ (63)
$`+`$ $`{\displaystyle \frac{\kappa }{2}}\left(\kappa _2\rho {\displaystyle \frac{}{\rho }}+\lambda _2\rho ^2{\displaystyle \frac{^2}{\rho ^2}}\right)\mathrm{ln}\left|{\displaystyle \frac{\theta _{\alpha ,\beta }(i\rho )}{\eta (i\rho )}}\right|`$ (64)
$``$ $`\kappa \pi \rho {\displaystyle \underset{p=2}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{\Omega }_{2p}}{p(2p)!}}\left({\displaystyle \frac{\pi ^2\rho }{S}}\right)^{p1}\mathrm{𝚁𝚎}\mathrm{K}_{2p}^{\alpha ,\beta }(i\lambda \rho )`$ (65)
### D Asymptotic Expansion of $`Z_{0,0}^{\prime \prime \prime }(0)`$

Let us now consider the case $`(\alpha ,\beta )=(0,0)`$. Eq. (12) implies immediately that

$$\underset{\mu 0}{lim}\frac{Z_{0,0}^{\prime \prime }(\mu )}{Z_{0,0}(\mu )}=\frac{Z_{0,0}^{\prime \prime \prime }(0)}{Z_{0,0}^{}(0)}$$
(66)
Taking the third derivative of Eq.(4) with respect to mass variable $`\mu `$ and then considering limit $`\mu 0`$ we obtain

$`{\displaystyle \frac{Z_{0,0}^{\prime \prime \prime }(0)}{Z_{0,0}^{}(0)}}`$ $`=`$ $`2M^21+3M{\displaystyle \underset{n=1}{\overset{N1}{}}}\omega _0^{\prime \prime }\left(\frac{\pi n}{N}\right)\mathrm{cth}\left[M\omega _0\left(\frac{\pi n}{N}\right)\right]`$ (67)
Asymptotic expansion of the $`Z_{0,0}^{\prime \prime \prime }(0)`$ can be found along the same lines as above. In terms of the Kronecker’s double series, the expansion can be written as

$`{\displaystyle \frac{Z_{0,0}^{\prime \prime \prime }(0)}{Z_{0,0}^{}(0)}}`$ $`=`$ $`{\displaystyle \frac{6\kappa S}{\pi }}\left[{\displaystyle \frac{1}{2\kappa }}{\displaystyle _0^\pi }f(x)dx+\mathrm{ln}N+C_E\mathrm{ln}\eta (i\rho )\right]`$ (68)
$`+`$ $`3\kappa \left(\kappa _2\rho {\displaystyle \frac{}{\rho }}+\lambda _2\rho ^2{\displaystyle \frac{^2}{\rho ^2}}\right)\mathrm{ln}\eta (i\rho )1`$ (69)
$``$ $`3\kappa \pi \rho {\displaystyle \underset{p=2}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{\Omega }_{2p}}{p(2p)!}}\left({\displaystyle \frac{\pi ^2\rho }{S}}\right)^{p1}\mathrm{K}_{2p}^{0,0}(i\lambda \rho )`$ (70)
Expansions (33), (35), (64) and (70) are the main results of the paper. Kronecker’s double series $`\mathrm{K}_p^{\alpha ,\beta }`$ with $`\alpha `$ and $`\beta `$ taking values $`0`$ and $`1/2`$ can all be expressed in terms of the elliptic $`\theta `$-functions only (Appendix F).

## IV Asymptotic Expansion of the Free Energy, <br>the Internal Energy and the Specific Heat

After reaching this point, one can easily write down all the terms of the exact asymptotic expansion (1) of the free energy, $`F=\mathrm{ln}Z`$, at the critical point for all three models under consideration. For the Ising model and the Gaussian model we have found that the exact asymptotic expansion of the internal energy, $`U=\mathrm{ln}Z(J)/J`$, and the specific heat, $`C=^2\mathrm{ln}Z(J)/J^2`$, at the critical point can be written in the following form

$`U`$ $`=`$ $`u_{\mathrm{bulk}}S+{\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}u_p(\rho )S^{p+\frac{1}{2}}`$ (71)
$`C`$ $`=`$ $`c_{\mathrm{bulk}}S+{\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}c_{2p}(\rho )S^{p+\frac{1}{2}}+{\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}c_{2p+1}(\rho )S^p`$ (72)
1. The asymptotic expansion of the free energy of the Ising model is

$`f_{\mathrm{bulk}}`$ $`=`$ $`\mathrm{ln}\sqrt{2}{\displaystyle \frac{2\gamma }{\pi }}`$ (73)
$`f_0(\rho )`$ $`=`$ $`\mathrm{ln}{\displaystyle \frac{\theta _2+\theta _3+\theta _4}{2\eta }}`$ (74)
$`f_1(\rho )`$ $`=`$ $`{\displaystyle \frac{\pi ^3\rho ^2}{180}}{\displaystyle \frac{\frac{7}{8}(\theta _2^9+\theta _3^9+\theta _4^9)+\theta _2\theta _3\theta _4\left[\theta _2^3\theta _4^3\theta _3^3\theta _2^3\theta _3^3\theta _4^3\right]}{\theta _2+\theta _3+\theta _4}}`$ (75)
$`f_2(\rho )`$ $`=`$ $`{\displaystyle \frac{\pi ^6\rho ^5}{18432}}{\displaystyle \frac{\theta _2\theta _3\theta _4\left\{\theta _3^7(\theta _2^4\theta _4^4)^2+\theta _2^7(\theta _3^4+\theta _4^4)^2+\theta _4^7(\theta _3^4+\theta _2^4)^2\right\}}{(\theta _2+\theta _3+\theta _4)^2}}`$ (76)
$``$ $`{\displaystyle \frac{31\pi ^5\rho ^4}{24192}}{\displaystyle \frac{\theta _3^4(\theta _4^9\theta _2^9)+\theta _2^4(\theta _4^9\theta _3^9)+\theta _4^4(\theta _3^9\theta _2^9)}{\theta _2+\theta _3+\theta _4}}\left(1+4\rho {\displaystyle \frac{}{\rho }}\mathrm{ln}\theta _2\right)`$ (77)
$``$ $`{\displaystyle \frac{\pi ^5\rho ^4}{1512}}{\displaystyle \frac{\theta _2\theta _3\theta _4\left[\theta _3^3(\theta _2^7\theta _4^7)+\theta _2^3(\theta _3^7\theta _4^7)+\theta _4^3(\theta _2^7\theta _3^7)\right]}{\theta _2+\theta _3+\theta _4}}\left(1+4\rho {\displaystyle \frac{}{\rho }}\mathrm{ln}\theta _2\right)`$ (78)
$``$ $`{\displaystyle \frac{31\pi ^6\rho ^5}{72576}}{\displaystyle \frac{\theta _3^4\theta _4^4\left[3\theta _3^4\theta _4^4(\theta _3+\theta _4)\theta _2^4(2\theta _2^5+\theta _3^5\theta _4^5)\right]}{\theta _2+\theta _3+\theta _4}}`$ (79)
$``$ $`{\displaystyle \frac{\pi ^6\rho ^5}{4536}}{\displaystyle \frac{\theta _3^4\theta _4^4\left[\theta _2\theta _3(\theta _2^7\theta _3^7)+\theta _2\theta _4(\theta _2^7\theta _4^7)+2\theta _2^4\theta _3\theta _4(\theta _3^3\theta _4^3)4\theta _2\theta _3^4\theta _4^4\right]}{\theta _2+\theta _3+\theta _4}}`$ (80)
$`\mathrm{}`$ (81)
$`u_{\mathrm{bulk}}`$ $`=`$ $`\sqrt{2}`$ (82)
$`u_0(\rho )`$ $`=`$ $`2\sqrt{\rho }{\displaystyle \frac{\theta _2\theta _3\theta _4}{\theta _2+\theta _3+\theta _4}}`$ (83)
$`u_1(\rho )`$ $`=`$ $`{\displaystyle \frac{\pi ^3\rho ^{5/2}}{48}}{\displaystyle \frac{\theta _2\theta _3\theta _4(\theta _2^9+\theta _3^9+\theta _4^9)}{(\theta _2+\theta _3+\theta _4)^2}}`$ (84)
$`u_2(\rho )`$ $`=`$ $`{\displaystyle \frac{\pi ^6\rho ^{9/2}}{4608}}{\displaystyle \frac{\theta _2^2\theta _3^2\theta _4^2\left[\theta _3^7(\theta _2^4\theta _4^4)^2+\theta _2^7(\theta _3^4+\theta _4^4)^2+\theta _4^7(\theta _3^4+\theta _2^4)^2\right]}{(\theta _2+\theta _3+\theta _4)^3}}`$ (85)
$`+`$ $`{\displaystyle \frac{\pi ^6\rho ^{9/2}}{9216}}{\displaystyle \frac{\theta _2\theta _3\theta _4\left[23(\theta _2^{17}+\theta _3^{17}+\theta _4^{17})+8\theta _2^4\theta _4^4(5\theta _2^9+5\theta _4^98\theta _3^9+2\theta _2^4\theta _4^5+2\theta _2^5\theta _4^4)\right]}{(\theta _2+\theta _3+\theta _4)^2}}`$ (86)
$``$ $`{\displaystyle \frac{\pi ^5\rho ^{7/2}}{192}}{\displaystyle \frac{\theta _2\theta _3\theta _4\left[\theta _3^9(\theta _2^4\theta _4^4)+\theta _2^9(\theta _3^4+\theta _4^4)\theta _4^9(\theta _3^4+\theta _2^4)\right]}{(\theta _2+\theta _3+\theta _4)^2}}\left(1+{\displaystyle \frac{\pi \rho }{2}}\theta _3^4+4\rho {\displaystyle \frac{}{\rho }}\mathrm{ln}\theta _2\right)`$ (87)
$`\mathrm{}`$ (88)
$`c_{\mathrm{bulk}}(\rho )`$ $`=`$ $`{\displaystyle \frac{8}{\pi }}\left(\mathrm{ln}\sqrt{{\displaystyle \frac{S}{\rho }}}+\mathrm{ln}{\displaystyle \frac{2^{5/2}}{\pi }}+C_E{\displaystyle \frac{\pi }{4}}\right)4\rho \left({\displaystyle \frac{\theta _2\theta _3\theta _4}{\theta _2+\theta _3+\theta _4}}\right)^2{\displaystyle \frac{16}{\pi }}{\displaystyle \frac{_{i=2}^4\theta _i\mathrm{ln}\theta _i}{\theta _2+\theta _3+\theta _4}}`$ (89)
$`c_0(\rho )`$ $`=`$ $`2\sqrt{2}\sqrt{\rho }{\displaystyle \frac{\theta _2\theta _3\theta _4}{\theta _2+\theta _3+\theta _4}}`$ (90)
$`c_1(\rho )`$ $`=`$ $`{\displaystyle \frac{\pi ^2\rho ^2}{6}}{\displaystyle \frac{(\theta _2^8\theta _3^8)\theta _2\theta _3\mathrm{ln}\frac{\theta _3}{\theta _2}+(\theta _4^8\theta _3^8)\theta _4\theta _3\mathrm{ln}\frac{\theta _3}{\theta _4}+(\theta _2^8\theta _4^8)\theta _4\theta _2\mathrm{ln}\frac{\theta _4}{\theta _2}}{(\theta _2+\theta _3+\theta _4)^2}}`$ (91)
$`+`$ $`{\displaystyle \frac{\pi ^2\rho ^2}{9}}{\displaystyle \frac{\theta _3^4\theta _4^4(2\theta _2\theta _3\theta _4)}{\theta _2+\theta _3+\theta _4}}+{\displaystyle \frac{\pi ^3\rho ^3}{12}}{\displaystyle \frac{\theta _2^2\theta _3^2\theta _4^2(\theta _2^9+\theta _3^9+\theta _4^9)}{(\theta _2+\theta _3+\theta _4)^3}}`$ (92)
$`+`$ $`{\displaystyle \frac{\pi \rho }{9}}{\displaystyle \frac{\theta _2^5\theta _4^5+\theta _3(\theta _2^4\theta _4^4)2\theta _2\theta _4(\theta _2^3\theta _4^3)}{\theta _2+\theta _3+\theta _4}}\left(1+4\rho {\displaystyle \frac{}{\rho }}\mathrm{ln}\theta _2\right)`$ (93)
$`c_2(\rho )`$ $`=`$ $`{\displaystyle \frac{\pi ^3\rho ^{5/2}}{24\sqrt{2}}}{\displaystyle \frac{\theta _2\theta _3\theta _4(\theta _2^9+\theta _3^9+\theta _4^9)}{(\theta _2+\theta _3+\theta _4)^2}}`$ (94)
$`\mathrm{}`$ (95)
We have also used the following relations between derivatives of the elliptic functions

$$\frac{}{\rho }\mathrm{ln}\theta _3=\frac{\pi }{4}\theta _4^4+\frac{}{\rho }\mathrm{ln}\theta _2\text{and}\frac{}{\rho }\mathrm{ln}\theta _4=\frac{\pi }{4}\theta _3^4+\frac{}{\rho }\mathrm{ln}\theta _2$$

Note, that with the help of the identities

$$\frac{}{\rho }\mathrm{ln}\theta _2=\frac{1}{2}\theta _3^2E,\text{and}\frac{E}{\rho }=\frac{\pi ^2}{4}\theta _3^2\theta _4^4\frac{\pi }{2}\theta _4^4E$$

one can express all derivatives of the elliptic functions in terms of the elliptic functions $`\theta _2,\theta _3,\theta _4`$ and the elliptic integral of the second kind $`E`$
2. Similar expansion of the free energy of the dimer model is

$`f_{\mathrm{bulk}}`$ $`=`$ $`{\displaystyle \frac{\gamma }{\pi }}`$ (96)
$`f_0(\rho )`$ $`=`$ $`\mathrm{ln}{\displaystyle \frac{\theta _2^2+\theta _3^2+\theta _4^2}{2\eta ^2}}`$ (97)
$`f_1(\rho )`$ $`=`$ $`{\displaystyle \frac{\pi ^3\rho ^2}{90}}{\displaystyle \frac{\frac{7}{8}(\theta _2^{10}+\theta _3^{10}+\theta _4^{10})+\theta _2^2\theta _3^2\theta _4^2\left[\theta _2^2\theta _4^2\theta _2^2\theta _3^2\theta _3^2\theta _4^2\right]}{\theta _2^2+\theta _3^2+\theta _4^2}}`$ (98)
$`f_2(\rho )`$ $`=`$ $`{\displaystyle \frac{\pi ^6\rho ^4}{16200}}{\displaystyle \frac{\theta _3^2\left(\frac{7}{8}\theta _3^8+\theta _2^4\theta _4^4\right)^2+\theta _2^2\left(\frac{7}{8}\theta _2^8\theta _3^4\theta _4^4\right)^2+\theta _4^2\left(\frac{7}{8}\theta _4^8\theta _3^4\theta _2^4\right)^2}{\theta _2^2+\theta _3^2+\theta _4^2}}`$ (99)
$`+`$ $`{\displaystyle \frac{\pi ^6\rho ^4}{16200}}\left[{\displaystyle \frac{\frac{7}{8}(\theta _2^{10}+\theta _3^{10}+\theta _4^{10})+\theta _2^2\theta _3^2\theta _4^2\left[\theta _2^2\theta _4^2\theta _2^2\theta _3^2\theta _3^2\theta _4^2\right]}{\theta _2^2+\theta _3^2+\theta _4^2}}\right]^2`$ (100)
$`+`$ $`{\displaystyle \frac{\pi ^6\rho ^4}{756}}{\displaystyle \frac{\theta _3^8(\theta _2^2\theta _4^2)+\theta _4^8(\theta _2^2\theta _3^2)+\frac{5}{8}\left[\theta _3^2(\theta _2^8\theta _4^8)+\theta _4^2(\theta _2^8\theta _3^8)\right]+\frac{5}{16}(2\theta _2^{10}\theta _3^{10}\theta _4^{10})}{\theta _2^2+\theta _3^2+\theta _4^2}}`$ (101)
$`+`$ $`{\displaystyle \frac{31\pi ^5\rho ^3}{12096}}{\displaystyle \frac{\theta _3^{10}(\theta _2^4\theta _4^4)+\theta _2^{10}(\theta _3^4+\theta _4^4)\theta _4^{10}(\theta _3^4+\theta _2^4)}{\theta _2^2+\theta _3^2+\theta _4^2}}\left(1+4\rho {\displaystyle \frac{}{\rho }}\mathrm{ln}\theta _2\right)`$ (102)
$``$ $`{\displaystyle \frac{\pi ^5\rho ^3}{756}}{\displaystyle \frac{\theta _2^2\theta _3^2\theta _4^2\left[\theta _3^6(\theta _2^2\theta _4^2)+\theta _2^6(\theta _3^2+\theta _4^2)\theta _4^6(\theta _3^2+\theta _2^2)\right]}{\theta _2^2+\theta _3^2+\theta _4^2}}\left(1+4\rho {\displaystyle \frac{}{\rho }}\mathrm{ln}\theta _2\right)`$ (103)
$`\mathrm{}`$ (104)
3. Finally, the asymptotic expansion of the free energy, the internal energy and the specific heat of the Gaussian model after subtraction of the zero modes, $`\mathrm{ln}\mu \sqrt{8S}`$, $`\sqrt{2}S/\mu `$ and $`S(1/\mu +1/\mu ^2)`$, respectively, can be written as

$`f_{\mathrm{bulk}}`$ $`=`$ $`{\displaystyle \frac{2\gamma }{\pi }}\mathrm{ln}\sqrt{2}`$ (105)
$`f_0(\rho )`$ $`=`$ $`\mathrm{ln}\sqrt{\rho }\eta ^2`$ (106)
$`f_1(\rho )`$ $`=`$ $`{\displaystyle \frac{\pi ^3\rho ^2}{180}}(\theta _2^4\theta _4^4\theta _2^4\theta _3^4\theta _3^4\theta _4^4)`$ (107)
$`f_2(\rho )`$ $`=`$ $`{\displaystyle \frac{\pi ^6\rho ^4}{1512}}\theta _3^4\theta _4^4\left(\theta _2^82\theta _3^4\theta _4^4\right)+{\displaystyle \frac{\pi ^5\rho ^3}{1512}}(\theta _3^4+\theta _4^4)(\theta _3^4+\theta _2^4)(\theta _2^4\theta _4^4)\left(1+4\rho {\displaystyle \frac{}{\rho }}\mathrm{ln}\theta _2\right)`$ (108)
$`\mathrm{}`$ (109)
For internal energy there are no finite-size correction terms ($`u_{bulk}=0,u_p=0`$ for all $`p`$).

$`c_{\mathrm{bulk}}(\rho )`$ $`=`$ $`{\displaystyle \frac{8}{\pi }}\left(\mathrm{ln}\sqrt{{\displaystyle \frac{S}{\rho }}}+\mathrm{ln}{\displaystyle \frac{\sqrt{2}}{\pi }}+C_E{\displaystyle \frac{\pi }{4}}+\mathrm{ln}\eta \right)`$ (110)
$`c_0(\rho )`$ $`=`$ $`0`$ (111)
$`c_1(\rho )`$ $`=`$ $`{\displaystyle \frac{2}{3}}\left(1+4\rho {\displaystyle \frac{}{\rho }}\mathrm{ln}\eta +4\rho ^2{\displaystyle \frac{^2}{\rho ^2}}\mathrm{ln}\eta \right)`$ (112)
$`c_2(\rho )`$ $`=`$ $`0`$ (113)
$`\mathrm{}`$ (114)

## V Summary

In this paper, we have derived exact asymptotic expansion of the partition function with twisted boundary conditions at the critical point, Eqs. (33) and (35). As an application of this result, we have obtained exact asymptotic expansion of the free energy, the internal energy and the specific heat for a class of free exactly solvable models of statistical mechanics on torus. Moreover, the partition function of the dimer model on the Klein bottle can also be expressed via the partition function with twisted boundary conditions, namely $`Z_{\mathrm{Klein}}=Z_{\frac{1}{4},\frac{1}{2}}(0)`$ . Exact asymptotic expansion of the latter can immediately be written down with the help of our general formula (33). An interesting question is: whether this is also the case for other free models? This, however, is the problem for the future.

## A Euler-Maclaurin Summation Formula.

Suppose that $`F(x)`$ together with its first $`2m`$ derivatives is continuos within the interval $`(a,b)`$. Then the general Euler-Maclaurin summation formula states

$$\underset{n=0}{\overset{N1}{}}F(a+nh+\alpha h)=\frac{1}{h}_a^bF(\tau )d\tau +\underset{k=1}{\overset{p}{}}\frac{h^{k1}}{k!}\mathrm{B}_k(\alpha )\left(F^{(k1)}(b)F^{(k1)}(a)\right)\mathrm{R}_p(\alpha )$$
(A1)
where $`p2m`$, $`0\alpha 1`$, $`h=(ba)/N`$ and reminder term $`\mathrm{R}_p(\alpha )`$ is given by

$$\mathrm{R}_p(\alpha )=\frac{h^p}{p!}_0^1\widehat{\mathrm{B}}_p(\alpha \tau )\left\{\underset{n=0}{\overset{N1}{}}F^{(p)}(a+nh+\tau h)\right\}d\tau $$
(A2)
$`\widehat{\mathrm{B}}_p(\alpha )`$ are so called periodic Bernoulli functions which are defined as follows

$$\widehat{\mathrm{B}}_p(\alpha )=\frac{p!}{(2\pi i)^p}\underset{n0}{}\frac{e^{2\pi in\alpha }}{n^p}$$
(A3)
These functions have singularities at the integer values of $`\alpha `$ and on the interval $`\alpha [0,1]`$ they coincide with so-called Bernoulli polynomials

$`\mathrm{B}_1(\alpha )`$ $`=`$ $`\alpha \frac{1}{2}`$ (A4)
$`\mathrm{B}_2(\alpha )`$ $`=`$ $`\alpha ^2\alpha +\frac{1}{6}`$ (A5)
$`\mathrm{B}_3(\alpha )`$ $`=`$ $`\alpha ^3\frac{3}{2}\alpha ^2+\frac{1}{2}\alpha `$ (A6)
$`\mathrm{B}_4(\alpha )`$ $`=`$ $`\alpha ^42\alpha ^3+\alpha ^2\frac{1}{30}`$ (A7)
The generating function for the Bernoulli polynomials is

$$\frac{\lambda e^{\lambda \alpha }}{e^\lambda 1}=1+\underset{p=1}{\overset{\mathrm{}}{}}\frac{\lambda ^p}{p!}\mathrm{B}_p(\alpha )$$

Fourier transform of the generating function gives the important identity

$$\frac{e^{2\pi iz\alpha }}{e^{2\pi iz}1}=\underset{n=0}{\overset{\mathrm{}}{}}e^{2\pi iz(n+\alpha )}=\frac{1}{2\pi i}\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}\frac{e^{2\pi in\alpha }}{z+n}$$
(A8)
It is well known that Euler-Maclaurin summation formula is closely related to (generally speaking divergent) asymptotic series. For further discussion of the properties of these series the interested reader is referred to the book of Hardy .

In this paper we are mainly interested in sums of the form

$$\frac{1}{N}\underset{n=0}{\overset{N1}{}}f\left(\frac{\pi (n+\alpha )}{N}\right)$$
(A9)
The asymptotic expansion of the sum (A9) in the limit $`N\mathrm{}`$ can be obtained from Eq. (A1) by setting $`(a=0,b=\pi )`$. If we assume that all the derivatives of $`f(x)`$ are integrable over the interval $`(0,\pi )`$, i.e., the integral in Eq. (A2) is finite, we can formally extend the sum in Eq. (A1) to $`k=\mathrm{}`$ and drop the reminder term $`\mathrm{R}_p(\alpha )`$. In this case we can write the asymptotic expansion of the sum (A9) as follows

$$\frac{1}{N}\underset{n=0}{\overset{N1}{}}f\left(\frac{\pi (n+\alpha )}{N}\right)=\frac{1}{\pi }_0^\pi f(\tau )d\tau +\frac{1}{\pi }\underset{k=1}{\overset{\mathrm{}}{}}\left(\frac{\pi }{N}\right)^k\frac{\mathrm{B}_k(\alpha )}{k!}\left(f^{(k1)}(\pi )f^{(k1)}(0)\right)$$
(A10)
Note, that we also use abbreviated notation for Bernoulli polynomials: $`\mathrm{B}_p(\alpha )=\mathrm{B}_p^\alpha `$.

## B Relation between moments and cumulants

Moments $`Z_k`$ and cumulants $`F_k`$ which enters the expansion of exponent

$$\mathrm{exp}\left\{\underset{k=1}{\overset{\mathrm{}}{}}\frac{x^k}{k!}F_k\right\}=1+\underset{k=1}{\overset{\mathrm{}}{}}\frac{x^k}{k!}Z_k$$

are related to each other as

$`Z_1`$ $`=`$ $`F_1`$ (B1)
$`Z_2`$ $`=`$ $`F_2+F_1^2`$ (B2)
$`Z_3`$ $`=`$ $`F_3+3F_1F_2+F_1^3`$ (B3)
$`Z_4`$ $`=`$ $`F_4+4F_1F_3+3F_2^2+6F_1^2F_2+F_1^4`$ (B4)
$`\mathrm{}`$ (B5)
$`Z_k`$ $`=`$ $`{\displaystyle \underset{r=1}{\overset{k}{}}}{\displaystyle \left(\frac{F_{k_1}}{k_1!}\right)^{i_1}\mathrm{}\left(\frac{F_{k_r}}{k_r!}\right)^{i_r}\frac{k!}{i_1!\mathrm{}i_r!}}`$ (B6)
where summation is over all positive numbers $`\{i_1\mathrm{}i_r\}`$ and different positive numbers $`\{k_1,\mathrm{},k_r\}`$ such that $`k_1i_1+\mathrm{}+k_ri_r=k`$.

## C Elliptic Theta Functions

We adopt the following definition of the elliptic $`\theta `$-functions:

$`\theta _{\alpha ,\beta }(z,\tau )`$ $`=`$ $`{\displaystyle \underset{nZ}{}}\mathrm{exp}\left\{\pi i\tau \left(n+\frac{1}{2}\alpha \right)^2+2\pi i\left(n+\frac{1}{2}\alpha \right)\left(z+\frac{1}{2}\beta \right)\right\}`$ (C1)
$`=`$ $`\eta (\tau )\mathrm{exp}\left\{\pi i\tau \left(\alpha ^2\alpha +\frac{1}{6}\right)+2\pi i\left(\frac{1}{2}\alpha \right)\left(z+\frac{1}{2}\beta \right)\right\}`$ (C2)
$`\times `$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left[\mathrm{\hspace{0.17em}1}e^{2\pi i\tau \left(n+\alpha \right)2\pi i\left(z\beta \right)}\right]\left[\mathrm{\hspace{0.17em}1}e^{2\pi i\tau \left(n+1\alpha \right)+2\pi i\left(z\beta \right)}\right]`$ (C3)
These should be compared with the notations of Mumford .

The elliptic $`\theta `$-functions satisfies the heat equation

$$\frac{}{\tau }\theta _{\alpha ,\beta }(z,\tau )=\frac{1}{4\pi i}\frac{^2}{z^2}\theta _{\alpha ,\beta }(z,\tau )$$
(C4)
Relation to standard notations is

$`\theta _{0,0}(z,\tau )`$ $`=`$ $`\theta _1(z,\tau )`$ (C5)
$`\theta _{0,\frac{1}{2}}(z,\tau )`$ $`=`$ $`\theta _2(z,\tau )`$ (C6)
$`\theta _{\frac{1}{2},0}(z,\tau )`$ $`=`$ $`\theta _4(z,\tau )`$ (C7)
$`\theta _{\frac{1}{2},\frac{1}{2}}(z,\tau )`$ $`=`$ $`\theta _3(z,\tau )`$ (C8)
The Dedekind $`\eta `$-function is usually defined as

$$\eta (\tau )=e^{\pi i\tau /12}\underset{n=1}{\overset{\mathrm{}}{}}\left[\mathrm{\hspace{0.17em}1}e^{2\pi i\tau n}\right]$$

Considering the functions $`\theta _{\alpha ,\beta }(\tau )=\theta _{\alpha ,\beta }(0,\tau )`$ and $`\eta (\tau )`$ of pure imaginary aspect ratio, $`\tau =i\rho `$, we obtain the identity

$`\mathrm{ln}\left|{\displaystyle \frac{\theta _{\alpha ,\beta }(i\rho )}{\eta (i\rho )}}\right|+\pi \rho \mathrm{B}_2^\alpha `$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{ln}\left|\mathrm{\hspace{0.17em}1}e^{2\pi \rho \left(n+\alpha \right)2\pi i\beta }\right|`$ (C10)
$`+{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{ln}\left|\mathrm{\hspace{0.17em}1}e^{2\pi \rho \left(n+1\alpha \right)2\pi i\beta }\right|`$
Taking the derivative of Eq. (C10) with respect to $`\rho `$ we can obtain the following useful identity

$`\mathrm{𝚁𝚎}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\left[(n+\alpha )e^{2\pi m\left[\rho \left(n+\alpha \right)+i\beta \right]}+(n+1\alpha )e^{2\pi m\left[\rho \left(n+1\alpha \right)+i\beta \right]}\right]`$ (C11)
$`={\displaystyle \frac{\mathrm{B}_2^\alpha }{2}}+{\displaystyle \frac{1}{2\pi }}{\displaystyle \frac{}{\rho }}\mathrm{ln}\left|{\displaystyle \frac{\theta _{\alpha ,\beta }(i\rho )}{\eta (i\rho )}}\right|`$ (C12)
Taking the second derivative of $`\mathrm{ln}\left|\theta _{\alpha ,\beta }(z,i\rho )\right|`$ with respect to $`z`$ at $`z=0`$ and using the heat equation Eq. (C4) we obtain

$`{\displaystyle \frac{}{\rho }}\left(\mathrm{𝚁𝚎}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\left\{{\displaystyle \frac{1}{n+\alpha }}e^{2\pi m\left[\rho \left(n+\alpha \right)+i\beta \right]}+{\displaystyle \frac{1}{n+1\alpha }}e^{2\pi m\left[\rho \left(n+1\alpha \right)+i\beta \right]}\right\}\right)=`$ (C13)
$`=2{\displaystyle \frac{}{\rho }}\mathrm{ln}\left|\theta _{\alpha ,\beta }(i\rho )\right|+{\displaystyle \frac{1}{2\pi }}\left({\displaystyle \frac{}{z}}\mathrm{ln}\left|\theta _{\alpha ,\beta }(i\rho )\right|\right)^2`$ (C14)
## D Kronecker’s Double Series

Kronecker’s double series can be defined as

$$\mathrm{K}_p^{\alpha ,\beta }(\tau )=\frac{p!}{(2\pi i)^p}\underset{\genfrac{}{}{0.0pt}{}{m,nZ}{(m,n)(0,0)}}{}\frac{e^{2\pi i(n\alpha +m\beta )}}{(n+\tau m)^p}$$

In this form, however, they cannot be directly applied to our analysis. We need to cast them in a different form. To this end, let us separate from the double series a subseries with $`m=0`$

$$\mathrm{K}_p^{\alpha ,\beta }(\tau )=\frac{p!}{(2\pi i)^p}\underset{n0}{}\frac{e^{2\pi in\alpha }}{n^p}\frac{p!}{(2\pi i)^p}\underset{m0}{}\underset{nZ}{}\frac{e^{2\pi i(n\alpha +m\beta )}}{(n+\tau m)^p}$$

Here the first sum gives nothing but Fourier representation of Bernoulli polynomials Eq.(A3). The second sum can be rearranged with the help of the identity

$$\frac{p!}{(2\pi i)^p}\underset{nZ}{}\frac{e^{2\pi in\alpha }}{(z+n)^p}=p\underset{n=0}{\overset{\mathrm{}}{}}(n+\alpha )^{p1}e^{2\pi iz(n+\alpha )}$$

which can easily be derived from Eq.(A8) differentiating it $`p`$ times. The final result of our resummation of double Kronecker sum is

$$\mathrm{K}_p^{\alpha ,\beta }(\tau )=\mathrm{B}_p^\alpha p\underset{m0}{}\underset{n=0}{\overset{\mathrm{}}{}}(n+\alpha )^{p1}e^{2\pi im(\tau (n+\alpha )\beta )}$$

Considering real part of the Kronecker sums with pure imaginary aspect ratio, $`\tau =i\rho `$, we can further rearrange this expression to get summation only over positive $`m1`$

$`\mathrm{B}_{2p}^\alpha \mathrm{𝚁𝚎}\mathrm{K}_{2p}^{\alpha ,\beta }(i\rho )`$ $`=`$ $`{\displaystyle \frac{(2p)!}{(4\pi ^2)^p}}\mathrm{𝚁𝚎}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{nZ}{}}{\displaystyle \frac{e^{2\pi i(n\alpha +m\beta )}+e^{2\pi i(n(1\alpha )m\beta )}}{(n+i\rho m)^{2p}}}`$ (D1)
$`=`$ $`2p\mathrm{𝚁𝚎}\{{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}(n+\alpha )^{2p1}e^{2\pi m(\rho (n+\alpha )+i\beta )}`$ (D3)
$`+{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}(n+1\alpha )^{2p1}e^{2\pi m(\rho (n+1\alpha )+i\beta )}\}`$
## E Asymptotic Expansion of the Digamma Function

Let us start with well known expansion of the digamma function $`\psi (N)`$

$`\psi (x)`$ $`=`$ $`\mathrm{ln}x{\displaystyle \frac{1}{2x}}{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{B}_{2p}}{2p}}{\displaystyle \frac{1}{x^{2p}}}`$ (E1)
$`=`$ $`\mathrm{ln}x{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}(1)^p{\displaystyle \frac{\mathrm{B}_p}{p}}{\displaystyle \frac{1}{x^p}}`$ (E2)
Plugging in the above expansion $`x=N+\alpha `$ and expand the resulting factors $`\mathrm{ln}(1+\alpha /N)`$, $`(1+\alpha /N)^p`$ in powers of $`N^1`$ we obtain

$`\psi (N+\alpha )`$ $`=`$ $`\mathrm{ln}N{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}(1)^p{\displaystyle \frac{\alpha ^p}{pN^p}}{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}(1)^{k+p}\mathrm{B}_p{\displaystyle \frac{(p+k1)!}{k!p!}}{\displaystyle \frac{\alpha ^k}{N^{p+k}}}`$ (E3)
$`=`$ $`\mathrm{ln}N{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}(1)^p{\displaystyle \frac{\alpha ^p}{pN^p}}{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{p=1}{\overset{l}{}}}(1)^l\mathrm{B}_p{\displaystyle \frac{(l1)!}{(lp)!p!}}{\displaystyle \frac{\alpha ^{lp}}{N^l}}`$ (E4)
$`=`$ $`\mathrm{ln}N{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{p=0}{\overset{l}{}}}(1)^l\mathrm{B}_p{\displaystyle \frac{(l1)!}{(lp)!p!}}{\displaystyle \frac{\alpha ^{lp}}{N^l}}`$ (E5)
Using the relation between Bernoulli polynomials $`\mathrm{B}_p^\alpha `$ and Bernoulli numbers $`\mathrm{B}_p`$

$$\mathrm{B}_l^\alpha =\underset{p=0}{\overset{l}{}}\mathrm{B}_p\frac{l!}{(lp)!p!}\alpha ^{lp}$$
(E6)
we finally obtain the Eq. (44)

$$\psi (N+\alpha )=\mathrm{ln}N\underset{p=1}{\overset{\mathrm{}}{}}(1)^p\frac{\mathrm{B}_p^\alpha }{p}\frac{1}{N^p}$$
(E7)
## F Reduction of Kronecker’s Double Series to Theta Functions

Let us consider Laurent expansion of the Weierstrass function

$`\mathrm{}(z)`$ $`=`$ $`{\displaystyle \frac{1}{z^2}}+{\displaystyle \underset{(n,m)(0,0)}{}}\left[{\displaystyle \frac{1}{(zn\tau m)^2}}{\displaystyle \frac{1}{(n+\tau m)^2}}\right]`$ (F1)
$`=`$ $`{\displaystyle \frac{1}{z^2}}+{\displaystyle \underset{p=2}{\overset{\mathrm{}}{}}}a_p(\tau )z^{2p2}`$ (F2)
The coefficients $`a_p(\tau )`$ of the expansion can all be written in terms of the elliptic $`\theta `$-functions with the help of the recursion relation

$$a_p=\frac{3}{(p3)(2p+1)}(a_2a_{p2}+a_3a_{p3}+\mathrm{}+a_{p2}a_2)$$

where first terms of the sequence are

$`a_2`$ $`=`$ $`\frac{\pi ^4}{15}(\theta _2^4\theta _3^4\theta _2^4\theta _4^4+\theta _3^4\theta _4^4)`$ (F3)
$`a_3`$ $`=`$ $`\frac{\pi ^6}{189}(\theta _2^4+\theta _3^4)(\theta _4^4\theta _2^4)(\theta _3^4+\theta _4^4)`$ (F4)
$`a_4`$ $`=`$ $`\frac{1}{3}a_2^2`$ (F5)
$`a_5`$ $`=`$ $`\frac{3}{11}(a_2a_3)`$ (F6)
$`a_6`$ $`=`$ $`\frac{1}{39}(2a_2^3+3a_3^2)`$ (F7)
$`\mathrm{}`$ (F8)
Kronecker functions $`\mathrm{K}_{2p}^{0,0}(\tau )`$ are related directly to the coefficients $`a_p(\tau )`$

$$\mathrm{K}_{2p}^{0,0}(\tau )=\frac{(2p)!}{(4\pi ^2)^p}\frac{a_p(\tau )}{(2p1)}$$

Kronecker functions $`\mathrm{K}_{2p}^{\alpha ,\beta }(\tau )`$ with $`\alpha `$ and $`\beta `$ taking values $`0`$ or $`1/2`$ can in their turn be related to the function $`\mathrm{K}_{2p}^{0,0}(\tau )`$ by means of simple resummation of Kronecker’s double series

$`\mathrm{K}_p^{0,\frac{1}{2}}(\tau )`$ $`=`$ $`2\mathrm{K}_p^{0,0}(2\tau )\mathrm{K}_p^{0,0}(\tau )`$ (F9)
$`\mathrm{K}_p^{\frac{1}{2},0}(\tau )`$ $`=`$ $`2^{1p}\mathrm{K}_p^{0,0}(\tau /2)\mathrm{K}_p^{0,0}(\tau )`$ (F10)
$`\mathrm{K}_p^{\frac{1}{2},\frac{1}{2}}(\tau )`$ $`=`$ $`(1+2^{2p})\mathrm{K}_p^{0,0}(\tau )2^{1p}\mathrm{K}_p^{0,0}(\tau /2)2\mathrm{K}_p^{0,0}(2\tau )`$ (F11)
Thus, Kronecker functions $`\mathrm{K}_{2p}^{\alpha ,\beta }(\tau )`$ with $`\alpha `$ and $`\beta `$ taking values $`0`$ or $`1/2`$ can all be expressed in terms of the elliptic $`\theta `$-functions only. For practical calculations the following identities are also helpful

| $`2\theta _2^2(2\tau )`$ | $`=`$ | $`\theta _3^2\theta _4^2`$ |  | $`\theta _2^2(\tau /2)`$ | $`=`$ | $`2\theta _2\theta _3`$ |
| --- | --- | --- | --- | --- | --- | --- |
| $`2\theta _3^2(2\tau )`$ | $`=`$ | $`\theta _3^2+\theta _4^2`$ |  | $`\theta _3^2(\tau /2)`$ | $`=`$ | $`\theta _2^2+\theta _3^2`$ |
| $`2\theta _4^2(2\tau )`$ | $`=`$ | $`2\theta _3\theta _4`$ |  | $`\theta _4^2(\tau /2)`$ | $`=`$ | $`\theta _3^2\theta _2^2`$ |

From the general formulas above we can easily write down all the Kronecker functions that have appeared in our asymptotic expansions

$`\mathrm{K}_4^{0,0}(\tau )`$ $`=`$ $`\frac{1}{30}(\theta _2^4\theta _4^4\theta _2^4\theta _3^4\theta _3^4\theta _4^4)`$ (F12)
$`\mathrm{K}_4^{0,\frac{1}{2}}(\tau )`$ $`=`$ $`\frac{1}{30}(\frac{7}{8}\theta _2^8\theta _3^4\theta _4^4)`$ (F13)
$`\mathrm{K}_4^{\frac{1}{2},0}(\tau )`$ $`=`$ $`\frac{1}{30}(\frac{7}{8}\theta _4^8\theta _2^4\theta _3^4)`$ (F14)
$`\mathrm{K}_4^{\frac{1}{2},\frac{1}{2}}(\tau )`$ $`=`$ $`\frac{1}{30}(\frac{7}{8}\theta _3^8+\theta _2^4\theta _4^4)`$ (F15)
$`\mathrm{K}_6^{0,0}(\tau )`$ $`=`$ $`\frac{1}{84}(\theta _2^4+\theta _3^4)(\theta _4^4\theta _2^4)(\theta _3^4+\theta _4^4)`$ (F16)
$`\mathrm{K}_6^{0,\frac{1}{2}}(\tau )`$ $`=`$ $`\frac{1}{84}(\theta _3^4+\theta _4^4)(\frac{31}{16}\theta _2^8+\theta _3^4\theta _4^4)`$ (F17)
$`\mathrm{K}_6^{\frac{1}{2},0}(\tau )`$ $`=`$ $`\frac{1}{84}(\theta _2^4+\theta _3^4)(\frac{31}{16}\theta _4^8+\theta _2^4\theta _3^4)`$ (F18)
$`\mathrm{K}_6^{\frac{1}{2},\frac{1}{2}}(\tau )`$ $`=`$ $`\frac{1}{84}(\theta _2^4\theta _4^4)(\frac{31}{16}\theta _3^8\theta _4^4\theta _2^4)`$ (F19)
$`\mathrm{}`$ (F20)
Note that when $`\rho \mathrm{}`$ we have limits $`\theta _20`$, $`\theta _41`$, $`\theta _31`$ and the Kronecker’s function reduce to the Bernoulli polynomials.

We are grateful to V.B. Priezzhev for critical remarks on the manuscript. This work was supported by the Russian Foundation for Basic Research through grant No. 99-01-00882, by the Swiss National Science Foundation through grant No. 7-SUPJ-62295 and by the National Science Counsil of the Republic of China (Taiwan) through grant No. NSC89-2112-M-001-084.