text:
# 1 Introduction

## 1 Introduction

Several years ago, the SO(5) model was proposed as a description of high-temperature superconductors . At low temperatures these materials exhibit superconductivity (SC) or antiferromagnetism (AF) depending on the doping; the SO(5) model attempts to unify these two phenomena (both of which involve spontaneous symmetry breaking) into a single symmetry group.

Not surprisingly, the SO(5) model is somewhat more complicated than the corresponding model describing conventional superconductivity (which is based on the smaller group SO(2)); correspondingly, richer behavior can be seen. In conventional superconductivity, one parameter, the Ginzburg-Landau (GL) parameter $`\kappa `$, determines the behavior. For instance, $`\kappa =1/\sqrt{2}\kappa _c`$ marks the boundary between type I and type II superconductors.

In the SO(5) model, there are two dimensionless parameters, which can be taken to be $`\kappa `$ and a second parameter, $`\beta `$ (to be defined below); the latter is related to the degree to which SO(5) is explicitly broken (by doping, for instance). As will be shown below, this second parameter has a profound influence on the magnetic behavior of an SO(5) superconductor; for instance, the critical value of $`\kappa `$ in the SO(5) model is given by $`\kappa _c(\beta )=(1/\sqrt{2})\sqrt{(1+\beta )/(1\beta )}`$. One sees that as $`\beta 1^{}`$, $`\kappa _c\mathrm{}`$, which is a dramatic departure from the conventional value. This is of some significance since high-temperature superconductors typically have $`\kappa 1`$, and therefore are normally thought to be extreme type II superconductors. Such a conclusion would be premature in the SO(5) model, however, since (as the above relation shows) one could have $`\kappa _c>\kappa 1`$.

In this paper we will first remind the reader of the magnetic properties of conventional superconductors. The standard approach based on the surface energy density of a boundary between a normal phase at critical magnetic field $`H_c`$ and a superconducting phase will be briefly reviewed. We will then discuss an alternative approach based on vortex energetics. Although much of the above is fairly familiar material, this discussion will establish notation and set the stage for the parallel discussion in the context of the SO(5) model.

Next, we will discuss the case of SO(5) superconductivity. After a brief review of the model itself, the critical fields $`H_c`$ and $`H_{c2}`$ will be calculated. From these, the critical value of the parameter $`\kappa `$ can be calculated explicitly as a function of $`\beta `$. This expression can be confirmed numerically, either by examining the surface energy density or vortex energetics. In the latter approach, one can see that the dramatically different behavior in the SO(5) model is due to the possibility of an antiferromagnetic core in the vortex .

## 2 Review of conventional superconductivity111For a more complete discussion on conventional SC, see  or

A conventional superconductor is described by the following Helmholtz free energy in the GL theory:

$$\widehat{F}=𝑑𝐱\left\{f_n\frac{a_1^2}{2}|\varphi |^2+\frac{b^2}{4}|\varphi |^4+\frac{1}{2m^{}}\left|\left(i\mathrm{}\frac{e^{}𝐀}{c}\right)\varphi \right|+\frac{𝐡^2}{8\pi }\right\},$$
(1)
where $`f_n`$ is a constant, $`\widehat{𝐡}=\times 𝐀`$ is the microscopic magnetic field and $`a_1`$, $`b`$ are parameters. The minimum of the potential is $`|\varphi |^2=a_1^2/bv^2`$.

The Helmholtz free energy is minimized when the temperature $`T`$ and interior magnetic induction $`B`$ are constant (i.e., thermodynamic equilibrium). One uses the Gibbs free energy for equilibrium where $`T`$ and the external magnetic field $`H`$ are constant. These two quantities are related by the following Legendre transformation:

$`g(H,T)=f(h,T){\displaystyle \frac{1}{4\pi }}hH.`$ (2)
The second term in this equation is the microscopic version of $`(4\pi )^1BH`$, since the magnetic induction $`B(𝐫)=h(𝐫)`$.

The critical magnetic field $`H_c`$ is defined as the field which establishes the following condition on the Gibbs energies in the normal and superconductor regions:

$`g_n(H_c,T)=g_s(H_c,T).`$ (3)
We also have the following relations, the first one due to the Meissner effect:

$`g_s(H,T)=f_s(T),`$
$`g_n(H,T)=f_n{\displaystyle \frac{H^2}{4\pi }}.`$
From Eq. (1), we can obtain a correspondence with the parameters of the theory:

$`f_n(T)f_s(T)={\displaystyle \frac{a_1^4}{4b}}.`$
One finds:

$`{\displaystyle \frac{H_c^2}{4\pi }}={\displaystyle \frac{a_1^4}{2b}}={\displaystyle \frac{\mathrm{\Phi }_0^2}{32\pi ^3\lambda ^2\xi ^2}},`$ (4)
where $`\mathrm{\Phi }_0=hc/e^{}`$ is the quantum of magnetic flux, $`\lambda =(m^{}c^2/4\pi e^^2v^2)^{1/2}`$ is the penetration depth and $`\xi =(\mathrm{}^2/m^{}a_1^2)^{1/2}`$ is the coherence length.

There is also another critical field, $`H_{c2}`$, which is the strongest field allowed in a superconducting region before the system goes to the normal state. It is given by:

$`H_{c2}={\displaystyle \frac{\mathrm{\Phi }_0}{2\pi \xi ^2}}=\sqrt{2}\kappa H_c,`$
where, as usual, $`\kappa =\lambda /\xi `$. We see that if $`\kappa `$ is equal to its critical value $`1/\sqrt{2}`$, $`H_c=H_{c2}`$ (more on this in Sections 2.1 and 3.2).

### 2.1 Surface Energy

As is well known, a superconductor will behave in one of two ways when placed in a magnetic field. In a type I SC, macroscopic normal regions where $`h=H_c`$ will form, while in a type II SC, a lattice of vortices of flux $`\mathrm{\Phi }_0`$ will appear. One way to determine under which circumstances a SC is of type I or II is by studying the energy per unit area of a boundary between a normal region and a SC region. If this energy is positive, the system will tend to minimize the interface’s surface area in order to lower this energy as much as possible, indicative of a type I SC. Similarly for negative surface energy, this surface area will be maximized, inducing type II behavior.

Analytically, the surface energy is given by:

$`\widehat{G}_s={\displaystyle \frac{\mathrm{\Phi }_0^2\kappa ^2}{32\pi ^3\lambda ^3}}I={\displaystyle \frac{H_c^2\lambda }{4\pi }}I,`$
where:

$`I={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑s\left\{{\displaystyle \frac{1}{2}}(1f^4)+h^2\sqrt{2}h\right\}.`$ (5)
Consider therefore a field configuration between the normal state for $`x\mathrm{}`$ and the SC state for $`x+\mathrm{}`$:

$$\varphi =\varphi (x),𝐀=A(x)\widehat{y}.$$

Let us go to a description in terms of dimensionless variables, setting:

$$x=\lambda s,\varphi =vf,A=\frac{c\mathrm{}}{e^{}\xi }a,\widehat{h}=\frac{c\mathrm{}}{e^{}\xi \lambda }h.$$
(6)
We obtain the following dimensionless free energy:

$$F=𝑑s\left\{f_n+\frac{1}{\kappa ^2}f^2+(a^2+f^21)f^2+h^2\right\}.$$
(7)
Minimizing F yields the following equations (the so-called 1D GL equations):

$$\frac{1}{\kappa ^2}\frac{d^2f}{ds^2}+(1a^2)ff^3=0,$$
(8)
$$\frac{d^2a}{ds^2}af^2=0.$$
(9)
The sign of the surface energy can be determined qualitatively in the limit of either small or large value of $`\kappa `$ in Eq. (8), by looking at the length scales of $`h`$ and $`f`$. By construction (see Eq. (6)), the length scale of $`h`$ is $`l_h=1`$, while that of $`f`$ turns out to be $`l_f=1/(\sqrt{2}\kappa ))`$.

From Fig. 1, we can see how the sign of the surface energy is determined by the behavior of $`f`$ and $`h`$ in those limiting cases. In the first limiting case (Fig. 1a) which corresponds to a type I sample ($`\kappa 1`$, so that $`l_fl_h`$), the energy comes from the expulsion of the magnetic field within the SC region. More explicitly, $`h`$ is rapidly expulsed at the interface while $`f`$ slowly assumes its asymptotic value. This gives rise to a positive scalar field energy density without a cancelling negative magnetic energy; the surface energy is positive. Fig. 1b shows the reverse case of a type II sample, where $`f`$ quickly reaches its asymptotic behavior while the magnetic field slowly decreases, resulting in a negative surface energy.

The critical value of $`\kappa `$ is the one for which the surface energy (or equivalently, the integral (5)) is zero. The well-known result $`\kappa _c=1/\sqrt{2}`$ can be determined either numerically 222 A numerical calculation of the surface energy is given in ., or by making the oracular observation that when $`f`$ and $`h`$ obey the equations of motion for $`\kappa =1/\sqrt{2}`$, then the integrand of (5) is zero333Note that there is no need for the integrand of (5) itself to be zero; the weaker condition that the integral is zero would suffice. Indeed, we will see in the next section that in the SO(5) model the analogous integrand is *not* zero when the integrated surface energy is zero.. This observation can be explained by a recasting of the work of Bogomol’nyi  (to be discussed in the next section) into a form appropriate to the surface energy.

### 2.2 Vortex Energetics

An alternative approach to determining the behavior of a superconductor when placed in a magnetic field is the study of the energy of vortices as a function of their winding number $`m`$. In this section, we will see what information can be drawn from the analytic expression, as well as an important limiting case, $`m\mathrm{}`$. In this limit, we will also show how one can establish a link between the vortex and surface energies.

As mentioned above, when placed in a magnetic field, the system will eventually find the configuration which minimizes the free energy. We can construct a function $`(m)=F(m)/m`$ which is the free energy of a vortex of winding number $`m`$ divided by the winding number, or in other words the free energy per flux quantum of an $`m`$-vortex. For large $`m`$, if $`(m)`$ is an increasing function of $`m`$, it is more energetically expensive for a magnetic field to penetrate in vortices of large winding number; rather, a lattice of vortices of winding number 1 will form. This, clearly, is what we expect of a type II superconductor. Similarly, if $`(m)`$ is a decreasing function of $`m`$, the least energetic configuration will be a macroscopic normal region containing a large number of flux quanta; this is what we expect of a type I superconductor.

For the $`m`$-vortex configuration, we use an ansatz that is generalized to winding number $`m`$:

$$\varphi =f(r)e^{im\theta },\widehat{A_i}=\frac{a_1c\sqrt{m^{}}}{e^{}}ϵ_{ij}\frac{r_j}{r}a(r).$$

Again we can use the dimensionless variables defined in (6). Here, however, $`s`$ denotes the radial variable in a 2D plane, as opposed to the cartesian coordinate orthogonal to the domain wall.

Using these changes in Eq. (1) yields the vortex energy:

$`F_v={\displaystyle _0^{\mathrm{}}}ds{\displaystyle \frac{s}{2}}\{{\displaystyle \frac{1}{2}}+(a^{}+{\displaystyle \frac{a}{s}})^2+{\displaystyle \frac{1}{\kappa ^2}}(f^2+({\displaystyle \frac{m}{s}}+\kappa a)^2f^2)`$
$`(1f^2)f^2\},`$ (10)
as well as the equations of motion:

$$\frac{1}{\kappa ^2}\left(f^{\prime \prime }+\frac{1}{s}f^{}\left(\frac{m}{s}+\kappa a\right)^2f\right)+ff^3=0,$$
(11)
$$h^{}+\left(\frac{m}{\kappa s}+a\right)f^2=0.$$
(12)
It is possible to rewrite Eq. (2.2) in the following form :

$`2\kappa ^2F(m)={\displaystyle _0^{\mathrm{}}}sds\{[f^{}({\displaystyle \frac{m}{s}}+\kappa a)f]^2`$
$`+[\kappa (a^{}+{\displaystyle \frac{a}{s}})+{\displaystyle \frac{1}{2}}(1f^2)]^2+{\displaystyle \frac{1}{2}}(\kappa ^2{\displaystyle \frac{1}{2}})(1f^2)^2+{\displaystyle \frac{d\chi }{ds}}+{\displaystyle \frac{1}{s}}\chi \},`$ (13)
where $`\chi =\kappa a(f^21)+mf^2/s`$. Integrating the last two terms yields the winding number $`m`$, independent of the details of the field. Clearly, if $`\kappa =1/\sqrt{2}`$, the third term vanishes and $`F`$ consists of the winding number $`m`$ plus two positive semidefinite terms; it will be minimized if they are zero. One can see that these terms are in fact zero if $`f`$ and $`a`$ satisfy the equations of motion. In this case, $`(m)=F(m)/m=1`$ independent of $`m`$, which indeed corresponds to critical behavior with regard to vortice lattice stability, as there is no preference for either type I (coalescing) or type II (lattice) behavior.

We can determine the behavior of $`(m)`$ numerically for a large value of the winding number $`m`$. Our numerical results already suggest that $`(m)`$ reaches asymptotically a constant value at $`m\mathrm{}`$. In this limiting case, the core of the vortex is very wide, and the fields fluctuate rapidly at only one spatial position– that is at the vortex perimeter. Let us define this position as $`s_m`$. Analytically, we will see this more clearly by inserting the limiting forms of the field into Eq. (2.2), which are:

$`s<s_m:f0,a{\displaystyle \frac{H_c}{2}}s,h=0,`$
$`s>s_m:f1,a{\displaystyle \frac{m}{\kappa s}},h=a^{}{\displaystyle \frac{a}{s}}{\displaystyle \frac{1}{\sqrt{2}}}.`$
After integration, the second and third terms of Eq. (2.2) become, respectively:

$$\left(\frac{1}{2}\kappa H_c\right)^2\frac{s_m^2}{2}\text{and}\frac{1}{4}(\kappa ^2\frac{1}{2})s_m^2.$$

We have $`m`$ and the remaining bracketed term, which does not contribute to the integral (this, of course, stems from our choice of step-functions as the role of fields $`f`$ and $`h`$, since we essentially remove the presence of the surface altogether). Let us rewrite the variable $`s_m`$ in terms of the quantum of magnetic flux $`\mathrm{\Phi }_0`$, using our $`m\mathrm{}`$ ansätze:

$`\mathrm{\Phi }_0={\displaystyle \frac{1}{m}}{\displaystyle _0^{\mathrm{}}}hd^2x={\displaystyle \frac{\pi s_m^2H_c}{m}}.`$
Therefore:

$`s_m^2={\displaystyle \frac{2m}{\kappa H_c}}.`$
We find:

$`(m)={\displaystyle \frac{1}{\sqrt{2}\kappa }}(m\mathrm{}).`$
If we add thickness to the interface, the first term of Eq. (2.2), which was vanishing, will instead be proportional to $`\sqrt{m}`$ (i.e., a “perimeter” term $`s_m`$).

We are now able to establish a correspondence between the free energies of the vortex and wall domain configurations. Indeed, for a large value of $`m`$, the perimeter of the vortex is large enough to significantly lower the curvature of the vortex locally, likening a given small region into the form of a domain wall. We should expect to see confirmation of this, using the appropriate approximations. We will do this here.

Let us start from the Gibbs energy of the vortex and derive an expression that clearly shows the relation with the domain wall Gibbs energy. We have:

$$\widehat{G}_v=\widehat{F}_v\frac{H_c}{4\pi }\mathrm{\Phi }_0m=\frac{\mathrm{\Phi }_0^2}{16\pi ^2\lambda ^2}\left(2\kappa ^2F(m)\sqrt{2}\kappa m\right).$$

Substituting $`F(m)`$ with Eq. (2.2) and expanding the terms in square parentheses:

$`\widehat{G}_v={\displaystyle \frac{\mathrm{\Phi }_0^2}{16\pi ^2\lambda ^2}}(m+{\displaystyle _0^{\mathrm{}}}sds\{{\displaystyle \frac{\kappa ^2}{2}}(1f^2)^2+\kappa ^2h^2\kappa h(1f^2)`$
$`+f^22ff^{}({\displaystyle \frac{m}{s}}+\kappa a)+({\displaystyle \frac{m}{s}}+\kappa a)^2f^2\}\sqrt{2}\kappa m).`$ (14)
We notice that the following terms are equal to a familiar quantity:

$`\kappa h(1f^2)2ff^{}\left({\displaystyle \frac{m}{s}}+\kappa a\right)={\displaystyle \frac{d\chi }{ds}},`$
which, once integrated, cancels the $`m`$ term. If we take Eq. (11) and multiply it by $`f\kappa ^2`$, we have:

$`\left({\displaystyle \frac{m}{s}}+\kappa a\right)^2f^2=ff^{\prime \prime }+\kappa ^2(f^2f^4).`$ (15)
Next, we change variables $`sz=ss_m`$, where it is understood that the region of interest is $`|z||s_m|`$. Therefore

$`{\displaystyle _0^{\mathrm{}}}s𝑑s{\displaystyle _{\mathrm{}}^{\mathrm{}}}(z+s_m)𝑑z{\displaystyle _{\mathrm{}}^{\mathrm{}}}s_m𝑑z.`$
Substituting (15) into (14) and simplifying, we have:

$`\widehat{G}_v={\displaystyle \frac{\kappa ^2\mathrm{\Phi }_0^2}{16\pi ^2\lambda ^2}}({\displaystyle _{\mathrm{}}^{\mathrm{}}}s_mdz\{{\displaystyle \frac{1}{2}}(1f^4)+h^2\sqrt{2}h\}`$
$`+{\displaystyle \frac{s_m}{\kappa ^2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d(ff^{})}{dz}}dz).`$ (16)
We recognize the first integral: it is $`s_m`$ times $`I`$, as defined by Eq. (5); the second integral is zero due to the boundary conditions. Therefore, the link between the Gibbs energies between the vortex and domain wall case is summed up by the following relation, valid for large $`m`$:

$`\widehat{G}_v={\displaystyle \frac{\kappa ^2\mathrm{\Phi }_0^2s_m}{16\pi ^2\lambda ^2}}I=\widehat{G}_s(2\pi \lambda s_m).`$
In this section, we have studied vortex energetics as an alternative approach to the determination of the type (I vs. II) of a given superconductor, in the conventional case. We now turn our attention to the case of SO(5) superconductivity.

## 3 SO(5) Superconductivity

Having briefly reviewed the basics of conventional SC, let us turn to the case of SO(5) superconductivity. As mentioned in the introduction, the SO(5) model attempts to unify the superconducting and antiferromagnetic aspects of high-temperature superconductors. The order parameters for these two phenomena form a 5-dimensional real field whose dynamics has an approximate SO(5) symmetry, $`(\varphi _1,\varphi _2,\eta _1,\eta _2,\eta _3)`$, where $`\varphi =\varphi _1+i\varphi _2`$ and $`\eta =(\eta _1,\eta _2,\eta _3)`$ are, respectively, the order parameters of superconductivity and antiferromagnetism. This so-called superspin is described by a free energy which will be given below; explicit SO(5) breaking enables spontaneous breaking of either the SO(2) symmetry of $`\varphi `$ or the SO(3) symmetry of $`\eta `$, describing SC and AF, respectively.

Mathematically, the low-energy effective theory is described in terms of the following free energy :

$$\widehat{F}=𝑑𝐱\left(\frac{\widehat{𝐡}^2}{8\pi }+\frac{\mathrm{}^2}{2m^{}}\left|\left(i\frac{e^{}}{\mathrm{}c}\widehat{𝐀}\right)\varphi \right|^2+\frac{\mathrm{}^2}{2m^{}}(\eta )^2+V(\varphi ,\eta )\right),$$
(17)
where $`\widehat{𝐡}=\times \widehat{𝐀}`$ is again the microscopic magnetic field.

The nature of the ground state can be determined by examining the potential. We consider the most general symmetry-breaking potential including even powers of the fields up to quartic terms,

$$V(\varphi ,\eta )=\frac{a_1^2}{2}\varphi ^2\frac{a_2^2}{2}\eta ^2+\frac{b_1\varphi ^4+2b_3\varphi ^2\eta ^2+b_2\eta ^4}{4},$$

where we have written $`\varphi =|\varphi |`$ and $`\eta =|\eta |`$. This potential is invariant under an $`SO(3)\times SO(2)`$ symmetry. SO(5) symmetry would be attained by setting the two mass parameters and the three quartic coupling constants to be equal; for an approximate symmetry their values are approximately equal. To simplify the analysis we will set the quartic parameters to a common value $`b_1=b_2=b_3=b`$; explicit symmetry breaking in the potential is found only in the quadratic terms.

Since we are interested in the magnetic properties of the superconducting state, we restrict ourselves to parameters which describe a SC ground state. This is done by requiring that the global minimum of $`V`$ be on the $`\varphi `$ axis. The ground state will then have the value $`(\varphi ,\eta )=(v,0)`$, where $`va_1/\sqrt{b}`$; the global condition is fulfilled if $`\beta a_2^2/a_1^2<1`$. Note that we recover SO(5) symmetry if $`\beta =1`$ and, therefore, this is the value corresponding to critical doping. For $`\beta >1`$ the ground state is antiferromagnetic.

The goal of this section is to calculate the surface energy, whose sign indicates the type of superconductor described by the model. We will first set up the physical context of the surface energy and derive the critical fields $`H_c`$ and $`H_{c2}`$. These enable us to determine the value of $`\kappa _c`$ as a function of $`\beta `$. Subsequently, we will give numerical calculations which confirm our analytical results.

### 3.1 One Dimensional Free Energy

As outlined in Section 2.1, the surface energy is the Gibbs free energy per unit area of a domain wall separating a normal (non-superconducting) phase in the thermodynamic critical field $`H_c`$ and a superconducting phase in the absence of a magnetic field. We use the same ansatz for $`\varphi `$ and $`A`$ as mentioned in the beginning of Section 2.1, and with the AF order parameter written $`\eta =\eta (x)\widehat{x}`$, the free energy (17) becomes:

$`\widehat{F}={\displaystyle }dx\{f_n+{\displaystyle \frac{\mathrm{}^2}{2m^{}}}\left({\displaystyle \frac{d\varphi }{dx}}\right)^2+{\displaystyle \frac{1}{2m^{}}}\left({\displaystyle \frac{e^{}}{c}}\right)^2A^2\varphi ^2+{\displaystyle \frac{1}{8\pi }}\left({\displaystyle \frac{dA}{dx}}\right)^2`$
$`+{\displaystyle \frac{\mathrm{}^2}{2m^{}}}\left({\displaystyle \frac{d\eta }{dx}}\right)^2+V(\varphi ,\eta )\}.`$ (18)
It will again be useful to use the dimensionless quantities defined in (6). Writing $`\eta (x)=vn(s)`$, we find, after some algebra, the following dimensionless free energy:

$`F={\displaystyle }dx\{f_n+{\displaystyle \frac{1}{\kappa ^2}}(f^2+n^2)+a^2f^2+h^2f^2\beta n^2`$
$`+{\displaystyle \frac{1}{2}}(f^4+2f^2n^2+n^4)\}.`$ (19)
Minimizing with respect to $`f`$, $`n`$ and $`a`$, we obtain the SO(5) GL equations:

$$\frac{1}{\kappa ^2}\frac{d^2f}{ds^2}+\left(1a^2\right)ff^3n^2f=0,$$
(20)
$$\frac{d^2a}{ds^2}=af^2,$$
(21)
$$\frac{1}{\kappa ^2}\frac{d^2n}{ds^2}+\beta nn^3nf^2=0.$$
(22)
Note that if we set $`n=0`$ in the first two equations, we recover the conventional GL Eqs. (8-9), as expected. Note also that we are now left with two parameters ($`\kappa `$, $`\beta `$) rather than three ($`a_1^2`$, $`a_2^2`$ and $`b`$). The first of these two is the traditional GL parameter $`\kappa =\lambda /\xi `$, which is usually quite large for high $`T_c`$ SC (typically between 15 and 100). The second fundamental parameter is $`\beta `$, which is related to the doping of the system away from critical doping (that corresponding to the AF/SC transition). $`\beta `$ can be written as a function of the chemical potential $`\mu `$ as $`\beta =18m^{}\xi ^2\chi (\mu ^2\mu _c^2)/\mathrm{}^2`$, where $`\chi `$ is the charge susceptibility and $`\mu _c`$ is the critical value of the chemical potential .

### 3.2  Calculation of $`H_c`$ and $`H_{c2}`$

When a critical magnetic field is applied in a region, the SC order parameter will be forced to zero, destroying superconductivity. In order to minimize the Gibbs free energy, a nonzero AF order parameter will be induced. This is because, $`f`$ being held to zero by the magnetic field, the potential in this region is minimized if $`n=\sqrt{\beta }`$, as can be seen from the potential terms in (3.1). The value of the critical field is determined by the competition between the positive potential energy of the AF state (relative to the SC ground state), on the one hand, and the negative magnetic energy ($`H_c^2/8\pi `$), on the other.

Let us compute, first of all, the critical magnetic field $`H_c`$. The result differs from the conventional one because an AF order parameter will now appear in the normal region.

Far into the superconducting region ($`x\mathrm{}`$), we get the following behavior for the fields: $`h=0`$, $`f=1`$, $`n=0`$. Thus from Eqs. (2) and (3.1),

$$g_s(H,T)g(h=0,T)=f_s(T)=f_n(T)\frac{a_1^4}{4b}.$$
(23)
In the normal region ($`x\mathrm{}`$), we include a nonzero AF order parameter, and set $`h=H`$, $`n=\sqrt{\beta }`$, $`f=0`$, giving,

$$g_n(H,T)=f_n(T)\frac{a_2^4}{4b}\frac{H^2}{8\pi }.$$
(24)
By equating the two Gibbs free energies, we find the following $`\beta `$-dependent expression for the critical magnetic field:

$$H_c(\beta )=H_c^0\sqrt{1\beta ^2},$$
(25)
where $`H_c^0`$ is the conventional expression (4). Note that $`H_c(\beta )`$ will refer to the critical field in the SO(5) model, and must be distinguished from the one in the GL theory (from now on, we will write $`H_c^0`$). Eq. (25) tells us that for $`\beta =1`$, $`H_c(\beta )`$ vanishes. In fact, this value of $`\beta `$ corresponds to the critical value, where the transition from SC ($`\beta <1`$) to AF ($`\beta >1`$) occurs; it is natural that any magnetic field, however small, will break SC, since SC is not energetically favored over AF: their free energies are equal. As $`\beta `$ decreases, $`H_c(\beta )`$ reaches a maximum when $`\beta =0`$, i.e., when the sample is at maximal doping. In this case, the value of $`H_c(\beta )`$ is the same as in conventional GL theory with the same $`\kappa `$. Thus, far in the normal state, the critical field will take the value $`H_c(\beta )=(1/\sqrt{2})\sqrt{1\beta ^2}`$ since $`H_c^0=1/\sqrt{2}`$ in the conventional SO(2) case (see Section 2.2).

Also of interest is the critical field $`H_{c2}`$. For a type I SC, at $`T<T_c`$, when $`H>H_c^0`$, the SC phase is destroyed. After that, if $`H`$ is reduced, the SC is recovered at $`H_{c2}<H_c^0`$ (phenomenon of supercooling). For a type II SC, $`H_{c2}`$ corresponds to the maximum magnetic field in which a vortex lattice can appear; it is always higher than $`H_c^0`$. Then, $`H_{c2}`$ is the strongest field which permits the SC phase to occur.

Near $`H_{c2}`$, the order parameter is small (the system is near the transition), and we can use the linearized SO(5) GL equations. On the other hand, we can show that, at first order, the microscopic field $`𝐡`$ can be taken equal to the applied field $`𝐇`$. Beginning with the free energy Eq. (17), we choose the Landau gauge $`𝐇=H\widehat{z}`$ and $`A\widehat{y}=Hx`$. Linearizing around $`\varphi =0`$ and $`\eta =a_2/\sqrt{b}`$ (the values of fields at the transition), the equation for the SC order parameter takes the form :

$$\left[^2+\frac{4\pi i}{\mathrm{\Phi }_0}\stackrel{~}{H}x\frac{}{y}\left(\frac{2\pi \stackrel{~}{H}x}{\mathrm{\Phi }_0}\right)^2\right]\psi =\kappa ^2(1\beta )\psi ,$$

where we have set $`\stackrel{~}{H}=\lambda ^2H`$, $`\varphi =v\psi `$ and $`\mathrm{\Phi }_0`$ is the flux quantum.

This equation can be put in an harmonic oscillator equation form for $`x`$, by seting $`\psi =f(x)\mathrm{exp}\left(i(k_yy+k_zz)\right)`$ and the quantization of the energy leads to a set of allowed values of the magnetic field. The highest possible value of $`H_n`$, then, corresponds to the upper critical field and we find

$$H_{c2}(\beta )=\frac{\mathrm{\Phi }_0}{2\pi \xi ^2}(1\beta )=\sqrt{2}\kappa H_c^0(1\beta ).$$
(26)
Note that it is the conventional critical field that appears in the right hand of Eq. (26).

In the GL theory of conventional SC, we have seen that the boundary between type I and type II SC occurs when $`\kappa =\kappa _c=1/\sqrt{2}`$; at this value, $`H_{c2}=H_c^0`$. We can already anticipate the value of $`\kappa _c`$ in the SO(5) model as the value for which these two critical fields are equal.

From (25) and (26) we find:

$$\kappa _c(\beta )=(1/\sqrt{2})\sqrt{(1+\beta )/(1\beta )}.$$
(27)
Thus we are able to obtain an analytical expression for $`\kappa _c`$ as a function of the doping. Two facts emerge from Eq. (27): firstly, we recover (as we must) the conventional result when $`\beta =0`$, and secondly, for any $`\beta `$, a value of $`\kappa _c`$ exists which separates type I and type II behavior (Figure 2).

As $`\beta `$ approaches 1 (the SO(5) symmetric case), the instability of a vortex lattice (in the type I region) will be amplified, and therefore, to reach the type II region, we will need a higher $`\kappa `$. In practice, it is easier to vary the doping $`\beta `$ (or $`\mu `$) than the GL parameter $`\kappa `$ in a sample.

We can reverse this point of view: if we consider a SC at a fixed value of $`\kappa `$, we will find that, for $`\beta >\beta _c`$ (defined implicitly in (27)), the sample will exhibit type I behavior. Thus, for high values of $`\kappa `$ (i.e., $`\kappa 1`$), we also need high $`\beta `$ (i.e., $`\beta `$ approaching 1) in order to see type I behavior. We will see in the next section that the numerical calculation of the surface energy agrees perfectly with these results.

### 3.3 Surface Energy

In this section, we will derive an expression for the surface energy in the SO(5) model, and rederive from it numerically the function $`\kappa _c(\beta )`$. We will begin by analyzing qualitatively the surface energy; we will see that for $`\kappa `$ sufficiently large (small), the surface energy is negative (positive). However, even this qualitative argument indicates that the critical value of $`\kappa `$ depends on $`\beta `$, as anticipated in the previous section. We will then present numerical results since, as in the conventional model, the surface energy cannot be evaluated analytically; these numerical results agree with the conclusion of the previous section based on the critical fields.

To analyze qualitatively the surface energy, consider Fig. 3. In the normal region as $`x\mathrm{}`$, the positive contribution to the energy coming from the scalar fields $`f`$ and $`n`$ exactly cancels the negative condensation energy of the magnetic field. In the SC region, as $`x\mathrm{}`$, the scalar field and magnetic energy densities individually are zero. Inside the surface, variations in the fields contribute to the surface energy. The details of these variations will determine the sign of the surface energy.

We can characterize the fields’ variations by their characteristic lengths. From Eqs. (20-22), these are seen to be

$$l_h=1,l_f=\frac{1}{\sqrt{2}\kappa },l_n=\frac{1}{\kappa \sqrt{1\beta }}.$$
(28)
At issue is whether in going from one asymptotic regime to the other, the scalar fields or the magnetic field are capable of changing more quickly. If the scalar fields vary faster, the magnetic field will penetrate into the superconducting region, and the negative magnetic energy density in this transition zone will not be accompanied by a positive scalar energy density, resulting in a negative surface energy. If the magnetic field varies faster, the situation is reversed and the surface energy will be positive. Thus we need to compare the larger of $`l_n`$ and $`l_f`$ with $`l_h`$. We see from Eq. (28) that $`l_n`$ is always greater than $`l_f`$, so $`l_n`$ determines how quickly the scalar fields can vary. Two cases are important to us when $`\kappa 1/\sqrt{2}`$. First, when $`\beta `$ is small (Fig. 3a), we find from (28) that $`l_h>l_n`$, and as described above, the surface energy will be negative. This corresponds to the behavior of a type II SC.

On the other hand, when $`\beta `$ is large (Fig. 3b), we have the opposite situation: $`l_h<l_n`$, and the surface energy will be positive. This is the case of a type I SC.

By continuity, for any $`\kappa `$, there must exist a $`\beta `$ where the surface energy is zero. This provides an alternative way of finding $`\kappa (\beta )`$, and must agree with Eq. (27).

We define the SO(5) surface energy in the same way as in Section 2.1. Subtracting Eq. (24) from Eq. (2) we are led to

$$\widehat{G}_s=_{\mathrm{}}^{\mathrm{}}𝑑x\left\{f_s\frac{\widehat{h}H_c}{4\pi }\left(f_n\frac{a_2^4}{4b}\frac{H_c^2}{8\pi }\right)\right\}.$$
(29)
Using (3.1) we find, after some algebra,

$$\widehat{G}_s=\frac{(H_c^0)^2}{4\pi }\lambda I,$$

where

$$I=_{\mathrm{}}^{\mathrm{}}𝑑s\left\{\frac{1}{2}(1f^4n^4)f^2n^2+h^2\sqrt{2}h\sqrt{1\beta ^2}\right\},$$
(30)
with $`\lambda `$ and $`H_c^0`$ defined in Section 2. Therefore, the $`\beta `$ dependence appears only in $`I`$, and the expression reduces to the conventional one, Eq. (5), when we set $`n=\beta =0`$.

Of interest are the conditions under which this integral is zero. It will be recalled that in the SO(2) case, when $`\kappa =1/\sqrt{2}`$, not only is the integral zero but the integrand is zero – a stronger statement, and one which is explained by Bogomol’nyi . One might hope that the same feature would be found in the SO(5) model, along with an analogous analytical explanation. Unfortunately, such is not the case, and the surface energy must be studied numerically to find $`\kappa _c(\beta )`$.

We have solved Eqs. (20-22) numerically, using a relaxation algorithm. The boundary conditions appropriate for the problem are:

$$\begin{array}{cccc}x\mathrm{}:& f0,& h\sqrt{1\beta ^2}/\sqrt{2},& n\sqrt{\beta };\\ x\mathrm{}:& f1,& h0,& n0.\end{array}$$
(31)
Having found the solution, it is easy to compute the dimensionless surface energy $`I`$ as a function of the doping $`\beta `$ by numerical integration of (30). We present, in Fig. 4, the results for $`\kappa =0.707`$ and $`\kappa =7.07`$ respectively. The case $`\kappa <1/\sqrt{2}`$ is of no interest since the surface energy is positive for all $`\beta `$.

Fig. 4a corresponds to $`\kappa _c`$ for conventional SC. First of all, we see that the surface energy is zero for $`\beta =0`$, as expected since, at this particular value, we recover the SO(2) case. For all other values of $`\beta `$, we have a positive surface energy. Thus, the sample remains in the type I region for any value of the doping. This is in agreement with Eq. (27), as seen in Fig. 2. Note also that the energy reaches a maximum and decreases until it reaches zero again at $`\beta =1`$. At this point, $`h=0`$ throughout the sample (see Eq. (31)), and the ground state is SO(5) symmetric. Hence, no broken symmetry appears and for topological reasons (since the surface energy can be treated like a soliton), the surface energy must be zero.

A more interesting case is presented in Fig. 4b, corresponding to $`\kappa =7.07`$ (in the conventional case, a type II superconductor). The curve starts at $`\beta =0`$ at the conventional value and passes through zero at $`\beta 0.98`$, in agreement with the analytical result (27). Below this critical value, the sample has type II behavior since the surface energy is negative. Above, the surface energy is positive, corresponding to type I behavior, i.e., vortices in this region are not energetically stable. Observe that the surface energy is again zero for $`\beta =1`$.

We can determine in this way $`\kappa _c`$ as a function of $`\beta `$; the results are reported in Table 1, and are in excellent agreement with Eq. (27).

| $`\beta `$ | $`\kappa _c(\beta )`$ (numerical) | From (27) |
| --- | --- | --- |
| 0.3945 | 1 | 1.07 |
| 0.5288 | 1.2 | 1.27 |
| 0.6620 | 1.5 | 1.57 |
| 0.7887 | 2 | 2.06 |
| 0.8976 | 3 | 3.04 |
| 0.9410 | 4 | 4.06 |
| 0.98 | 7.07 | 7.04 |
| 0.9998 | 70.7 | 70.71 |

Table 1: Comparison of $`\kappa _c`$ obtained numerically via vortex energetics and surface energy, and analytically via Eq. (27).

### 3.4 Vortex Energetics in the SO(5) Model

As discussed in Section 2.2, the magnetic behavior of a superconductor can also be determined by studying vortex energetics . Thus, one can calculate the free energy per quantum of magnetic flux of an $`m`$-vortex $`(m)=F(m)/m`$; if this function is of positive or negative slope, the superconductor is of type II or I, respectively; zero slope indicates the critical case, which can be expressed in terms of $`\kappa _c(\beta )`$. In the SO(5) model, the numerical calculation of $`(m)`$ for $`\kappa =7.07`$ and $`70.7`$ are given, respectively, in Fig. 2c and Fig. 2d in . From the first of these, we can see that for $`\beta <0.98`$ $`(m)`$ is of positive slope, while for $`\beta >0.98`$ it is of negative slope; to a very good approximation, the slope is zero at $`\beta =0.98`$. Thus, $`\kappa _c(\beta =0.98)=7.07`$. Similarly, Fig. 2d of the same reference indicates $`\kappa _c(\beta =0.9998)=70.7`$. A comparison of this method with the analytical method discussed in Section 3.2 is shown in Table 1; the agreement is excellent.

## 4 Summary and Conclusions

In this paper, magnetic properties of superconductors were discussed, with an eye towards analyzing the SO(5) model of high-temperature superconductivity. The Ginzburg-Landau theory of conventional superconductivity was discussed in Section 2, mainly to review some fundamental and well-known results and also to establish notation. The energy density of a surface separating superconducting and normal regions at critical applied field, and its importance for determining the behavior of the superconductor, were discussed. In addition, we reviewed the approach of Bogomol’nyi wherein the vortex energy is written in a form that shows clearly that at $`\kappa _c`$, the energy is proportional to the winding number of the vortex. Finally, we discussed an alternative approach to determining the critical field which is suggested by Bogomol’nyi’s work, namely, studying the energy of an $`m`$-vortex as a function of its winding number. Both analytically and numerically, this approach can be shown to be equivalent to the surface energy approach.

In the second part, we introduced briefly the SO(5) model, and derived, in the context of that model, expressions for the critical magnetic fields, $`H_c(\beta )`$ and $`H_{c2}(\beta )`$. These enable us to calculate an analytic expression for $`\kappa _c(\beta )`$. Next, we studied the surface energy of a boundary between superconducting and non-superconducting regions in an applied critical field. Finally, we discussed briefly an alternative approach based on vortex energetics.

It is worth comparing the two numerical approaches to calculating $`\kappa _c`$, namely, by calculating the surface energy or by vortex energetics. In the former approach, one must solve the field equations (for instance, (20-22) in the SO(5) model), evaluate the surface integral (30), and vary $`\beta `$ for fixed $`\kappa `$ (or vice-versa) until the surface integral vanishes. In the latter approach, one must solve the vortex field equations (the generalization of (11) and (12) to SO(5)), evaluate the energy as a function of winding number, and determine for which $`\kappa ,\beta `$ the function $`(m)`$ is independent of $`m`$ for large $`m`$. The latter method is somewhat more demanding numerically, essentially since several values of $`m`$ must be considered for each $`\kappa ,\beta `$, though it is perhaps more intuitive, since one can think of type I behavior as due to an energetic preference for the coalescence of many single vortices into one large $`m`$-vortex.

The main new result in this paper is Eq. (27), which expresses the value of $`\kappa `$ as a function of $`\beta `$ corresponding to the boundary between type I and type II superconductors. This result is rather important: it states that as $`\beta 1`$ (i.e., as the doping is reduced to its critical value), then even with a large value of $`\kappa `$, the material would be a type I superconductor, contrary to what one would naively expect. Some implications of this have been discussed in .


format: text

text:
# Exact finite-size corrections for the square lattice Ising model with Brascamp-Kunz boundary conditions

## I Introduction

In the study of phase transitions and critical phenomena, it is extremely important to understand finite-size corrections to thermodynamical quantities. In experiments and in numerical studies of model systems, it is essential to take into account finite size effects in order to extract correct infinite-volume predictions from the data. Finite-size scaling  concerns the critical behavior of systems in which one or more directions are finite, even though microscopically large and is valuable in the analysis of experimental and numerical data in many situations, for example, for films of finite thickness. As soon as one has a finite system one must consider the question of boundary conditions on the outer surfaces or “walls” of the system. The systems under various boundary conditions have the same per-site free energy, internal energy, specific heat, etc, in the bulk limit, whereas the finite size corrections are different. To understand the effects of boundary conditions on finite-size scaling and finite-size corrections, it is valuable to study model systems, such as percolation model  and the Ising model . Therefore, in recent decades there have many investigations on finite-size scaling, finite-size corrections, and boundary effects for critical model systems . Of particular importance in such studies are exact results where the analysis can be carried out without numerical errors.

The Ising model has exact solutions on finite lattices with many kinds of boundary conditions, including cylindrical , toroidal , and Mobius strip and Klein bottle . This class also includes the special boundary conditions introduced by Brascamp and Kunz . The calculation of the exact partition function of the two-dimensional Ising model in the zero field wrapped on the cylinder was performed by Onsager in 1944 . Exploiting the exactly known partition function of the two dimensional Ising model on finite $`M\times N`$ square lattice with toroidal boundary conditions , Ferdinand and Fisher  computed the finite-size corrections to the free energy, the internal energy, and the specific heat up to order $`N^1`$. Recently, there has been much effort in understanding the behavior of finite-size corrections of the free energy, internal energy and specific heat. Izmailian and Hu  and Salas  extended the results of  for the free energy and the internal energy up to order $`N^5`$ and for the specific heat up to order $`N^3`$. Lu and Wu  obtained expressions for the partition function of the Ising model on an quadratic lattice embedded on a Mobius strip and a Klein bottle. They find finite-size corrections for free energy to order $`N^1`$. Brascamp and Kunz  calculated the partition function of the Ising model on the $`M\times 2N`$ square lattice for special boundary conditions shown in Fig. 1. Recently Janke and Kenna  has calculated the finite-size corrections of the specific heat for this boundary condition up to $`M^3`$ order. Very recently, Ivashkevich, Izmailian and Hu  provided a systematic method to compute finite-size corrections to the partition function and their derivatives of the Ising model on torus. Their approach is based on an intimate relation between the terms of the asymptotic expansion and the so-called Kronecker’s double series  which are directly related to elliptic theta functions. Expressing the final result in terms of theta functions avoids messy sums (as in some earlier works) and greatly simplifies the task of verifying the behavior of the different terms in the asymptotic expansion under duality transformation $`MN`$. Using this approach, Salas  computed the finite-size corrections to the free energy, internal energy and specific heat of the critical Ising model on a triangular and honeycomb lattices wrapped on a torus.

Using the exact partition of Ref.  and the method of Ref. , in the present paper we derive exact finite-size corrections for the free energy $`F`$ and and the specific heat $`C`$ of the critical ferromagnetic Ising model on the $`\times 2𝒩`$ square lattice with Brascamp-Kunz (BK) boundary conditions  and compared such results with those under toroidal boundary conditions . We find that when the ratio $`\xi /2=(+1)/2𝒩`$ is smaller than 1 the behaviors of finite-size corrections for $`C`$ are quite different for BK and toroidal boundary conditions; when $`\mathrm{ln}(\xi /2)`$ is larger than 3, finite-size corrections for $`C`$ in two boundary conditions approach the same values. In the limit $`𝒩\mathrm{}`$ we obtain the expansion of the free energy for infinitely long strip with BK boundary conditions. Our results are consistent with the conformal field theory prediction for the mixed boundary conditions by Cardy  although the definitions of boundary conditions in two cases are different in one side of the long strip.

This paper is organized as follows. In Sec. II we show how to lead the partition function of Ising model under Brascamp-Kunz (BK) boundary conditions to the form of partition function with twisted boundary conditions. In Sec. III asymptotic expansions of the free energy is presented. In Sec. IV expansions of the specific heat is presented. Our results are summarized and discussed in Sec. V.

## II Ising model under Brascamp-Kunz <br>boundary conditions

For the Ising model on a lattice $`G`$ of $`N`$ sites, the $`i`$-th site of the lattice for $`1iN`$ is assigned a classical spin variable $`s_i`$, which has values $`\pm 1`$. The spins interact according to the Hamiltonian

$$\beta H=J\underset{<ij>}{}s_is_j$$
(1)
where $`J`$ is exchange energy, the sum runs over the nearest neighbor pairs of spins, and $`\beta =1/k_BT`$ is the inverse temperature. The partition function of the Ising model is given by the sum over all spin configurations on the lattice

$$Z_{Ising}(J)=\underset{s}{}e^{\beta H(s)}.$$
(2)
As is mentioned in introduction there are a few boundary conditions for which the Ising model has been solved exactly. Among them is the special boundary conditions studied by Brascamp and Kunz (BK) . They considered a lattice with $`2𝒩`$ sites in the $`x`$ direction and $``$ sites in the $`y`$ direction. The boundary conditions are periodic in the $`x`$ direction; in the $`y`$ direction, the spins are up (+1) along the upper border of the resulting cylinder and have the alternative values along the lower border of the resulting cylinder as is shown in Fig. 1. Recently Janke and Kenna  has analyzed the Ising model in two dimensions with these boundary conditions. They have derived exact expressions for the finite-size scaling of the specific heat up to the $`^3`$ order. In this paper we obtain all exact finite-size corrections for the free energy and the specific heat. Moreover in our case the terms of asymptotic expansion are analytical functions. They are related to Kronecker’s double series  which in turn can be expressed by elliptic $`\theta `$-functions .

For the BK boundary conditions, the Ising partition function given in Ref.  can be rewritten as

$$Z_{,2𝒩}=(\sqrt{2}e^\mu )^{2𝒩}\underset{i=1}{\overset{𝒩}{}}\underset{j=1}{\overset{}{}}F(i,j)$$
(3)
where $`\mu =1/2\mathrm{ln}\mathrm{sinh}2J`$ and

$$F(i,j)=4\left[2\mathrm{sinh}^2\mu +\mathrm{sin}^2\left(\frac{\pi (i1/2)}{2𝒩}\right)+\mathrm{sin}^2\left(\frac{\pi j}{2(+1)}\right)\right].$$
(4)
Now we try to express the partition function $`Z_{,2𝒩}`$ given by Eq. (3) to the form of partition function with twisted boundary conditions $`Z_{\alpha ,\beta }(\mu )`$

$$Z_{\alpha ,\beta }^2(\mu )=\underset{i=0}{\overset{N1}{}}\underset{j=0}{\overset{M1}{}}4\left[\mathrm{sin}^2\left(\frac{\pi (i+\alpha )}{N}\right)+\mathrm{sin}^2\left(\frac{\pi (j+\beta )}{M}\right)+2\mathrm{sinh}^2\mu \right]$$
(5)
for which a general theory about its asymptotic expansion has been given in Ref. . For this purpose we can express the double products
$`_{i=0}^{2𝒩1}_{j=0}^{2+1}F(i+1,j)`$ through $`_{i=1}^𝒩_{j=1}^{}F(i,j)`$ as

$$\underset{i=0}{\overset{2𝒩1}{}}\underset{j=0}{\overset{2+1}{}}F(i+1,j)=\left(\underset{i=1}{\overset{𝒩}{}}\underset{j=1}{\overset{}{}}F(i,j)\right)^4\underset{i=0}{\overset{2𝒩1}{}}F(i+1,0)F(i+1,+1).$$
(6)
Here we use the properties of function $`F(i,j)`$

$$F(2𝒩+1k,j)=F(k,j)\text{and}F(i,2+2k)=F(i,k).$$
(7)
This transformation leads the rectangular lattice $`\times 2𝒩`$ under consideration to the lattice $`2(+1)\times 2𝒩`$. In what follows we will use for convenience the definition of the aspect ratio as $`\xi =\frac{+1}{𝒩}`$ instead of conventional one ($`\xi =\frac{}{2𝒩}`$).

The left-hand side of Eq. (6) is nothing but the partition function with twisted boundary conditions $`Z_{1/2,0}^2(\mu )`$ given by Eq. (5) with $`N=2𝒩`$ and $`M=2(+1)`$. With the help of the identity

$$\underset{i=0}{\overset{2𝒩1}{}}4\left[\mathrm{sinh}^2\omega +\mathrm{sin}^2\left(\frac{\pi (i+1/2)}{2𝒩}\right)\right]=4\mathrm{cosh}^2(2𝒩\omega )$$

the second product in the right-hand side of Eq. (6) can be transformed into the form

$$\underset{i=0}{\overset{2𝒩1}{}}F(i+1,0)F(i+1,+1)=\left[4\mathrm{cosh}\left(2𝒩\omega _\mu (0)\right)\mathrm{cosh}\left(2𝒩\omega _\mu (\pi /2)\right)\right]^2$$
(8)
where

$$\omega _\mu (k)=\mathrm{arcsinh}\sqrt{\mathrm{sin}^2k+2\mathrm{sh}^2\mu }$$
(9)
is a lattice dispersion relation.

By the aid of Eqs. (4) - (6) and (8) the partition function $`Z_{,2𝒩}`$ can be expressed as

$$Z_{,2𝒩}^2=\frac{\left(\sqrt{2}e^\mu \right)^{4𝒩}}{4\mathrm{cosh}\left[2𝒩\omega _\mu (0)\right]\mathrm{cosh}\left[2𝒩\omega _\mu (\pi /2)\right]}Z_{1/2,0}(\mu ).$$
(10)
For our further purposes we transform the partition function $`Z_{1/2,0}`$ into the simpler form

$$Z_{1/2,0}(\mu )=\underset{n=0}{\overset{N1}{}}2\mathrm{sinh}\left[M\omega _\mu \left(\frac{\pi (n+1/2)}{N}\right)\right]$$
(11)
where $`N=2𝒩`$ and $`M=2(+1)`$.

## III Asymptotic Expansion of the Free Energy

In the previous section it was shown that the partition function of the $`\times 2𝒩`$ Ising model with BK boundary conditions can be expressed in terms of the partition function with twisted boundary conditions $`Z_{1/2,0}`$, which has been well studied in Ref. . Farther we will use it and for simplicity we will remind some necessary parts from there. For reader’s convenience, all the technical details of our calculations and definitions of the special functions are summarized in the appendices at the end of this paper. Considering the logarithm of the partition function with twisted boundary conditions, Eq. (11), we note, that it can be transformed as

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

$$M\underset{n=0}{\overset{N1}{}}\omega _0\left(\frac{\pi (n+1/2)}{N}\right)=\frac{S}{\pi }_0^\pi \omega _0(x)dx\pi \xi \mathrm{B}_2^{1/2}2\pi \xi \underset{p=1}{\overset{\mathrm{}}{}}\left(\frac{\pi ^2\xi }{S}\right)^p\frac{\lambda _{2p}}{(2p)!}\frac{\mathrm{B}_{2p+2}^{1/2}}{2p+2}.$$
(13)
Here $`S=NM=4𝒩(+1)`$, $`\mathrm{B}_p^{1/2}`$ are the Bernoulli polynomials $`\mathrm{B}_p^\alpha `$ at $`\alpha =1/2`$ which are related to the Bernoulli numbers $`B_pB_p^0`$ as $`B_p^{1/2}=(2^{1p}1)B_p`$ and $`\xi =\frac{M}{N}=\frac{+1}{𝒩}`$. We have also used the symmetry property, $`\omega _0(k)=\omega _0(\pi k)`$, of the lattice dispersion relation given by Eq. (9) and its Taylor expansion

$$\omega _0(k)=k\left(\lambda +\underset{p=1}{\overset{\mathrm{}}{}}\frac{\lambda _{2p}}{(2p)!}k^{2p}\right),$$
(14)
where $`\lambda =1`$, $`\lambda _2=2/3`$, $`\lambda _4=4`$, etc.

We may transform the second term in Eq. (12) as

$`{\displaystyle \underset{n=0}{\overset{N1}{}}}\mathrm{ln}\left(\mathrm{\hspace{0.17em}1}e^{2M\omega _0\left(\frac{\pi (n+1/2)}{N}\right)}\right)=\mathrm{ln}{\displaystyle \frac{\theta _4}{\eta }}+\pi \xi \mathrm{B}_2^{1/2}`$ (15)
$`2\pi \xi {\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\pi ^2\xi }{S}}\right)^p{\displaystyle \frac{\mathrm{\Lambda }_{2p}}{(2p)!}}{\displaystyle \frac{\mathrm{𝚁𝚎}\mathrm{K}_{2p+2}^{1/2,0}(i\lambda \xi )\mathrm{B}_{2p+2}^{1/2}}{2p+2}},`$ (16)
where $`\eta =(\theta _2\theta _3\theta _4/2)^{1/3}`$ is the Dedekind-$`\eta `$ function; $`\theta _2`$, $`\theta _3`$, $`\theta _4`$ are elliptic $`\theta `$-functions and $`\mathrm{K}_{2p+2}^{1/2,0}(i\lambda \xi )`$ are Kronecker’s double series  (see also Appendix A). Taking into account the relation between moments and cumulants (Appendix B), 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

$`\mathrm{\Lambda }_2`$ $`=`$ $`\lambda _2,`$ (17)
$`\mathrm{\Lambda }_4`$ $`=`$ $`\lambda _4+3\lambda _2^2{\displaystyle \frac{}{\lambda }},`$ (18)
$`\mathrm{\Lambda }_6`$ $`=`$ $`\lambda _6+15\lambda _4\lambda _2{\displaystyle \frac{}{\lambda }}+15\lambda _2^3{\displaystyle \frac{^2}{\lambda ^2}},`$ (19)
$`\mathrm{}`$ (20)
$`\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}}.`$ (21)
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`$.

Substituting Eqs. (13) and (16) into Eq. (12) 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_{1/2,0}(0)`$ $`=`$ $`{\displaystyle \frac{S}{\pi }}{\displaystyle _0^\pi }\omega _0(x)dx+\mathrm{ln}{\displaystyle \frac{\theta _4}{\eta }}`$ (23)
$`2\pi \xi {\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\pi ^2\xi }{S}}\right)^p{\displaystyle \frac{\mathrm{\Lambda }_{2p}}{(2p)!}}{\displaystyle \frac{\mathrm{𝚁𝚎}\mathrm{K}_{2p+2}^{1/2,0}(i\lambda \xi )}{2p+2}},`$
where $`_0^\pi \omega _0(x)dx=2G`$ and $`G=0.915966`$ is the Catalan’s constant.

After reaching this point, one can easily write down the exact asymptotic expansion of the free energy, $`F=\mathrm{ln}Z_{,2𝒩}`$, at the critical point. Plugging Eq. (23) back in Eq. (10) we have finally obtain

$`F`$ $`=`$ $`2𝒩\left({\displaystyle \frac{1}{2}}\mathrm{ln}2+{\displaystyle \frac{2G}{\pi }}\right)+2𝒩\left[{\displaystyle \frac{1}{2}}\mathrm{ln}(1+\sqrt{2}){\displaystyle \frac{2G}{\pi }}\right]{\displaystyle \frac{1}{2}}\mathrm{ln}{\displaystyle \frac{\theta _4}{2\eta }}`$ (24)
$`+`$ $`\pi \xi {\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\pi }{2𝒩}}\right)^{2p}{\displaystyle \frac{\mathrm{\Lambda }_{2p}}{(2p)!}}{\displaystyle \frac{\mathrm{𝚁𝚎}\mathrm{K}_{2p+2}^{1/2,0}(i\lambda \xi )}{2p+2}}.`$ (25)
Note that the Kroneker’s functions $`K_p^{\frac{1}{2},0}(i\lambda \xi )`$ can be expressed in terms of the elliptic $`\theta `$-function only. Thus, Eq. (25) can be rewritten in the following form

$$F=2𝒩f_{bulk}+2𝒩f_1+f_0+\underset{p=1}{\overset{\mathrm{}}{}}\frac{f_{2p}}{\left(2𝒩\right)^{2p}},$$
(26)
where

$`f_{bulk}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}2{\displaystyle \frac{2G}{\pi }}=0.929695\mathrm{},`$ (27)
$`f_1`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}(1+\sqrt{2}){\displaystyle \frac{2G}{\pi }}=0.142435\mathrm{},`$ (28)
$`f_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}{\displaystyle \frac{\theta _4}{2\eta }},`$ (29)
$`f_2`$ $`=`$ $`{\displaystyle \frac{\pi ^3\xi }{360}}\left({\displaystyle \frac{7}{8}}\theta _4^8\theta _2^4\theta _3^4\right),`$ (30)
$`f_4`$ $`=`$ $`{\displaystyle \frac{\pi ^5\xi }{48384}}[\pi \xi \theta _3^4\theta _4^4(\theta _3^8+{\displaystyle \frac{5}{4}}\theta _3^4\theta _4^4{\displaystyle \frac{5}{16}}\theta _4^8)`$ (32)
$`+(\theta _2^4+\theta _3^4)({\displaystyle \frac{31}{16}}\theta _4^8+\theta _2^4\theta _3^4)(1+4\xi {\displaystyle \frac{\theta _2^{}}{\theta _2}})],`$
$`f_6`$ $`=`$ $`{\displaystyle \frac{\pi ^7\xi }{87091200}}[70\pi ^2\xi ^2\theta _3^4\theta _4^4(\theta _2^{16}+{\displaystyle \frac{\theta _2^{12}\theta _4^4}{2}}{\displaystyle \frac{81\theta _2^8\theta _4^8}{8}}`$ (34)
$`{\displaystyle \frac{295\theta _2^4\theta _4^{12}}{16}}{\displaystyle \frac{635\theta _4^{16}}{64}})`$
$`+`$ $`630\pi \xi \theta _3^4\theta _4^4\left(\theta _2^{12}{\displaystyle \frac{3\theta _2^8\theta _4^4}{4}}{\displaystyle \frac{21\theta _2^4\theta _4^8}{8}}{\displaystyle \frac{127\theta _4^{12}}{64}}\right)\left(1+4\xi {\displaystyle \frac{\theta _2^{}}{\theta _2}}\right)`$ (35)
$`+`$ $`\left(\theta _2^{16}+2\theta _2^{12}\theta _4^4{\displaystyle \frac{3}{4}}\theta _2^8\theta _4^8{\displaystyle \frac{7}{4}}\theta _2^4\theta _4^{12}{\displaystyle \frac{127}{128}}\theta _4^{16}\right)`$ (37)
$`(711+5040\xi {\displaystyle \frac{\theta _2^{}}{\theta _2}}+8400\xi ^2{\displaystyle \frac{\theta _2^2}{\theta _2^2}}+560\xi ^2{\displaystyle \frac{\theta _2^{\prime \prime }}{\theta _2}})],`$
$`f_8`$ $`=`$ $`\pi ^{12}\xi ^4{\displaystyle \frac{\theta _4^4\theta _3^4}{33634123776}}(1280\theta _2^{24}+20224\theta _2^{20}\theta _4^4+83664\theta _2^{16}\theta _4^8`$ (38)
$`+`$ $`210496\theta _2^{12}\theta _4^{12}+361115\theta _2^8\theta _4^{16}+323910\theta _2^4\theta _4^{20}+107310\theta _4^{24})`$ (39)
$``$ $`\pi ^{11}\xi ^3{\displaystyle \frac{\theta _4^4\theta _3^4\left(1+4\xi \frac{\theta _2^{}}{\theta _2}\right)}{1868562432}}(1280\theta _2^{20}+8832\theta _2^{16}\theta _4^4+19056\theta _2^{12}\theta _4^8`$ (40)
$`+`$ $`33568\theta _2^8\theta _4^{12}+38655\theta _2^4\theta _4^{16}+15330\theta _4^{20})`$ (41)
$``$ $`{\displaystyle \frac{\pi ^{10}\xi ^2\theta _4^4\theta _3^4}{1177194332160}}\left(3789+27720\xi {\displaystyle \frac{\theta _2^{}}{\theta _2}}+48720\xi ^2{\displaystyle \frac{\theta _2^2}{\theta _2^2}}+2240\xi ^2{\displaystyle \frac{\theta _2^{\prime \prime }}{\theta _2}}\right)`$ (43)
$`\left(1280\theta _2^{16}+3136\theta _2^{12}\theta _4^4+3216\theta _2^8\theta _4^8+5176\theta _2^4\theta _412+2555\theta _4^{16}\right)`$
$``$ $`{\displaystyle \frac{\pi ^9\xi \left(\theta _2^4+\theta _3^4\right)}{235438866432}}(1479+15156\xi {\displaystyle \frac{\theta _2^{}}{\theta _2}}+47880\xi ^3{\displaystyle \frac{\theta _2^2}{\theta _2^2}}+47880\xi ^3{\displaystyle \frac{\theta _2^3}{\theta _2^3}}`$ (44)
$`+`$ $`2520\xi ^2{\displaystyle \frac{\theta _2^{\prime \prime }}{\theta _2}}+7980\xi ^3{\displaystyle \frac{\theta _2^{}}{\theta _2}}{\displaystyle \frac{\theta _2^{\prime \prime }}{\theta _2}}+140{\displaystyle \frac{\theta _2^{\prime \prime \prime }}{\theta _2}}).`$ (45)
$`\mathrm{}`$ (46)
The free energy per unit length of an infinitely long strip of width $`L`$ at critically has the finite-size scaling form

$$F=fL+f^{}+\frac{\mathrm{\Delta }}{L}+\mathrm{},$$
(47)
where $`f`$ is the bulk free energy per unit area, $`\frac{1}{2}f^{}`$ is the surface energy, $`L^1`$ is a scaling field, and $`\mathrm{\Delta }`$ is an universal constant which depends only on the type of boundary conditions ,

$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{\pi }{12}}\text{periodic boundary conditions},`$ (48)
$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{\pi }{6}}\text{antiperiodic boundary conditions},`$ (49)
$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{\pi }{48}}\text{free boundary conditions},`$ (50)
$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{\pi }{48}}\text{fixed }++\text{ boundary conditions},`$ (51)
$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{23\pi }{48}}\text{fixed }+\text{ boundary conditions},`$ (52)
$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{\pi }{24}}\text{mixed boundary conditions}.`$ (53)
For fixed $`++`$ (or $`+`$) boundary conditions the spins are fixed to the same (or opposite) values on two sides of the strip. The mixed boundary conditions corresponds to free boundary conditions on one side of the strip, and fixed boundary conditions on the other. Therefore, BK and the mixed boundary conditions are the same in one side of the long strip (fixed to + for all spins) and they are different in another side of the long strip (fixed to $`++\mathrm{}`$ for BK boundary conditions and free boundary conditions for the mixed boundary conditions).

Using Kronecker’s functions asymptotic form when $`\xi 0`$ and $`\xi \mathrm{}`$ we can obtain from Eq. (25) the free energy per unit length of an infinitely long strip of finite width. In the limit $`\xi \mathrm{}`$ (i.e. $`\mathrm{}`$) for fixed $`2𝒩`$ from Eq. (25) one obtains the free energy expansion for infinitely long cylinder of circumference $`2𝒩`$

$`\underset{\mathrm{}}{lim}{\displaystyle \frac{F}{}}`$ $`=`$ $`2𝒩f_{bulk}{\displaystyle \frac{\pi }{24𝒩}}+2{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\pi }{2𝒩}}\right)^{2p+1}{\displaystyle \frac{\lambda _{2p}B_{2p+2}^{1/2}}{(2p)!(2p+2)}}`$ (54)
$`=`$ $`2𝒩\left({\displaystyle \frac{2G}{\pi }}{\displaystyle \frac{1}{2}}\mathrm{ln}2\right){\displaystyle \frac{\pi }{12}}\left({\displaystyle \frac{1}{2𝒩}}\right){\displaystyle \frac{7\pi ^3}{1440}}\left({\displaystyle \frac{1}{2𝒩}}\right)^3`$ (55)
$``$ $`{\displaystyle \frac{31\pi ^5}{24192}}\left({\displaystyle \frac{1}{2𝒩}}\right)^5{\displaystyle \frac{10033\pi ^7}{9676800}}\left({\displaystyle \frac{1}{2𝒩}}\right)^7\mathrm{}.`$ (56)
This result coincides with that obtained in  with the leading finite-size correction to free energy $`\frac{\pi }{12}(2𝒩)^1`$. In the limit $`\xi 0`$ (i.e. $`𝒩\mathrm{}`$) for fixed $``$ we obtain the expansion of free energy of infinitely long strip with BK boundary condition of the width $``$

$`\underset{𝒩\mathrm{}}{lim}{\displaystyle \frac{F}{2𝒩}}`$ $`=`$ $`f_{bulk}+f_1+{\displaystyle \frac{\pi }{24(+1)}}`$ (57)
$`+`$ $`{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\left[{\displaystyle \frac{\pi }{2(+1)}}\right]^{2p+1}{\displaystyle \frac{\lambda _{2p}B_{2p+2}}{(2p)!(2p+2)}}`$ (58)
$`=`$ $`\left({\displaystyle \frac{2G}{\pi }}{\displaystyle \frac{1}{2}}\mathrm{ln}2\right)+{\displaystyle \frac{1}{2}}\mathrm{ln}(1+\sqrt{2}){\displaystyle \frac{2G}{\pi }}+{\displaystyle \frac{\pi }{24}}{\displaystyle \frac{1}{}}`$ (59)
$``$ $`{\displaystyle \frac{\pi }{24}}{\displaystyle \frac{1}{^2}}+\left({\displaystyle \frac{\pi }{24}}+{\displaystyle \frac{\pi ^3}{2800}}\right){\displaystyle \frac{1}{^3}}\left({\displaystyle \frac{\pi }{24}}+{\displaystyle \frac{\pi ^3}{960}}\right){\displaystyle \frac{1}{^4}}+\mathrm{}.`$ (60)
Here the leading finite-width correction to free energy is $`\frac{\pi }{24}^1`$. ¿From Eqs. (47) - (58) one can see, that

$`1.`$ $`L=2𝒩,f=f_{bulk},f^{}=0,\mathrm{\Delta }={\displaystyle \frac{\pi }{12}},`$ (61)
$`2.`$ $`L=,f=f_{bulk},f^{}=f_1,\mathrm{\Delta }={\displaystyle \frac{\pi }{24}}.`$ (62)
Our results are consistent with the conformal field theory prediction for the mixed boundary condition (see Eq. (53)) although the mixed boundary condition and the BK boundary condition are different in one side of the long strip.

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

The internal energy per spin and the specific heat per spin can be obtained from the partition function $`Z_{,2𝒩}`$

$`U`$ $`=`$ $`{\displaystyle \frac{1}{2𝒩}}{\displaystyle \frac{d}{dJ}}\mathrm{ln}Z_{,2𝒩}={\displaystyle \frac{\sqrt{1+e^{4\mu }}}{2𝒩}}{\displaystyle \frac{d}{d\mu }}\mathrm{ln}Z_{,2𝒩},`$ (63)
$`C`$ $`=`$ $`{\displaystyle \frac{1}{2𝒩}}{\displaystyle \frac{d^2}{dJ^2}}\mathrm{ln}Z_{,2𝒩}`$ (64)
$`=`$ $`{\displaystyle \frac{e^{4\mu }}{𝒩}}\left({\displaystyle \frac{1+e^{4\mu }}{2}}{\displaystyle \frac{d^2}{d\mu ^2}}\mathrm{ln}Z_{,2𝒩}{\displaystyle \frac{d}{d\mu }}\mathrm{ln}Z_{,2𝒩}\right).`$ (65)
Let us first consider the internal energy. At the critical point $`T=T_c`$ ($`\mu =0`$) the internal energy is given by

$$U=\sqrt{2}+\sqrt{2}\frac{d}{d\mu }\mathrm{ln}Z_{1/2,0}(0).$$
(66)
One can note that $`Z_{1/2,0}(\mu )`$ is an even function with respect to its argument $`\mu `$, which imply immediately that $`\left(dZ_{1/2,0}(\mu )/d\mu \right)_{\mu =0}=0`$. Thus we find that internal energy for the finite system is equal to its bulk values without any finite-size corrections, namely $`U=\sqrt{2}`$.

At the critical point $`T=T_c`$ ($`\mu =0`$) the specific heat is given by

$$C=2\frac{4𝒩}{}\frac{\sqrt{2}}{}\mathrm{tanh}\left[2𝒩\mathrm{ln}(1+\sqrt{2})\right]+\frac{1}{2𝒩}\frac{d^2}{d\mu ^2}\mathrm{ln}Z_{1/2,0}(0).$$
(67)
The analysis of the $`Z_{1/2,0}^{\prime \prime }(0)`$ is a little more involved. Taking the second derivative of Eq.(11) with respect to mass variable $`\mu `$ and then considering the limit $`\mu 0`$, we obtain

$`{\displaystyle \frac{Z_{1/2,0}^{\prime \prime }(0)}{Z_{1/2,0}(0)}}`$ $`=`$ $`M{\displaystyle \underset{n=0}{\overset{N1}{}}}\omega _0^{\prime \prime }\left(\frac{\pi (n+1/2)}{N}\right)\mathrm{cth}\left[M\omega _0\left(\frac{\pi (n+1/2)}{N}\right)\right]`$ (68)
$`=`$ $`M{\displaystyle \underset{n=0}{\overset{N1}{}}}\omega _0^{\prime \prime }\left(\frac{\pi (n+1/2)}{N}\right)+2M{\displaystyle \underset{n=0}{\overset{N1}{}}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\omega _0^{\prime \prime }\left(\frac{\pi (n+1/2)}{N}\right)e^{2m\left[M\omega _0\left(\frac{\pi (n+1/2)}{N}\right)\right]},`$ (69)
where $`M=2(+1)`$, $`N=2𝒩`$, and $`\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}}}.`$ (70)
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{2}{x}}\left\{1+{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\kappa _{2p}}{(2p)!}}x^{2p}\right\},`$ (71)
where $`\kappa _2=2/3`$, $`\kappa _4=172/15`$, etc. The first sum in Eq. (68) we may transform as

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

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

$$\underset{n=0}{\overset{N1}{}}\frac{1}{n+1/2}=\psi (N+1/2)\psi (1/2).$$
(75)
The asymptotic expansion of the digamma function $`\psi (x)`$ is given by (see Appendix D)

$$\psi (N+1/2)=\mathrm{ln}N\underset{p=1}{\overset{\mathrm{}}{}}(1)^p\frac{\mathrm{B}_p^{1/2}}{p}\frac{1}{N^p}.$$
(76)
Using the property of the Bernoulli polynomials $`\mathrm{B}_p^{1/2}`$, namely, $`\mathrm{B}_{2p+1}^{1/2}=0`$, Eq. (75) can be rewritten as

$$\underset{n=0}{\overset{N1}{}}\frac{1}{n+1/2}=\mathrm{ln}N\underset{p=1}{\overset{\mathrm{}}{}}\frac{\mathrm{B}_{2p}^{1/2}}{2p}\frac{1}{N^{2p}}\psi (1/2).$$
(77)
Plugging Eqs. (74) and (77) back in Eq. (72), we have finally obtain

$`M{\displaystyle \underset{n=0}{\overset{N1}{}}}\omega _0^{\prime \prime }\left(\frac{\pi (n+1/2)}{N}\right)`$ $`=`$ $`{\displaystyle \frac{4S}{\pi }}\left\{\mathrm{ln}N+{\displaystyle \frac{1}{2}}\mathrm{ln}2\mathrm{ln}\pi \psi (1/2)\right\}`$ (78)
$``$ $`2\pi \xi {\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\pi ^2\xi }{S}}\right)^{p1}{\displaystyle \frac{\kappa _{2p}\mathrm{B}_{2p}^{1/2}}{p(2p)!}}.`$ (79)
Let us now consider the second sum in Eq. (68). Note that function $`\omega _0^{\prime \prime }(x)`$ can be represented as

$$\omega _0^{\prime \prime }(x)=\frac{2}{x}\mathrm{exp}\left\{\underset{p=1}{\overset{\mathrm{}}{}}\frac{\epsilon _{2p}}{(2p)!}x^{2p}\right\},$$
(80)
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), the second sum in Eq. (68) can be written as

$`2M{\displaystyle \underset{n=0}{\overset{N1}{}}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\omega _0^{\prime \prime }\left(\frac{\pi (n+1/2)}{N}\right)e^{2mM\omega _0\left(\frac{\pi (n+1/2)}{N}\right)}=`$ (81)
$`={\displaystyle \frac{4S}{\pi }}\left\{R_{1/2,0}(\xi )+\psi (1/2)\right\}+\left(\kappa _2\xi {\displaystyle \frac{}{\xi }}+\lambda _2\xi ^2{\displaystyle \frac{^2}{\xi ^2}}\right)\mathrm{ln}{\displaystyle \frac{\theta _4(\xi )}{\eta (\xi )}}`$ (82)
$`2\pi \xi {\displaystyle \underset{p=2}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{\Omega }_{2p}}{p(2p)!}}\left({\displaystyle \frac{\pi ^2\xi }{S}}\right)^{p1}\mathrm{𝚁𝚎}\mathrm{K}_{2p}^{1/2,0}(i\lambda \xi )+2\pi \xi {\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\kappa _{2p}\mathrm{B}_{2p}^{1/2}}{p(2p)!}}\left({\displaystyle \frac{\pi ^2\xi }{S}}\right)^{p1},`$ (83)
where

$`R_{1/2,0}(\xi )=2\mathrm{ln}\theta _4(\xi )+C_E+2\mathrm{ln}2`$ (84)
and $`C_E`$ is the Euler constant. 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,`$ (85)
$`\mathrm{\Omega }_4`$ $`=`$ $`\omega _4+3\omega _2^2,`$ (86)
$`\mathrm{}`$ (87)
Substituting Eqs. (79) and (83) into Eq. (68), we obtain exact asymptotic expansion of $`Z_{1/2,0}^{\prime \prime }(0)`$

$`{\displaystyle \frac{Z_{1/2,0}^{\prime \prime }(0)}{Z_{1/2,0}(0)}}`$ $`=`$ $`{\displaystyle \frac{4S}{\pi }}\left(\mathrm{ln}N+C_E+\mathrm{ln}{\displaystyle \frac{2^{5/2}}{\pi }}2\mathrm{ln}\theta _4(\xi )\right)`$ (88)
$`+`$ $`\left(\kappa _2\xi {\displaystyle \frac{}{\xi }}+\lambda _2\xi ^2{\displaystyle \frac{^2}{\xi ^2}}\right)\mathrm{ln}{\displaystyle \frac{\theta _4(\xi )}{\eta (\xi )}}`$ (89)
$``$ $`2\pi \xi {\displaystyle \underset{p=2}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{\Omega }_{2p}}{p(2p)!}}\left({\displaystyle \frac{\pi ^2\xi }{S}}\right)^{p1}\mathrm{𝚁𝚎}\mathrm{K}_{2p}^{1/2,0}(i\lambda \xi ).`$ (90)
Plugging Eq. (89) back in Eq. (67) we finally obtain exact asymptotic expansion of the specific heat

$$C=\frac{8}{\pi }\left(1+\frac{1}{\xi 𝒩1}\right)\mathrm{ln}2𝒩+\underset{p=0}{\overset{\mathrm{}}{}}\frac{C_p}{(2𝒩)^p},$$
(91)
where

$`C_0`$ $`=`$ $`{\displaystyle \frac{8}{\pi }}\left(C_E+\mathrm{ln}{\displaystyle \frac{2^{5/2}}{\pi }}{\displaystyle \frac{\pi }{4}}\right){\displaystyle \frac{4}{\xi }}{\displaystyle \frac{16}{\pi }}\mathrm{ln}\theta _4(\xi ),`$ (92)
$`C_1`$ $`=`$ $`{\displaystyle \frac{2}{\xi }}(C_0+2\sqrt{2}),`$ (93)
$`C_2`$ $`=`$ $`{\displaystyle \frac{2}{\xi }}C_1{\displaystyle \frac{\pi }{9}}\left\{\pi \xi \theta _3^4\theta _4^4+(\theta _2^4+\theta _3^4)\left(1+4\xi {\displaystyle \frac{\theta _2^{}}{\theta _2}}\right)\right\},`$ (94)
$`C_3`$ $`=`$ $`{\displaystyle \frac{2}{\xi }}C_2,`$ (95)
$`C_4`$ $`=`$ $`{\displaystyle \frac{2}{\xi }}C_3+{\displaystyle \frac{\pi ^5\xi ^4\theta _4^4\theta _3^4}{270}}\left(\theta _2^8{\displaystyle \frac{3}{2}}\theta _4^4\theta _2^4{\displaystyle \frac{21}{4}}\theta _4^8\right)`$ (96)
$`+`$ $`{\displaystyle \frac{\pi ^4\xi ^3\theta _4^4\theta _3^4}{54}}\left(\theta _2^4{\displaystyle \frac{7}{4}}\theta _4^4\right)\left(1+4\xi {\displaystyle \frac{\theta _2^{}}{\theta _2}}\right)`$ (97)
$`+`$ $`{\displaystyle \frac{4\pi ^3\xi ^2}{135}}\left(\theta _2^4\theta _3^4{\displaystyle \frac{7}{8}}\theta _4^8\right)\left({\displaystyle \frac{43}{40}}+5\xi {\displaystyle \frac{\theta _2^{}}{\theta _2}}+7\xi ^2{\displaystyle \frac{\theta _2^2}{\theta _2^2}}+\xi ^2{\displaystyle \frac{\theta _2^{\prime \prime }}{\theta _2}}\right),`$ (98)
$`C_5`$ $`=`$ $`{\displaystyle \frac{2}{\xi }}C_4,`$ (99)
$`\mathrm{}`$ (100)
The $`\frac{1}{}`$ expansion of the specific heat has a form

$$C=\frac{8}{\pi }\left(1+\frac{1}{}\right)\mathrm{ln}+\underset{p=0}{\overset{\mathrm{}}{}}\frac{c_p}{^p},$$
(101)
where

$`c_0`$ $`=`$ $`C_0{\displaystyle \frac{8}{\pi }}\mathrm{ln}{\displaystyle \frac{\xi }{2}},`$ (102)
$`c_1`$ $`=`$ $`{\displaystyle \frac{\xi }{2}}C_1+{\displaystyle \frac{8}{\pi }}(1\mathrm{ln}{\displaystyle \frac{\xi }{2}}),`$ (103)
$`c_2`$ $`=`$ $`{\displaystyle \frac{\xi ^2}{4}}C_2{\displaystyle \frac{\xi }{2}}C_1+{\displaystyle \frac{4}{\pi }},`$ (104)
$`c_3`$ $`=`$ $`{\displaystyle \frac{\xi ^2}{4}}C_2+{\displaystyle \frac{\xi }{2}}C_1{\displaystyle \frac{4}{3\pi }},`$ (105)
$`c_4`$ $`=`$ $`{\displaystyle \frac{\xi ^4}{16}}C_4{\displaystyle \frac{\xi }{2}}C_1+{\displaystyle \frac{2}{3\pi }},`$ (106)
$`c_5`$ $`=`$ $`{\displaystyle \frac{3\xi ^4}{16}}C_4+{\displaystyle \frac{\xi ^2}{2}}C_2+{\displaystyle \frac{\xi }{2}}C_1{\displaystyle \frac{2}{5\pi }},`$ (107)
$`\mathrm{}`$ (108)
Typical value of the constants $`c_0`$ \- $`c_3`$ are given in Table 1, in which the coefficients are consistent with those obtained in .

Using Kronecker’s functions asymptotic form when $`\xi 0`$ and $`\xi \mathrm{}`$ we can obtain from Eqs. (67) and (89) the specific heat per unit length of an infinitely long strip of finite width. In the limit $`\xi \mathrm{}`$ (i.e. $`\mathrm{}`$) for fixed $`2𝒩`$ the specific heat expansion for infinitely long cylinder of circumference $`2𝒩`$ can be written as

$`C`$ $`=`$ $`{\displaystyle \frac{8}{\pi }}\mathrm{ln}2𝒩+{\displaystyle \frac{8}{\pi }}\left(C_E+\mathrm{ln}{\displaystyle \frac{2^{5/2}}{\pi }}{\displaystyle \frac{\pi }{4}}\right){\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{4\pi ^{2p1}\mathrm{\Omega }_{2p}B_{2p}^{1/2}}{p(2p)!}}{\displaystyle \frac{1}{(2𝒩)^{2p}}}`$ (109)
$`=`$ $`{\displaystyle \frac{8}{\pi }}\mathrm{ln}2𝒩+{\displaystyle \frac{8}{\pi }}\left(C_E+\mathrm{ln}{\displaystyle \frac{2^{5/2}}{\pi }}{\displaystyle \frac{\pi }{4}}\right){\displaystyle \frac{\pi }{9}}\left({\displaystyle \frac{1}{2𝒩}}\right)^2`$ (110)
$``$ $`{\displaystyle \frac{301\pi ^3}{10800}}\left({\displaystyle \frac{1}{2𝒩}}\right)^4{\displaystyle \frac{29419\pi ^5}{1905120}}\left({\displaystyle \frac{1}{2𝒩}}\right)^6{\displaystyle \frac{2759329\pi ^7}{145152000}}\left({\displaystyle \frac{1}{2𝒩}}\right)^8\mathrm{}.`$ (111)
In the limit $`\xi 0`$ (i.e. $`𝒩\mathrm{}`$) for fixed $``$ we obtain the expression of specific heat of infinitely long strip with BK boundary condition of width $``$

$`C`$ $`=`$ $`{\displaystyle \frac{8}{\pi }}{\displaystyle \frac{+1}{}}\left(\mathrm{ln}(+1)+C_E+\mathrm{ln}{\displaystyle \frac{2^{3/2}}{\pi }}\right)2{\displaystyle \frac{\sqrt{2}}{}}`$ (112)
$``$ $`{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{2k_{2p}B_{2p}}{p(2p)!}}\left({\displaystyle \frac{\pi }{2}}\right)^{2p1}{\displaystyle \frac{1}{(+1)^{2p1}}}`$ (113)
$`=`$ $`{\displaystyle \frac{8}{\pi }}\left(1+{\displaystyle \frac{1}{}}\right)\mathrm{ln}+{\displaystyle \frac{8}{\pi }}\left(C_E+\mathrm{ln}{\displaystyle \frac{2^{3/2}}{\pi }}{\displaystyle \frac{\pi }{4}}\right)`$ (114)
$`+`$ $`{\displaystyle \frac{8}{\pi }}\left(C_E+\mathrm{ln}{\displaystyle \frac{2^{3/2}}{\pi }}+1{\displaystyle \frac{\sqrt{2}\pi }{8}}\right){\displaystyle \frac{1}{}}+\left({\displaystyle \frac{4}{\pi }}+{\displaystyle \frac{\pi }{18}}\right)\left({\displaystyle \frac{1}{}}\right)^2`$ (115)
$``$ $`\left({\displaystyle \frac{4}{3\pi }}+{\displaystyle \frac{\pi }{18}}\right)\left({\displaystyle \frac{1}{}}\right)^3+\left({\displaystyle \frac{2}{3\pi }}+{\displaystyle \frac{\pi }{18}}+{\displaystyle \frac{43\pi ^3}{21600}}\right)\left({\displaystyle \frac{1}{}}\right)^4`$ (116)
$``$ $`\left({\displaystyle \frac{2}{5\pi }}+{\displaystyle \frac{\pi }{18}}+{\displaystyle \frac{43\pi ^3}{7200}}\right)\left({\displaystyle \frac{1}{}}\right)^5+\mathrm{}.`$ (117)
Note that the specific heat expansion for infinitely long cylinder contains only even powers of $`𝒩^1`$ (except of course the leading logarithmic term), while in the specific heat expansion for infinitely long strip with BK boundary condition any integer powers of $`^1`$ can occur.

In Fig. 2 we plot the aspect-ratio $`(\xi )`$ dependence of the finite-size specific heat correction terms $`C_0,C_1,C_2,`$ and $`C_3`$ for the Ising model with BK boundary condition and those of the torus . We use the logarithmic scales for the horizontal axis. For large enough $`\xi `$ ($`1`$), the finite size properties of the Ising model with BK boundary condition and those of the torus become the same because the boundaries along the shorter direction determine the finite size properties of the system; for both BK boundary condition and the torus, the boundary condition along the $`y`$ axis is the periodic one.

## V Summary and discussion

In this paper, we have used the method of  to derive exact finite-size corrections for the free energy $`F`$ and the specific heat $`C`$ of the critical ferromagnetic Ising model on the $`\times 2𝒩`$ square lattice with Brascamp-Kunz (BK) boundary conditions . We find that the finite-size corrections to the free energy and the specific heat are always integer powers of $`𝒩^1`$ ($`^1`$) except of course the leading logarithmic term in the specific heat. In the finite-size expansion of the free energy given by Eq. (26), only even power of $`𝒩^1`$ occur, except for the term $`𝒩`$. In the finite-size expansion of the specific heat given by Eqs. (91) and (101), any integer powers of $`𝒩^1`$ ($`^1`$) can occur.

We have compared our results with those under toroidal boundary conditions. When the ratio $`\xi /2=(+1)/2𝒩`$ is smaller than 1 the behaviors of finite-size corrections for $`C`$ are quite different for BK and toroidal boundary conditions; when $`\mathrm{ln}(\xi /2)`$ is larger than 3, finite-size corrections for $`C`$ in two boundary conditions approach the same values. In the limit $`𝒩\mathrm{}`$ we obtain the expansion of the free energy for infinitely long strip with BK boundary conditions. Our results are consistent with the conformal field theory prediction for the mixed boundary conditions by Cardy  although the definitions of boundary conditions in two cases are different in one side of the long strip. It is of interest to know under what conditions different boundary conditions could still give the same finite-size corrections.

The results of this paper shows that the method of Ref.  is quite useful for calculating exact finite-size corrections for critical systems. It is of interest to apply this method to calculate exact finite-size corrections for the Ising model and other free models  on various lattices with various boundary conditions so that some general features of such finite-size corrections could be found.

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.

APPENDIX

## A Kronecker’s Double Series

Kronecker’s double series can be defined as

$$\mathrm{K}_p^{1/2,0}(\tau )=\frac{p!}{(2\pi i)^p}\underset{\genfrac{}{}{0.0pt}{}{m,nZ}{(m,n)(0,0)}}{}\frac{e^{\pi in}}{(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^{1/2,0}(\tau )=\frac{p!}{(2\pi i)^p}\underset{n0}{}\frac{e^{\pi in}}{n^p}\frac{p!}{(2\pi i)^p}\underset{m0}{}\underset{nZ}{}\frac{e^{\pi in}}{(n+\tau m)^p}.$$

Here the first sum gives nothing but Fourier representation of Bernoulli polynomials

$$\mathrm{B}_p^\alpha =\frac{p!}{(2\pi i)^p}\underset{n0}{}\frac{e^{2\pi in\alpha }}{n^p}.$$
(A1)
The second sum can be rearranged with the help of the identity

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

which can easily be derived from following equation

$$\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}$$
(A2)
by differentiating it $`p`$ times. The final result of our resummation of double Kronecker sum is

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

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

$$\mathrm{B}_{2p}^{1/2}\mathrm{K}_{2p}^{1/2,0}(i\xi )=4p\underset{m=1}{\overset{\mathrm{}}{}}\underset{n=0}{\overset{\mathrm{}}{}}(n+1/2)^{2p1}e^{2\pi m\xi (n+1/2)}.$$
(A3)
## B Relation between moments and cumulants

Moments $`Z_{2k}`$ and cumulants $`F_{2k}`$ which enters the expansion of exponent

$$\mathrm{exp}\left\{\underset{k=1}{\overset{\mathrm{}}{}}\frac{x^{2k}}{(2k)!}F_{2k}\right\}=1+\underset{k=1}{\overset{\mathrm{}}{}}\frac{x^{2k}}{(2k)!}Z_{2k}$$

are related to each other as

$`Z_2`$ $`=`$ $`F_2,`$ (B1)
$`Z_4`$ $`=`$ $`F_4+3F_2^2,`$ (B2)
$`Z_6`$ $`=`$ $`F_6+15F_2F_4+15F_2^3,`$ (B3)
$`Z_8`$ $`=`$ $`F_8+28F_2F_6+35F_4^2+210F_2^2F_4+105F_2^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 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]`$ (C1)
$`=`$ $`{\displaystyle \frac{1}{z^2}}+{\displaystyle \underset{p=2}{\overset{\mathrm{}}{}}}a_p(\tau )z^{2p2}.`$ (C2)
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),`$ (C3)
$`a_3`$ $`=`$ $`\frac{\pi ^6}{189}(\theta _2^4+\theta _3^4)(\theta _4^4\theta _2^4)(\theta _3^4+\theta _4^4),`$ (C4)
$`a_4`$ $`=`$ $`\frac{1}{3}a_2^2,`$ (C5)
$`a_5`$ $`=`$ $`\frac{3}{11}(a_2a_3),`$ (C6)
$`a_6`$ $`=`$ $`\frac{1}{39}(2a_2^3+3a_3^2),`$ (C7)
$`\mathrm{}`$ (C8)
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}^{1/2,0}(\tau )`$ 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^{\frac{1}{2},0}(\tau )=2^{1p}\mathrm{K}_p^{0,0}(\tau /2)\mathrm{K}_p^{0,0}(\tau ).`$ (C9)
Thus, Kronecker functions $`\mathrm{K}_{2p}^{1/2,0}(\tau )`$ 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,`$ (C10)
$`\theta _2^2(\tau /2)`$ $`=`$ $`2\theta _2\theta _3,`$ (C11)
$`2\theta _3^2(2\tau )`$ $`=`$ $`\theta _3^2+\theta _4^2,`$ (C12)
$`\theta _3^2(\tau /2)`$ $`=`$ $`\theta _2^2+\theta _3^2,`$ (C13)
$`2\theta _4^2(2\tau )`$ $`=`$ $`2\theta _3\theta _4,`$ (C14)
$`\theta _4^2(\tau /2)`$ $`=`$ $`\theta _3^2\theta _2^2.`$ (C15)
¿From the formulas above we can easily write down the Kronecker functions that have appeared in our asymptotic expansions

$`\mathrm{K}_4^{\frac{1}{2},0}(\tau )`$ $`=`$ $`\frac{1}{30}(\frac{7}{8}\theta _4^8\theta _2^4\theta _3^4),`$ (C16)
$`\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).`$ (C18)
Note that when $`\xi \mathrm{}`$ we have limits $`\theta _20`$, $`\theta _41`$, $`\theta _31`$ . The case $`\xi 0`$ can be obtained by using Jacobi’s imaginary transformation of the $`\theta `$ \- functions. In this case $`\theta _2\frac{1}{\sqrt{\xi }}`$, $`\theta _40`$ and $`\theta _3\frac{1}{\sqrt{\xi }}`$ and the Kronecker’s function can again be reduce to the Bernoulli polynomials.

## D Asymptotic Expansion of the Digamma Function $`\psi (N+\alpha )`$

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}}}`$ (D1)
$`=`$ $`\mathrm{ln}x{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}(1)^p{\displaystyle \frac{\mathrm{B}_p}{p}}{\displaystyle \frac{1}{x^p}}.`$ (D2)
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}}}`$ (D3)
$`=`$ $`\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}}`$ (D4)
$`=`$ $`\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}}.`$ (D5)
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},$$
(D6)
we finally obtain Eq. (76)

$$\psi (N+\alpha )=\mathrm{ln}N\underset{p=1}{\overset{\mathrm{}}{}}(1)^p\frac{\mathrm{B}_p^\alpha }{p}\frac{1}{N^p}.$$
(D7)
FIGURE CAPTIONS

FIG. 1. The Brascamp - Kunz boundary conditions for $`M\times 2N`$ lattice. Here $`M=7`$ and $`2N=8`$.

FIG. 2. Aspect-ratio $`(\xi )`$ dependence of finite-size correction terms for the specific heat of the square lattice Ising model with Brascamp-Kunz boundary conditions (solid lines) and toroidal boundary conditions (dashed lines): (a) $`C_0`$, (b) $`C_1`$, (c) $`C_2`$, and (d) $`C_3`$.


text:
# Quantization of adiabatic pumped charge in the presence of superconducting lead

## Acknowledgments

We gratefully acknowledge support by a RGC grant from the SAR Government of Hong Kong under grant number HKU 7215/99P.