text:
# Global solution of viscous accretion disk around rotating compact objects: a pseudo-general-relativistic study

## 1 Introduction

For over last three decades, fluid dynamical study of accretion disks around compact object is indeed a topic extensively discussed, mostly for black holes (Abramowicz et al. 1988, 1996; Abramowicz & Kato 1989; Chakrabarti 1990, 1996a,b,c; Narayan & Yi 1994; Narayan et al. 1997, 1998; Kato et al. 1998; Sponholz & Molteni 1994; Peitz & Appl 1997; Gammie & Popham 1998; Popham & Gammie 1998; Miwa et al. 1998; Lu & Yuan 1998; Manmoto 2000) and a few for neutron stars (Chakrabarti & Sahu 1997; Mukhopadhyay 2000, 2002a; Popham & Sunyaev 2001). Though in some of the cases authors use full general relativistic approach, in most of the cases the approaches are pseudo-Newtonian. In the pseudo-Newtonian method, they use simple non-relativistic equations to describe the disk with the appropriate choice of modified gravitational force which takes care about the geometry of back-ground space-time. There are a number of basic differences between the accretion disk around black hole and any other gravitating objects namely neutron star, white dwarf, strange star and other gravitating stars particularly at the inner edge of the disk. (1) While the black hole has event horizon, other gravitating objects mentioned above have hard surface. Thus, whereas in the case of other gravitating object close to the surface the matter speed must be subsonic, for black holes the matter speed is supersonic. (2) In case of other gravitating object that has hard surface, the matter comes into contact with its outer surface. Thus, there is an obvious effect of dragging to the exterior space-time due to its rotation. In terms of general relativity, corresponding space-time metric should be the solution for interior as well as exterior to it. Therefore, it is necessary to consider the ’dragging effect’ of rotating gravitating star (particularly for neutron star, where the relativistic effect is very important) to the space-time. In case of black hole, which does not have hard surface, this effect comes to the space-time automatically through the ’Lense-Thirring drag’ from the Kerr geometry which does not describe any physics interior to the horizon. (3) Further, the overall temperature as well as the luminosity of the disk is expected to be high in the cases of other gravitating star particularly for neutron star with respect to that of black hole (Mukhopadhyay 2000, 2002a) for a particular accretion rate. If the same parameters are chosen for the infalling matter towards black hole and other gravitating objects separately, at same accretion rate, the maximum temperature in the disk around black hole will be lower at least a few factors times. Thus, in the disk around neutron star, expected generation/absorption of nuclear energy is higher than that of black hole (Mukhopadhyay 2002a). As there is a clear distinction between accretion disks around black hole and other gravitating objects, hereinafter, we will call the class of gravitating objects having hard surface as GOHS to distinguish from black hole. GOHS can be meant as neutron star, strange star, white dwarf and any other highly gravitating objects where disk may form. As the radius of white dwarf is very large, the general relativity is not much important in the accretion disk around it. But our discussions will be generally true for disk around any kind of GOHS.

Mukhopadhyay (2002a) has studied the accretion disk around weakly magnetized slowly rotating neutron stars, and found that the temperature of the inner disk is significantly higher than the case of black holes. Prasanna & Mukhopadhyay (2003) have described the accretion disk with the inclusion of ’Coriolis’ effect. They have shown that the rotational effect of compact object to the accretion disk can be approximately reproduced with the use of Coriolis effect. They have done perturbative analysis and shown that the self-similar solutions are well behaved for most of the parameter space. As no cosmic object is static, in this paper, our interest is to include the rotation of the neutron star, globally speaking rotation of GOHS, in a more proper manner and see its effect onto the disk. The recent observation tells about the evidence of millisecond pulser (Strohmayer & Markwardt 2002) of frequency 582Hz for 4U 1636-53 which is the accreting LMXB. Therefore, it is important to consider the rotational effect of the GOHS to the disk. We approach our study in pseudo-Newtonian manner. So far, there is no pseudo-potential in the literature which can describe the accretion disk around GOHS. Thus, first of all we should establish such a potential useful for our work.

Shakura & Sunyaev (1973) first initiated the accretion disk modeling using simple Newtonian potential meant for non-rotating black holes. As the relativistic effects are extremely important near the compact objects, this potential can not describe the essential inner properties of the disk. Paczyński & Wiita (1980) modified this Newtonian potential in conformity with the Schwarzschild geometry, which can naturally reveal approximately the properties of disk around non-rotating black holes and GOHSs, even that of the inner most part of the disk. The potential can reproduce the radius of marginally stable orbit $`(x_s)`$ and marginally bound orbit $`(x_b)`$, mechanical energy per unit mass at the last stable circular orbit $`(E_s)`$ and viscous efficiency of the disk. The last two parameters agree within a $`10\%`$ error but the first two cases exactly match with that of the Schwarzschild metric. Nowak & Wagoner (1991) proposed another potential for accretion disk around non-rotating black holes. This potential can also mimic approximately most of the properties of the disk governed by Schwarzschild geometry. However, angular and epicyclic frequency of the disk can be best analyzed with this potential. Artemova et al. (1996) proposed a couple of pseudo-potentials to describe accretion disk for rotating and non-rotating black holes. These potentials are well analyzed by a number of astrophysicists later (e.g. Lovas 1998; Semerák & Karas 1999). Very recently, Mukhopadhyay (2002b) and Mukhopadhyay & Misra (2003) prescribed some new potentials with which properties of accretion disk can be very well described in Kerr geometry. The more interesting fact lies in the methodology adopted by Mukhopadhyay (2002b) to describe a potential in the accretion disk, which can be used to derive the pseudo-potential for any metric according to the physics concerned. The potentials proposed by Mukhopadhyay & Misra (2003) can be used for the time-dependent simulation of accretion disk. But all the potentials are described either for the fluid dynamical study of disk around black hole or the study of temporal effect around compact objects. As the fluid properties of inner accretion disk are very much influenced by the nature of central gravitating object as mentioned above, it is essential to establish a pseudo-potential which can describe an accretion disk around GOHS. Although in case of black hole, the space-time can be described by Kerr geometry, but in the case of GOHS with slow rotation, the corresponding metric which can continuously describe the space-time inside and outside the star is Hartle-Thorne (hereinafter HT, Hartle & Thorne 1968). Thus, we should establish our potential according to HT metric.

Thus, we shall consider the general solutions for fluid disks around rotating GOHS. In order to consider a solution procedure, we will follow Chakrabarti (1996b) where the global solution for transonic accretion flows is studied around rotating black holes but in the case of weak viscosity. Here, we will consider the set of disk equations in terms of general viscous prescription and the accretion around GOHS. As we want to concentrate upon the features of the inner region, we shall restrict the study to ’sub-Keplerian’ flow. We study, how does the rotation of GOHS affect on the global parameter space of accretion disk. We know that for the case of non-rotating compact object, the stable accretion disk forms only for a certain set of physical parameter. We would like to check how this parameter region gets affected, how does the steady disk structure change, with different values of the rotation to GOHS. Once we have the clear picture about the parameter space, we will concentrate upon the fluid properties of the disk for different viscosity and how does the viscous fluid properties change with different rotations of GOHS.

Therefore, we are going to present a complete set of work on accretion disk around GOHS. In order to do so, in the next section, we will derive the required pseudo-Newtonian potential and establish its efficiency. In §3, we will present the basic equations of viscous accretion disk around GOHS. In §4 and §5, we will analyse the parameter space of the accretion disk around GOHS and corresponding properties of accreting fluid respectively. Finally we shall present our discussion and overall summary in §6.

## 2 Pseudo-Newtonian potential and its efficiency

Here, following Mukhopadhyay (2002b), we will derive the pseudo-Newtonian potential for accretion disk modeling around GOHS according to HT metric. The Lagrangian density for a particle in the HT geometry (neglecting the higher quadruple terms) at the equatorial plane ($`\theta =\pi /2`$) can be written as

$`2=\left(1{\displaystyle \frac{2M}{r}}{\displaystyle \frac{2J^2}{r^4}}\right)\dot{t}^2+\left(1{\displaystyle \frac{2M}{r}}+{\displaystyle \frac{2J^2}{r^4}}\right)^1\dot{r}^2+r^2\dot{\varphi }^2{\displaystyle \frac{4J}{r}}\dot{\varphi }\dot{t},`$ (1)
where over-dots denote the derivative with respect to the proper-time $`\tau `$ and $`J`$ denotes the angular momentum of GOHS. As the metric is valid only for slowly rotating stars, the rotation parameter $`J`$ is restricted as $`J0.5`$. The geodesic equations of motion are

$`E=\mathrm{constant}=\left(1{\displaystyle \frac{2M}{r}}{\displaystyle \frac{2J^2}{r^4}}\right)\dot{t}+{\displaystyle \frac{2J}{r}}\dot{\varphi },`$ (2)
$`\lambda =\mathrm{constant}=r^2\dot{\varphi }{\displaystyle \frac{2J}{r}}\dot{t}.`$ (3)
For the particle with non-zero rest mass, $`g_{\mu \nu }p^\mu p^\nu =m^2`$ (where $`p^\mu `$ is the momentum of the particles and $`g_{\mu \nu }`$ is the metric). Replacing the solution for $`\dot{t}`$ and $`\dot{\varphi }`$ from (2) and (3) into (1) gives a differential equation for $`r`$

$`\left({\displaystyle \frac{dr}{d\tau }}\right)^2`$ $`=`$ $`\left({\displaystyle \frac{4h(r)J^2}{g(r)r^6}}+{\displaystyle \frac{16J^4}{g(r)r^{10}}}{\displaystyle \frac{g(r)}{r^2}}\right)\lambda ^2+\left({\displaystyle \frac{h(r)}{g(r)}}+{\displaystyle \frac{4J^2}{g(r)r^4}}\right)E^2`$ (4)
$`+`$ $`\left({\displaystyle \frac{4h(r)J}{g(r)r^3}}{\displaystyle \frac{16J^3}{g(r)r^7}}\right)E\lambda g(r)m^2=\mathrm{\Psi },`$
where, $`g(r)=\left(1\frac{2m}{r}+\frac{2J^2}{r^4}\right)`$, $`h(r)=\left(1\frac{2m}{r}\frac{2J^2}{r^4}\right)`$. Here, $`\mathrm{\Psi }`$ can be identified as an effective potential for the radial geodesic motion. The conditions for circular orbits are

$`\mathrm{\Psi }=0,{\displaystyle \frac{d\mathrm{\Psi }}{dr}}=0.`$ (5)
Solving for $`E`$ and $`\lambda `$ from (5) we get

$`E=X(r),`$ (6)
$`\lambda `$ $`=`$ $`X(r)\sqrt{{\displaystyle \frac{m^4r^6(4J^2+Mr^3)^2}{[2J^3m^2+Jm^2r^3(4m+3r)+\sqrt{m^4(5J^2+Mr^3)(2J^2+r^3(2M+r))^2}]^2}}},`$
where,

$`X(r)`$ $`=`$ $`{\displaystyle \frac{1}{m^2r^3(4J^2+Mr^3)}}[2J^3m^2+Jm^2(4M3r)r^3`$ (8)
$`\sqrt{m^4(5J^2+Mr^3)(2J^2+r^3(2M+r))^2}]`$
$`\{(12J^4m^2r^2+2J^2m^2(9M7r)r^5+m^2M(3Mr)r^8`$
$`+`$ $`6Jr^2\sqrt{m^4(5J^2+Mr^3)(2J^2+r^3(2M+r))^2})/`$
$`(36J^4r^6(3M+r)^2+12J^2r^3(3M+2r))\}^{1/2}.`$
Now as standard practice, we can define the Keplerian angular momentum distribution, $`\lambda _K=\frac{\lambda }{E}`$. Therefore, corresponding centrifugal force in HT geometry can be written as

$`{\displaystyle \frac{\lambda _K^2}{r^3}}={\displaystyle \frac{m^4r^3(4J^2+Mr^3)^2}{[2J^3m^2+Jm^2r^3(4m+3r)+\sqrt{m^4(5J^2+Mr^3)(2J^2+r^3(2M+r))^2}]^2}}=F_r.`$ (9)
Thus from above, $`F_r`$ can be identified as the gravitational force of rotating GOHS at the Keplerian orbit in an equatorial plane. If we choose $`x=r/M`$ (as $`G=c=1`$) and $`m=1`$, above expression reduces to

$`F={\displaystyle \frac{x^3(4J^2+x^3)^2}{[2J^3+Jx^3(4+3x)+\sqrt{(5J^2+x^3)(2J^2+x^3(2+x))^2}]^2}}.`$ (10)
Also the above expression reduces to Paczyński-Wiita (1980) form for $`J=0`$. Thus we propose that (10) is the most general form of the gravitational force corresponding to the pseudo-potential in accretion disk around GOHS at equatorial plane.

### 2.1 Comparison of the Results for Hartle-Thorne geometry and Pseudo-Potential

Now we will establish the validity of this potential. We would like to check, whether this potential can reproduce the values of marginally bound ($`x_b`$), marginally stable ($`x_s`$) orbits and mechanical energy per unit mass at $`x_s`$ ($`E_s`$), as same as HT geometry or not.

In terms of the peudo-potential, one can write the radial velocity ($`v`$), angular velocity ($`\mathrm{\Omega }`$) and angular momentum per unit mass ($`\lambda `$) for a Keplerian orbit in an accretion disk as (Lynden-Bell 1969, Pringle & Rees 1973, Shakura & Sunyaev 1973, Novikov & Thorne 1973, Paczyński & Wiita 1980),

$`v=\sqrt{x{\displaystyle \frac{dV}{dx}}},\mathrm{\Omega }=\sqrt{{\displaystyle \frac{1}{x}}{\displaystyle \frac{dV}{dx}}},\lambda =\sqrt{x^3{\displaystyle \frac{dV}{dx}}}.`$ (11)
At the marginally bound orbit, mechanical energy $`E`$ reduces to zero and we get

$`{\displaystyle \frac{v^2}{2}}+V={\displaystyle \frac{x}{2}}{\displaystyle \frac{dV}{dx}}+V=0,`$ (12)
where, $`V=F𝑑x`$, which can be used to calculate $`x_b`$ for the potential. For the last stable orbit ($`x_s`$), $`d\lambda /dx=0`$, for our pseudo-potential which gives

$`x^5[8J^4+2J^2(52x)x^3+(x2)x^6][240J^6+(x6)x^9+8J^4x^3(12+5x)`$
$`+2J^2x^6(20x33)+(24J^3+12Jx^3)\sqrt{(5J^2+x^3)(2J^2+(x2)x^3)^2}]=0.`$ (13)
The solution of (13) gives the location of the last stable circular orbit ($`x_s`$). For any $`J`$, the value of $`x_s`$ computed from the above equation (13), matches exactly with the radius of last stable circular orbit in HT geometry. In Table-1, 2 and 3, we list $`x_b`$, $`x_s`$ and $`E_s`$ respectively for the potential $`V`$ and HT geometry for various values of $`J`$.

Table-1

Values of $`x_b`$

| $`J`$ | $`0`$ | $`0.1`$ | $`0.2`$ | $`0.3`$ | $`0.4`$ | $`0.5`$ |
| --- | --- | --- | --- | --- | --- | --- |
| $`V`$ | $`4`$ | $`3.7862`$ | $`3.5619`$ | $`3.3255`$ | $`3.0762`$ | $`2.8135`$ |
| HT | $`4`$ | $`3.7954`$ | $`3.5810`$ | $`3.3557`$ | $`3.118`$ | $`2.8689`$ |
| $`J`$ | $`0`$ | $`0.1`$ | $`0.2`$ | $`0.3`$ | $`0.4`$ | $`0.5`$ |
| $`V`$ | $`4`$ | $`4.2042`$ | $`4.4`$ | $`4.5885`$ | $`4.77`$ | $`4.9462`$ |
| HT | $`4`$ | $`4.1957`$ | $`4.3838`$ | $`4.5684`$ | $`4.7397`$ | $`4.9089`$ |

Table-2

Values of $`x_s`$

| $`J`$ | $`0`$ | $`0.1`$ | $`0.2`$ | $`0.3`$ | $`0.4`$ | $`0.5`$ |
| --- | --- | --- | --- | --- | --- | --- |
| $`V`$ & HT | $`6.0`$ | $`5.6648`$ | $`5.3107`$ | $`4.9339`$ | $`4.5299`$ | $`4.0934`$ |
| $`J`$ | $`0`$ | $`0.1`$ | $`0.2`$ | $`0.3`$ | $`0.4`$ | $`0.5`$ |
| $`V`$ & HT | $`6.0`$ | $`6.3189`$ | $`6.6240`$ | $`6.9170`$ | $`7.1994`$ | $`7.4724`$ |

Table-3

Values of $`E_s`$

| $`J`$ | $`0`$ | $`0.1`$ | $`0.2`$ | $`0.3`$ | $`0.4`$ | $`0.5`$ |
| --- | --- | --- | --- | --- | --- | --- |
| $`V`$ | $`0.0625`$ | $`0.0664`$ | $`0.0711`$ | $`0.0769`$ | $`0.0843`$ | $`0.0943`$ |
| HT | $`0.0572`$ | $`0.0607`$ | $`0.0649`$ | $`0.0701`$ | $`0.0767`$ | $`0.0857`$ |
| $`J`$ | $`0`$ | $`0.1`$ | $`0.2`$ | $`0.3`$ | $`0.4`$ | $`0.5`$ |
| $`V`$ | $`0.0625`$ | $`0.0592`$ | $`0.0563`$ | $`0.0538`$ | $`0.0516`$ | $`0.0497`$ |
| HT | $`0.0572`$ | $`0.0542`$ | $`0.0517`$ | $`0.0494`$ | $`0.0474`$ | $`0.0457`$ |

From Table-1, it is clear that for all values of $`J`$, $`V`$ can reproduce the value of $`x_b`$ in very good agreement with general relativistic (HT) results. The maximum error in $`x_b`$ is $`4\%`$. Table-2 indicates that $`V`$ reproduces $`x_s`$ exactly as that of HT geometry. Table-3 shows a maximum possible error in $`E_s`$ is $`10\%`$. Thus, the potential $`V`$ will produce a slightly larger luminosity than the general relativistic one in the accretion disk for a particular value of $`J`$. It is to be noticed that for counter-rotating GOHS, that is for retrograde orbits, the errors are less than those of co-rotating ones. Eventually, we can tell that our potential can describe approximately all the phenomena of that of HT geometry and can claim it as a good pseudo-potential to describe the relativistic accretion disk around the rotating GOHS, particularly close to the equatorial plane. Thus, for a thin disk or vertically averaged disk, this is an ideal pseudo-potential.

## 3 Basic equations of accretion disk and their consequences

Here, throughout our calculations, we express the radial coordinate in unit of $`GM/c^2`$, where $`M`$ is mass of the compact object, $`G`$ is the gravitational constant and $`c`$ is the speed of light. We also express the velocity in unit of speed of light and the angular momentum in unit of $`GM/c`$. The equations to be solved are given below.
(1) Equation of continuity:

$$4\pi x\mathrm{\Sigma }v=\dot{M},$$
(14)
where $`\mathrm{\Sigma }`$ is the vertically integrated density and $`\dot{M}`$ is the accretion rate. Following Matsumoto et al. (1984), we can calculate the $`\mathrm{\Sigma }`$ as

$$\mathrm{\Sigma }=I_n\rho _eh(x),$$
(15)
where $`\rho _e`$ is the density at equatorial plane, $`h(x)`$ is the half-thickness of the disk and $`I_n=\frac{(2^nn!)^2}{(2n+1)!}`$, where $`n`$ is the polytropic index. From the vertical equilibrium assumption, the half-thickness can be written as

$$h(x)=c_sx^{1/2}F^{1/2},$$
(16)
where $`c_s`$ is the speed of sound.
(2) Radial momentum balance equation:

$$v\frac{dv}{dx}+\frac{1}{\rho }\frac{dP}{dx}\frac{\lambda ^2}{x^3}+F(x)=0,$$
(17)
where the flow is chosen adiabatic with the equation of state to be $`P=k\rho ^\gamma `$ and $`c_s^2=\frac{\gamma P}{\rho }`$.
(3) Azimuthal momentum balance equation: Here, we will follow Chakrabarti (1996a) to express the viscous dissipation $`Q^+`$ in terms of shear stress $`W_{x\varphi }`$. Because of this viscous dissipation, the angular momentum varies in accretion disk. Thus, $`W_{x\varphi }=\alpha (I_{n+1}P+I_nv^2\rho )h(x)`$ and $`Q^+=\frac{W_{x\varphi }^2}{\eta }`$, $`\eta `$ is coefficient of viscosity and $`\alpha `$ is Shakura-Sunyaev (1973) viscosity parameter. Thus the equation to be

$$v\frac{d\lambda }{dx}=\frac{1}{\rho h(x)x}\frac{d}{dx}\left[x^2\alpha \left(\frac{I_{n+1}}{I_n}P+v^2\rho \right)h(x)\right].$$
(18)
(4) Entropy equation: According to the mixed shear stress (Chakrabarti 1996a), $`Q^+=\alpha (I_{n+1}P+I_nv^2\rho )h(x)x\frac{d\mathrm{\Omega }}{dx}`$. For simplicity, we also consider the heat lost is proportional to the heat gained by the flow. Thus the equation to be

$$\mathrm{\Sigma }vT\frac{ds}{dx}=\frac{vh(x)}{\mathrm{\Gamma }_31}\left(\frac{dP}{dx}\mathrm{\Gamma }_1\frac{P}{\rho }\frac{d\rho }{dx}\right)=Q^+Q^{}=fQ^+,$$
(19)
where $`s`$ is entropy density and $`f`$ is cooling factor which is close to $`0`$ and $`1`$ for the flow with efficient and inefficient cooling respectively. Following Cox & Giuli (1968), we can define

$$\mathrm{\Gamma }_3=1+\frac{\mathrm{\Gamma }_1\beta }{43\beta },\mathrm{\Gamma }_1=\beta +\frac{(43\beta )^2(\gamma 1)}{\beta +12(\gamma 1)(1\beta )},\beta =\frac{\rho kT/\mu m_p}{\overline{a}T^4/3+\rho kT/\mu m_p}.$$

Here, $`\beta `$ (ratio of gas pressure to total pressure) close to $`0`$ for radiation dominated flow (highly relativistic flow) and close to $`1`$ for gas dominated flow. Now combining (14),(17), (18) and (19), we get

$$\frac{dv}{dx}=\frac{f_1(x,v,c_s)}{f_2(v,c_s)},$$
(20)
where,

$`f_1(x,v,c_s)`$ $`=`$ $`\left[c_s^2\left({\displaystyle \frac{3}{2x}}{\displaystyle \frac{1}{2F}}{\displaystyle \frac{dF}{dx}}\right)\gamma \left(F{\displaystyle \frac{\lambda ^2}{x^3}}\right)\right]\left[v(\mathrm{\Gamma }_1+1){\displaystyle \frac{4\alpha \mathrm{\Lambda }}{v(3\gamma 1)}}\right]`$ (21)
$``$ $`\left[{\displaystyle \frac{2\alpha c_s^2}{(3\gamma 1)xv}}+{\displaystyle \frac{\alpha v}{x}}{\displaystyle \frac{2\lambda }{x^2}}\right]\mathrm{\Lambda }+\left({\displaystyle \frac{3}{2x}}{\displaystyle \frac{1}{2F}}{\displaystyle \frac{dF}{dx}}\right)vc_s^2(\mathrm{\Gamma }_11),`$
$`f_2(v,c_s)=\left[1{\displaystyle \frac{2c_s^2}{(3\gamma 1)v^2}}\right]\alpha \mathrm{\Lambda }c_s^2(\mathrm{\Gamma }_11)+\left(\gamma v{\displaystyle \frac{c_s^2}{v}}\right)\left[v(\mathrm{\Gamma }_1+1){\displaystyle \frac{4\alpha \mathrm{\Lambda }}{(3\gamma 1)v}}\right]`$ (22)
and $`\mathrm{\Lambda }=\gamma (\mathrm{\Gamma }_31)f\alpha \left(\frac{I_{n+1}}{I_n}\frac{c_s^2}{\gamma }+v^2\right)`$.
As far away from the black hole, $`v<c_s`$ and close to it, $`v>c_s`$, there must be an intermediate location, where the denominator of (20) must be vanished. To have a smooth solution, at that location, numerator has to be zero. This location is called the sonic point or critical point ($`x_c`$). The existence of sonic location plays an important role in accretion phenomena. Although, in case of accretion flow of matter around black hole sonic point must exist, for that around GOHS (say, neutron star) it is not always necessary (Mukhopadhyay 2002a). When shock forms in an accretion disk around a neutron star (Chakrabarti & Sahu 1997, Mukhopadhyay 2002a), sonic points play an important role. From the global study of the sonic points in an accretion disk, one can understand about the stability of physical parameter region, which we will discuss in the next section.

Now, when the sonic point exists in an accretion disk around GOHS, $`f_2(v_c,c_{sc})=0`$ at $`x=x_c`$. Thus at that location Mach number can be written as

$`M_c^2={\displaystyle \frac{+\sqrt{^24𝒜𝒞}}{2𝒜}},`$ (23)
where

$`𝒜`$ $`=`$ $`\alpha ^2\gamma f(\mathrm{\Gamma }_31)+\gamma (\mathrm{\Gamma }_1+1){\displaystyle \frac{4\alpha ^2f(\mathrm{\Gamma }_31)\gamma ^2}{3\gamma 1}},`$
$``$ $`=`$ $`2\mathrm{\Gamma }_1{\displaystyle \frac{4\alpha ^2\gamma f(\mathrm{\Gamma }_31)}{3\gamma 1}}\left(1{\displaystyle \frac{I_{n+1}}{I_n}}\right),`$
$`𝒞`$ $`=`$ $`{\displaystyle \frac{2\alpha ^2f(\mathrm{\Gamma }_31)}{3\gamma 1}}{\displaystyle \frac{I_{n+1}}{I_n}}.`$ (24)
Subsequently, at $`x=x_c`$, $`f_1(x_c,v_c,c_{sc})=0`$. Thus, using (23) we can eliminate $`v_c`$ from $`f_1(x_c,v_c,c_{sc})=0`$ and get an algebraic equation for $`c_{sc}`$ and $`x_c`$, which we can solve to find out sound speed at sonic location. To find out $`\frac{dv}{dx}|_c`$, we apply l’Hospital’s rule to (20). In order to understand the fluid properties in accretion disk, we have to solve (20) with an appropriate boundary condition. Also, integrating (18), we get the angular momentum of the accreting matter as

$`\lambda ={\displaystyle \frac{\alpha xc_s^2}{v}}\left({\displaystyle \frac{2}{3\gamma 1}}+M^2\right)+\lambda _{in},`$ (25)
where $`M`$ is the Mach number of the flow and $`\lambda _{in}`$ is the angular momentum at the inner edge of accretion disk.

Now integrating (14) and (17), we can write the entropy and energy of the flow at sonic point as

$$\dot{}_c=x_c^{3/2}F_c^{1/2}(\gamma +1)^{q/2}\left(\frac{\frac{\lambda ^2}{x_c^3}F_c}{\frac{1}{F_c}\frac{dF}{dx}|_c\frac{3}{x_c}}\right)^{\frac{\gamma }{\gamma 1}}$$
(26)
and

$$E_c=\frac{2\gamma }{(\gamma 1)}\left(\frac{\frac{\lambda ^2}{x_c^3}F_c}{\frac{1}{F_c}\frac{dF_c}{dx}|_c\frac{3}{x_c}}\right)+V_c+\frac{\lambda ^2}{2x_c^2},$$
(27)
where $`V_c=(F𝑑x)|_c`$ and $`q=\frac{(\gamma +1)}{2(\gamma 1)}`$. It can be mentioned that, disk entropy and accretion rate are related by a simple relation as $`\dot{}=(\gamma K)^n\dot{M}`$. While $`\dot{M}`$ is a conserved quantity for the particular accretion flow, $`\dot{}`$ is not, as it contains $`K^n`$ which carries the entropy information that is not conserved in a dissipative system. As a boundary condition for the particular flow, we have to supply the sonic energy $`E_c`$. Then from (27), we can find out the sonic location $`x_c`$. Therefore, knowing $`x_c`$ one can easily find out the fluid velocity and sound speed at the sonic point from $`f_1=f_2=0`$ with the help of (23) for a particular accretion flow. These have to be supplied as further boundary conditions of the flow.

Now let us come to the issue of the formation of shock. Mukhopadhyay (2002a) showed that in accretion disk around neutron star, double shock is very natural in certain physical situations. Here, whenever we mention the shock, we will mean Rankine-Hugoniot shock (Landau & Lifshitz 1987). If we generalise the conditions to form a shock in an accretion disk given by Chakrabarti (1989) for rotating GOHS, we get

$$\frac{1}{2}M_+^2c_{s+}^2+nc_{s+}^2=\frac{1}{2}M_{}^2c_s^2+nc_s^2,$$
(28)
$$\frac{c_{s+}^\nu }{\dot{}_+}\left(\frac{2\gamma }{3\gamma 1}+\gamma M_+^2\right)=\frac{c_s^\nu }{\dot{}_{}}\left(\frac{2\gamma }{3\gamma 1}+\gamma M_{}^2\right),$$
(29)
$$\dot{}_+>\dot{}_{},$$
(30)
where

$$\dot{}=Mc_s^{2(n+1)}\frac{x_s^{3/2}}{\sqrt{F_s}}.$$
(31)
Here, subscript ’$``$’ and ’$`+`$’ indicate the quantities just before and after the shock respectively, $`x_s`$ indicates the shock location and $`\nu =\frac{3\gamma 1}{\gamma 1}`$. From (28)-(31), it is very clear that except the entropy expression, all remain unchanged with respect to the case of a non-rotating central object. Also from (29) and (31), we get the shock invariant quantity as

$$C=\frac{\left(\frac{2}{M_+}+(3\gamma 1)M_+\right)^2}{M_+^2(\gamma 1)+2}=\frac{\left(\frac{2}{M_{}}+(3\gamma 1)M_{}\right)^2}{M_{}^2(\gamma 1)+2},$$
(32)
which remains unchanged with respect to a non-rotating case. If only all the conditions, (28)-(30) and (32) are simultaneously satisfied by the matter, shock will be generated in an accretion disk.

## 4 Global analysis of the parameter space

One of our aim is to check, how does the rotation of GOHS affect the disk parameter region known for non-rotating case. More precisely, we will check, how does the rotation of GOHS affect the sonic location, structure as well as the stability of disk. As mentioned in §1, we will follow the methodology adopted by Chakrabarti (1996b) for global solution and Mukhopadhyay & Chakrabarti (2001) for stability analysis. In Fig. 1, we show the variation of disk entropy as a function of sonic location. The intersection to the horizontal line (which indicates a constant entropy curve) with all the curves indicate the sonic points of the accretion disk for that particular entropy and rotation to the GOHS. It is clearly seen that at a particular entropy, if the co-rotation of GOHS increases, the inner sonic points shift to more inner edge of the disk, which are comparatively unstable region (because, beyond of $`x_s`$, disk orbit is not stable and decrement of the value of $`x_s`$ is slower than that of $`x_c`$ with $`J`$). For $`J=0.5`$, the inner sonic point disappears, as we see that dashed curve terminates on the radial axis at $`x_c6`$ and it does not turn up to give rise the inner sonic point. The reason behind it is that, for higher the co-rotation, the angular momentum of the system will become higher and radial matter speed will be unable to overcome the centrifugal barrier close to the GOHS and that will make the disk unstable at the inner edge. For the counter-rotating cases, if $`J`$ increases in magnitude for a particular entropy, disk will become more stable, as the radial speed will be able to overcome the centrifugal barrier at a larger radius compared to a co-rotating case and matter may attain supersonic speed at a larger distance, which is more stable region of the disk. However, for a high counter rotation, angular momentum of the system will become too low and the accretion will tend to ’Bondi-like flow’ (Bondi 1952) with the existence of single sonic point.

It is seen that for $`\alpha =00.1`$, the variation of curves are almost same and only significant changes come into the picture for $`\alpha >0.1`$. Also for the lower viscosity, if the co-rotation is high ($`J`$ is high), the possibility to have the all three sonic points in the disk is high. As the viscosity increases, the possibility for the existence of all three sonic points diminishes. In order to preserve the existence of sonic point for the high viscous cases, $`J`$ must reduce to zero and subsequently, more negative. As a whole, from Fig. 1a-c we can conclude that for lower $`\alpha `$, multiple sonic points may exist, but for high $`\alpha `$, there is only one sonic point which is inner one in the accretion disk. Similar features come out from Fig. 2 but in terms of energy of the sonic point in accretion disk. As we know that the observed angular frequency for the candidate 4U 1636-53 is 582Hz, which is equivalent to $`J=0.2877`$, one of the value of $`J`$ chosen for our discussion is $`0.2877`$. For both Fig. 1 and 2, sonic points with negative slope of the curve indicate the locations of ’saddle-type’ or ’X-type’ sonic point and positive slopes indicate the ’centre-type’ or ’O-type’ sonic point. Thus, the rotation of GOHS and viscosity of the disk are the important factors to the formation and location of sonic points, which are related to the structure of accretion disk.

Figure 3 indicates the variation of energy with entropy at sonic locations, when $`J`$ is parameter for different viscosity of the fluid. It is seen that, with the increase of $`J`$ (co-rotation), inner sonic points of the disk shift towards the lower entropy region keeping the outer sonic points unchanged and the disk becomes unstable. We know that the shock in accretion disk can be formed only if there is a possibility of matter in the outer sonic point branch of lower entropy to jump to the inner sonic point branch of higher entropy. As because, with the increase of $`J`$, inner sonic point branch itself shifts to the lower entropy region which is comparatively unstable, the shock as well as disk also become unstable. For $`J=0.5`$, the inner sonic point branch disappears totally and only outer and O-type sonic point branches remain. Clearly, for a highly rotating GOHS, the possibility of formation of shock reduces as well as the shock becomes unstable, as the inner sonic point branch itself disappears or tends to disappear. The physical reason behind it is that, with the increase of value of $`J`$, the angular momentum of the system increases, which helps the disk to maintain a high azimuthal speed of matter upto a very close radius to the GOHS, as a result the radial speed of the matter can overcome the corresponding sound speed at a very inner radius only. On the otherhand, the inner edge of the accretion disk is comparatively less stable, as the entropy decreases at lower radii. As the sonic points form at a more inner radii, with the increase of $`J`$, disk tends to an unstable situation.

In case of the counter-rotation of GOHS, as $`J`$ increases in magnitude, the inner sonic point branch shifts towards a higher entropy region, and the disk as well as the shock become more stable. The physical reason to have a more stable inner sonic point branch for the counter-rotating case is the following. As the counter rotation of central object increases, that is the angular momentum of the system reduces, the matter falls freely towards GOHS at larger radii and it becomes supersonic at a comparatively outer edge of the disk. Thus, the disk as well as shock become more stable. However, for a very high counter-rotation of GOHS, this inner sonic point branch itself tends to merge to the outer sonic point branch, and the entire system shift towards the outer side which is more stable. As in case of higher counter-rotation, the angular momentum of the system becomes very small, the corresponding centrifugal pressure onto the matter becomes insignificant, as a result the possibility of shock diminishes again. Thus, if the angular momentum of GOHS increases or decreases significantly, the shock wave in accretion disk becomes unstable and the disk itself tends to an unstable situation. In Fig. 3, ’C’ is the point, where the flow can pass through the inner and outer sonic point simultaneously. The parameters around C are important to study the shock in accretion disk.

As the viscosity of flow increases, the possibility of intersection between inner and outer sonic point branch decreases. With the increase of viscosity, more regions containing the X-type sonic points become O-type. The outer X-type sonic points recedes in further out and inner one proceeds in further inwards. Also for higher viscosity, the outer sonic points may no longer remain as X-type and only the inner sonic points may exist. Thus, the possibility of shock will reduce.

## 5 Fluid Properties of Accretion Disk

Now we will discuss the behaviour of viscous accreting fluid around rotating GOHS. In the previous section, we have already confirmed that the rotation significantly affects the global parameter space of accretion disk. Here, our intention is to see, how the various fluid dynamical results of an accretion disk are affected for various values of $`J`$ and $`\alpha `$. For a certain choice of physical parameter set, if the existence of sonic radii are possible in the accretion disk around GOHS, the flow structures can be studied for the corresponding choice of sonic radii (or energy) of the accretion flow.

In Fig. 4a, we show the variation of Mach number for three different values of $`J`$ at a particular viscosity. As the viscosity is low, the rate of energy momentum transfer in matter is less. It is known that for an accretion around neutron star, in a certain cases, matter may be always subsonic (Chakrabarti & Sahu 1997; Mukhopadhyay 2002a) and at the stellar surface matter speed reduces to zero. Owing to this fact, here such a situation is considered where the speed of the accreting fluid is low, which accounts for the high residence time of matter in the disk. The angular momentum of the disk is high at $`x10`$ (see Fig. 4d) and thus radial matter speed is low. Subsequently, though the disk angular momentum goes down, as the matter comes closer to the stellar surface, it starts to decelerate and the radial matter speed still remains low. This results to a higher the possibility of cooling, and thus cooling factor $`f`$ is chosen as $`0.1`$. As the energy momentum transfer rate of matter is low, the accretion rate is chosen intermediate. We have chosen a standard mass of GOHS as $`2M_{}`$. If the GOHS is chosen to be non-rotating, angular momentum of the system is less and the centrifugal barrier is minimum. With the increase of $`J`$, this centrifugal barrier increases, and it becomes very high particularly for $`J=0.5`$. Mach number profiles clearly indicate that at the surface of GOHS, accreting fluid stops. Also from Fig. 4b and 4c, it comes out that sudden deceleration of the accreting fluid gives rise to a sudden enhancement of density and temperature close to GOHS. Here, the temperature is thought to be cooled down by a factor of $`1/30`$ through the inverse-Compton effect in the radiation pressure dominated relativistic flow. Thus, we choose $`\beta 0.03`$. Nevertheless, the temperature is still very high, which was also reported earlier (Mukhopadhyay 2002a) in the context of non-rotating neutron star. The density profiles are plotted in unit of $`\frac{c^6}{G^3M^2}`$. Figure 4d shows the ratio of sub-Keplerian to Keplerian angular momentum variation of the accreting fluid as a function of radial coordinate. It is seen that, at the boundary between Keplerian and sub-Keplerian disk, $`\lambda _k/\lambda 1`$, and from that radius ($`x_K`$) we start our calculation as our discussion is for sub-Keplerian accretion disk. From the $`\lambda _k/\lambda `$ profile, it is again clear that for higher $`J`$, centrifugal barrier becomes stronger. Far away from GOHS, the effect of centrifugal barrier is weak, but as the matter approaches towards GOHS, it starts to dominate and at $`x10`$ it becomes stronger. Then again it goes off. However, at very close to the GOHS, again this centrifugal effect comes into the picture as the matter falls into the strong rotational field of GOHS. This effect appears very prominently for high value of $`J`$, which is here chosen to be $`0.5`$. Actually, this effect comes into the picture, once the matter enters into the corresponding ergosphere of HT geometry. But the important point to be noted that if the outer radius of GOHS is greater than the radius of ergosphere, this feature is not possible to occur.

In Fig. 5, we show the results of accretion phenomena for 4U 1636-53, whose angular momentum $`J=0.2877`$. The examples are shown, where the sonic point(s) exist in the flow and shock forms in the accretion disk around 4U 1636-53 for different viscosity parameters. Figure 5a shows the variation of Mach numbers as a function of radial coordinate. It reflects that, for a small $`\alpha `$ ($`0`$) only shock is possible at $`x=19.02`$ in the accretion disk. As $`\alpha `$ increases, the energy momentum transfer rate increases, rate of infalling the matter enhances and velocity of the accreting fluid becomes high at the inner edge of the disk. But, because of the inner hard surface, matter has to stop at close to GOHS and naturally another shock forms at the inner edge of the accretion disk. The shock locations for $`\alpha =0.02`$ and $`0.05`$, are at $`x=14.27,4.38`$ and $`x=15.37,4.31`$ respectively. If the viscosity is high, the cooling factor is intermediate (see Fig. 5 caption). Also, the angular momentum of matter at the sonic point is chosen as, $`\lambda _c=3`$, and for the radiation dominated, inverse-Comptonised flow, $`\beta `$ is chosen to be $`0.03`$. Figures 5b and 5c show the variations of density and temperature. The regions, where velocity decelerates abruptly, temperature and density also jumps up in a significant manner. Here also, the inverse-Comptonised temperature is very high as the cases of Fig. 4.

## 6 Summary and Discussion

Here, we have studied the viscous accretion phenomena around rotating gravitational object with hard surface (GOHS), i.e., mostly around compact object like neutron star and strange star. We started with a proposal of a new pseudo-Newtonian potential by means of Hartle-Thorne geometry, as this metric can describe the space-time of interior as well as exterior to a slowly rotating star. We have shown that this potential can describe all the essential relativistic properties like, radius of marginally bound and stable orbit, specific mechanical energy etc., of accretion disk within $`10\%`$ error. Then, along with this potential, we have described the accretion disk by using non-relativistic equations. Therefore, we can say, our method is a ’pseudo-general-relativistic’. Once the pseudo-Newtonian potential is proposed, we have applied it to the study of relativistic fluid properties of accretion disk around rotating GOHS. We have prescribed the set of basic equations for vertically averaged, thin, viscous accretion disk. Then, we have analysed the viscous parameter space globally. Subsequently, we have discussed about the properties of viscous fluid and the effect of rotation on the fluid behaviour. Examples are shown for the accretion disk around an observed candidate 4U 1636-53, which is one of the fast rotating compact object with angular frequency 582Hz (equivalent $`J=0.2877`$). As starting from a new pseudo-potential, we have discussed upto the solution of viscous accretion flow, this is one of the most self consistent paper on accretion disk.

We have studied the stability of accretion disk. Most of the earlier analysis of structure and stability of the accretion disk have been done around non-rotating central object. As no cosmic object is static, it is very important to consider the effect of rotation, particularly for the discussion of an inner edge of the disk. We have found that, when the rotation to a GOHS is incorporated, the locations of sonic point (if any) get shifted to the unstable region and the shock (if any) gets affected in the disk, and thus the disk structure gets influenced. Therefore, we can conclude that the rotation should be incorporated in the study of an accretion disk. Any related conclusion is dependent on the rotation parameter of GOHS. The mentioned disk properties and phenomena not only get modified and/or shifted in location, sometimes disappeared completely. Prasanna & Mukhopadhyay (2003), worked on the stability analysis of accretion disk around rotating compact objects in a perturbative manner. They incorporated the rotational effect of central object in a disk indirectly by the inclusion of Coriolis acceleration term. Here, this rotational effect is brought directly from the metric itself.

From the global analysis of sonic points, we have seen that for a higher co-rotation of GOHS, the disk becomes unstable for a particular angular momentum of the accreting matter. As the formation of shock needs a stable inner sonic point in the accretion disk, for the higher value of $`J`$, shock is unstable as the inner region of disk is unstable. On the other hand, for counter-rotating cases, the angular momentum of system reduces and the matter falls more strongly to a GOHS. As there is no significant centrifugal barrier to slow down the matter in disk, the possibility of shock reduces again. Also the branch of inner sonic point merges or tends to merge to that of outer sonic point. Thus, we can conclude that the parameter region, where the shock is expected to form for non-rotating GOHS (Chakrabarti & Sahu 1997), is affected for rotating ones and the possibility of shock is reduced for both the co-rotating and counter-rotating GOHS. Also for the higher viscosity, both the inner and outer sonic point branch merge each other and the possibility of shock reduces again. Actually, with the increase of viscosity, the region containing the X-type sonic points tends to become O-type. The outer X-type sonic points recede in more outwards and inner ones proceed in inwards more.

As the incoming matter slows down abruptly at the surface of GOHS as well as at two shock locations (if two shocks form), the overall density becomes higher for an accretion disk around GOHS compared to that for a black hole. Similarly, the temperature in the disk is higher for a GOHS. Thus, one can conclude that the accretion disk around a GOHS is very favourable for nucleosynthesis. We know that the accretion disk around a black hole is enough hot for nucleosynthesis (Mukhopadhyay 1999; Mukhopadhyay & Chakrabarti 2000), which are different from that in star. As the density and temperature may be much higher around a GOHS, more efficient nucleosynthesis is expected particularly at an inner region of the disk. In early, we saw that the high temperature of accretion disk around a black hole is very favourable for the photo-dissociations and the proton capture reactions (Chakrabarti & Mukhopadhyay 1999). As the accretion disk around a GOHS is hotter than that around a black hole, even $`{}_{}{}^{4}He`$ which has a high binding energy, may dissociate into deuterium and then into proton and neutron. If we consider the accreting matter to come from the nearby ’Sun like’ companion star, the initial abundance of $`{}_{}{}^{4}He`$ in the accreting matter is supposed to be about $`25\%`$. Therefore, by the dissociation of this $`{}_{}{}^{4}He`$, neutron may produce in a large scale, which could give rise to the neutron rich elements. Guessoum & Kazanas (1999) showed that the profuse neutron may be produced in the accretion disk and through the spallation reactions lithium may be produced in the atmosphere of the star. When the neutron comes out from the accretion disk by the formation of an outflow, in that comparatively cold environment, $`{}_{}{}^{7}Li`$ may be produced which can be detected on the stellar surface. In early, it was shown that the metalicity of the galaxy may be influenced when outflows form in the hot accretion disk around black holes (Mukhopadhyay & Chakrabarti 2000). In case of the lighter galaxy, the average abundance of the isotopes of $`Ca`$, $`Cr`$ may significantly change. Also the abundance of lighter elements, like the isotopes of $`C`$, $`O`$, $`Ne`$, $`Si`$ etc. may increase significantly. As the temperature of the accretion disk around GOHS is higher, the expected change of abundance of these elements and the corresponding influence on the metalicity of the galaxy is expected to be high.

We know that, out of an observed pair of kilohertz QPO frequencies for a particular candidate, one is the oscillation due to its Keplerian motion (Osherovich & Titarchuk 1999). For 4U 1636-53, lower frequency is 950Hz (Wijnands et al. 1997) which may be the Keplerian one. According to our pseudo-potential, if we calculate this Keplerian frequency for 4U 1636-53, it comes out to be 1027Hz for $`M=1.4M_{}`$ and $`x_K=18`$, which is not much over estimated with respect to the observed one. Thus, we can give a theoretical prediction of QPO, at least for one, out of a pair. Similarly, for other candidates, one can calculate Keplerian frequency according to our potential and compare with observation. The detailed theoretical study of QPO based on this present scheme will be pursued elsewhere.

## Acknowledgments

SG is thankful to Inter-University Centre for Astronomy and Astrophysics for Summer Training Programme, 2002, where this work was initiated.


format: text

text:
# Low-energy singlet dynamics of spin-1/2 kagomé Heisenberg antiferromagnets

## I Introduction

Unusual low-temperature properties of kagomé antiferromagnets (KAFs) attracted much attention of both theorists and experimenters in last decade. Seemingly the most intriguing features were observed in specific heat $`C`$ measurements of SrCr9pGa12-9pO19 (spin-$`\frac{3}{2}`$ kagomé material) Ramirez et al. (2000). It has been revealed a peak at $`T5`$ K which was practically independent of magnetic field up to 12 T and $`C`$ appeared to be proportional to $`T^2`$ at $`T5`$ K.

There was no appropriate theory describing KAF low-energy sector so far. Qualitative understanding of the KAF low-temperature physics was based mostly on results of numerous finite cluster investigations Lecheminant et al. (1997); Zeng and Elser (1990); Leung and Elser (1993); Waldtmann et al. (1998); Zeng and Elser (1995). They revealed a gap separating the ground state from the upper triplet levels and a band of nonmagnetic singlet excitations with a very small or zero gap inside the spin gap. The number of states in the band increases with the number of sites $`N`$ as $`\alpha ^N`$. It was obtained for samples with up to 36 sites that $`\alpha =1.15`$ and 1.18 for even and odd $`N`$, respectively Lecheminant et al. (1997); Waldtmann et al. (1998). It is believed now that this wealth of singlets is responsible for specific heat low-$`T`$ peak and explains its field independence Sindzingre et al. (2000); Ramirez et al. (2000).

The origin of the singlet band as well as the nature of the ground state were unclear until now. Previous exact diagonalization studies of clusters with $`N36`$ Chalker and Eastmond (1992); Leung and Elser (1993) revealed exponential decay of the spin-spin and dimer-dimer correlation functions. So the point of view that KAF is a spin liquid is widely accepted Sindzingre et al. (2000); Lecheminant et al. (1997); Sachdev (1992); Zeng and Elser (1995); Yang et al. (1993); Kalmeyer and Laughlin (1987); Marston and Zeng (1991); Waldtmann et al. (1998); Leung and Elser (1993); Chalker and Eastmond (1992).

It seems the best candidate for description of KAF low-energy properties is a quantum dimer model Sachdev (1992); Zeng and Elser (1995); Mila (1998). In Ref. Mambrini and Mila (2000) an approach firstly pioneered by Subrahmanyam in Ref. Subrahmanyam (1995) was developed in which a spin-$`\frac{1}{2}`$ kagomé lattice is considered as a set of interactive triangles with a spin in each apex. It was suggested there to work in the subspace where the total spin of each triangle is 1/2 investigating low-lying excitations. In spite of certain success of this approach (coincidence of the low-energy spectrum and the number of lower singlet excitations in samples with up to 36 sites with exact diagonalization results) a further development of this approach is required to get the full physical description of KAF.

Recently in Ref. Kotov et al. (2001) it was suggested a model of frustrated antiferromagnet which low-energy properties can be generic for KAF as well as for two other types of frustrated magnets which possess a similar behavior as KAF and have many singlet states inside the triplet gap: pyrochlore Canals and Lacroix (1998) and $`\mathrm{CaV}_4\mathrm{O}_9`$ Albrecht et al. (1996). Weakly interactive plaquettes on the square lattice were considered there. Each plaquette has two almost degenerate singlet ground states, so a band of singlet excitations arises if the inter-plaquette interaction is taken into account. It is shown that there is a quantum phase transition in the model at a critical value of frustration separating a disorder plaquette phase and a columnar dimer one. In the proximity of this transition the specific heat has a low-$`T`$ peak below which it possesses a power low-temperature dependence.

In our recent paper Ref. Syromyatnikov and Maleyev (2002) we have shown that a similar picture is relevant for spin-$`\frac{1}{2}`$ KAF. It was proposed to consider a kagomé lattice as a set of stars with 12 spins arranged in a triangular lattice (see Fig. 1). Numerical diagonalization has shown that star has two degenerate singlet ground states separated from the upper triplet levels by a gap. These states form a singlet energy band as a result of inter-star interaction. It was assumed that this band determines the KAF low-energy singlet sector. We have shown that it is described by Hamiltonian of a magnet in the external magnetic field where degenerate states of the stars are represented in terms of two projections of pseudospin-$`\frac{1}{2}`$.

The main advantage of a star over a triangle proposed in Refs. Mambrini and Mila (2000); Subrahmanyam (1995) as a starting point to study KAF is that consideration of star two singlet ground states explicitly leads to effective Hamiltonian which describes singlet low-temperature dynamics.

In the present more comprehensive paper we develop this approach. It is proved using symmetry considerations presented in Sec. II that the singlet band above discussed determines the KAF lower singlet sector.

This band is studied in Sec. III. The general form of the effective Hamiltonian is established there. The Hamiltonian parameters are calculated up to the third order of perturbation theory. The ground state energy calculated in the considered model is lower than that evaluated numerically in the previous finite clusters studies. We demonstrate that our approach can not be expanded directly on KAFs with larger spin values. A way of experimental verification of this picture using neutron scattering is discussed.

Finally we summarize our results in Sec. IV.

## II Symmetry consideration

We start with the Hamiltonian of the spin-$`\frac{1}{2}`$ kagomé Heisenberg antiferromagnet:

$$_0=J_1\underset{i,j}{}𝐒_i𝐒_j+J_2\underset{(i,j)}{}𝐒_i𝐒_j,$$
(1)
where $`i,j`$ and $`(i,j)`$ denote nearest and next-nearest neighbors on the kagomé lattice, respectively, shown in Fig. 1. The case of $`|J_2|J_1`$ is considered in this paper. We discuss a possibility of both signs of next-nearest-neighbor interactions — a ferromagnetic and an antiferromagnetic one. As is shown below, in spite of the smallness the second term in Eq. (1) can be of importance for the low-energy properties.

Kagomé lattice can be represented as a set of stars arranged in a triangular lattice (see Fig. 1). To begin with we neglect interaction between stars and put $`J_2=0`$ in Eq. (1). A star is a system of 12 spins. Let us consider its properties in detail.

As the Hamiltonian Eq. (1) commutates with all the projections of the total spin operator, all the star levels are classified by the values $`S`$, irreducible representations (IRs) of its symmetry group and are degenerated with $`S_z`$. The star symmetry group $`C_{6v}`$ contains six rotations and reflections with respect to six lines passing through the center. There are four one-dimensional and two two-dimensional IRs which are presented in Appendix A. In their basis the matrix of the Hamiltonian has a block structure. Each block has been diagonalized numerically.

As a result it was obtained that the star has doubly degenerated singlet ground state separated from the lower triplet level by a gap $`\mathrm{\Delta }0.26J_1`$. Ground state wave functions can be represented as follows:

$`\mathrm{\Psi }_+`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2+1/16}}}(\varphi _1+\varphi _2),`$ (2)
$`\mathrm{\Psi }_{}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{21/16}}}(\varphi _1\varphi _2),`$ (3)
where functions $`\varphi _1`$ and $`\varphi _2`$ are shown schematically in Fig. 2. A bold line there represents the singlet state of the corresponding two spins, i.e. $`(_i_j_i_j)/\sqrt{2}`$.

It can be shown that $`\varphi _1`$ and $`\varphi _2`$ are not orthogonal: their scalar product is $`(\varphi _1\varphi _2)=1/32`$. They contain six singlets each having the energy $`S(S+1)J_1=3/4J_1`$. One can obtain that interaction between singlets does not contribute to the energy of the ground states which is consequently equal to $`4.5J_1`$.

Functions $`\varphi _1`$ and $`\varphi _2`$ are invariant with respect to rotations of the star and they transform into each other under reflections. Hence $`\mathrm{\Psi }_+`$ is invariant under all the symmetry group transformations. In contrast, function $`\mathrm{\Psi }_{}`$ is invariant under rotations, changes the sign with reflections and is transformed under the representation (13). So the ground state has accidental degeneracy. As it is shown in the next section, the next-nearest neighbor interaction, which has the same symmetry as the original Hamiltonian, lifts this degeneracy.

KAF containing $`𝒩`$ non-interactive stars has an energy spectrum with great levels degeneracy when $`𝒩1`$. For example, the ground state degeneracy is $`2^𝒩`$ and that of the lowest triplet level is $`3𝒩2^{(𝒩1)}`$. Interaction between stars gives rise to an energy band from every such group of levels and it is a very difficult task to follow their evolution. Meanwhile group theory allows one to make some conclusions about the KAF low-energy sector. We show now that the singlet band stemmed from the ground state can not be overlapped by those originated from the upper singlet levels.

Let us consider a cluster with seven stars shown in the Fig. 3(a) neglecting for a beginning interaction between them. The symmetry group of the cluster is $`C_{6v}`$ also. The ground state has degeneracy $`2^7=128`$. Corresponding wave functions transformed under IRs of $`C_{6v}`$ are constructed as linear combinations of products containing $`\mathrm{\Psi }_+`$ or $`\mathrm{\Psi }_{}`$ for each star. It is easy to show using a standard procedure of irreducible representations basis construction Landau and Livshitz ; Petrashen and Trifonou  and Appendix A that there is at least two ground state wave functions of the cluster discussed transformed under any given IR.

It is important to mention that operator of inter-star interaction in the cluster has the same symmetry as the intra-star one which is a sum of the stars Hamiltonians. Then the interaction between stars commutates with the square of the total spin operator. So if we begin to increase the inter-star interaction from zero, all the levels would move by energy but their classification can not be changed. Levels can cross with each other as the interaction arises from 0 to $`J_1`$ but the crossing is forbidden for those of the same symmetry. It is a consequence of the symmetry theorem proved in Refs. Petrashen and Trifonou ; Landau and Livshitz . Hence, one can conclude that the lower singlet sector of the cluster is formed by states which stem from the origin 128 lower levels.

We can lead to the same conclusions considering clusters of the symmetry $`C_{6v}`$ with the greater number of stars. So we will assume in the following that the KAF low-energy singlet sector is formed by the states originated from those in which each star is in one of the states $`\mathrm{\Psi }_+`$ or $`\mathrm{\Psi }_{}`$. Because bands with $`S0`$ can overlap the singlet one under discussion we have to suppose that the KAF low-energy properties are determined by the lowest states in this singlet band.

As interaction between stars commutates with the square of the total spin operator, bands of the different $`S`$ can be studied independently. The KAF lower singlet sector is considered in detail in the next section. Investigation of states with $`S0`$ is out of the scope of this paper.

## III Singlet dynamics

In this section we derive the general form of the effective Hamiltonian described the lower singlet sector. Inter-star interaction is considered as a perturbation. Despite of its not being small comparing to the intra-star one there are reasons presented below to use perturbation expansion here.

Two-stars coupling.—We start with consideration of the interaction between two nearest stars still neglecting the second term in Eq. (1). Initially there is four times degenerated ground state with wave functions $`\{\mathrm{\Psi }_{n_1}^{(1)}\mathrm{\Psi }_{n_2}^{(2)}\}`$ (here $`n_i=+,`$ and the upper index labels the stars) and the energy $`E_{n_1n_2}^{(0)}=E_{n_1}^{(0)}+E_{n_2}^{(0)}=9J_1`$. As it is seen from Fig. 4, the interaction has the form

$$V=J_1(𝐒_1^{(1)}𝐒_1^{(2)}+𝐒_3^{(1)}𝐒_3^{(2)}).$$
(4)
According to the standard theory Landau and Livshitz  the following conditions should be fulfilled to consider $`V`$ as a perturbation:

$$\left|C_{m_1m_2}^{n_1n_2}\right|=\left|\frac{V_{n_1n_2;m_1m_2}}{E_{n_1n_2}^{(0)}E_{m_1m_2}^{(0)}}\right|1,$$
(5)
where $`V_{n_1n_2;m_1m_2}=\mathrm{\Psi }_{n_1}^{(1)}\mathrm{\Psi }_{n_2}^{(2)}V\mathrm{\Psi }_{m_1}^{(1)}\mathrm{\Psi }_{m_2}^{(2)}`$, $`m_1m_2`$ denotes excited singlet levels of the two stars and $`n_i=+,`$. We have calculated $`C_{m_1m_2}^{n_1n_2}`$ for $`n_i=+,`$ using wave functions obtained numerically and found that all of these coefficients do not exceed $`0.09`$. So conditions Eq. (5) are fulfilled. Then the maximum value of the sum $`_{m_1m_2}|C_{m_1m_2}^{n_1n_2}|^2`$ is $`0.28`$ which is also small enough. Thus interaction between stars will be considered as a perturbation in the following.

We proceed with calculations of corrections to the initial ground state energy of two stars. For this purpose, as the state is fourfold degenerated, one has to solve a secular equation Landau and Livshitz . The corresponding matrix elements in the third order of perturbation theory are given by Landau and Livshitz

$$\begin{array}{c}H_{n_1n_2;k_1k_2}=V_{n_1n_2;k_1k_2}+\underset{m_1,m_2}{}\frac{V_{n_1n_2;m_1m_2}V_{m_1m_2;k_1k_2}}{E_{n_1n_2}^{(0)}E_{m_1m_2}^{(0)}}\hfill \\ \hfill +\underset{m_1,m_2}{}\underset{q_1,q_2}{}\frac{V_{n_1n_2;m_1m_2}V_{m_1m_2;q_1q_2}V_{q_1q_2;k_1k_2}}{(E_{n_1n_2}^{(0)}E_{m_1m_2}^{(0)})(E_{n_1n_2}^{(0)}E_{q_1q_2}^{(0)})},\end{array}$$
(6)
where $`n_i,k_i=+,`$. Obviously the first term in Eq. (6) is zero and the second one can be represented as follows:

$$\begin{array}{c}H_{n_1n_2;k_1k_2}=i_0^{\mathrm{}}𝑑te^{\delta t+iE_{n_1n_2}^{(0)}t}\hfill \\ \hfill \times \mathrm{\Psi }_{n_1}^{(1)}\mathrm{\Psi }_{n_2}^{(2)}Ve^{it(_0^{(1)}+_0^{(2)})}V\mathrm{\Psi }_{k_1}^{(1)}\mathrm{\Psi }_{k_2}^{(2)},\end{array}$$
(7)
where $`_0^{(i)}`$ are Hamiltonians of the corresponding stars. The third term in Eq. (6) will be considered later. Using the symmetry of the functions $`\mathrm{\Psi }_+`$ and $`\mathrm{\Psi }_{}`$ discussed above and invariance of the system under reflection with respect to the dotted line plotted in Fig. 4, one can show that only nonzero elements belong to the first and to the second diagonals (i.e. with $`n_1=k_1`$, $`n_2=k_2`$ and with $`n_1k_1`$, $`n_2k_2`$). We have calculated them numerically with a very high precision by expansion of the operator exponent up to the power 150. The results can be represented in the following form:

$`H_{++;++}`$ $`=`$ $`a_1+a_2a_3,`$ (8a)
$`H_{+;+}`$ $`=`$ $`a_1+a_3,`$ (8b)
$`H_{+;+}`$ $`=`$ $`a_1+a_3,`$ (8c)
$`H_;`$ $`=`$ $`a_1a_2a_3,`$ (8d)
where $`a_1=0.256J_1`$, $`a_2=0.015J_1`$ and $`a_3=0.0017J_1`$. Terms of the second diagonal $`H_{++;}=H_{;++}=H_{+;+}=H_{+;+}=a_4=0.0002J_1`$ are much lesser then $`a_1`$, $`a_2`$ and $`a_3`$. So the interaction shifts all the levels on the value $`a_1`$ and lifts their degeneracy. Constants $`a_2`$, $`a_3`$ and $`a_4`$ determine the levels splitting. It is seen that the splitting is very small comparing to the shift.

So KAF appears to be a set of two-levels interacting systems and one can naturally represent the low-energy singlet sector of Hilbert space in terms of pseudospins: $`=\mathrm{\Psi }_{}`$ and $`=\mathrm{\Psi }_+`$. It is seen from Eqs. (8) that in these terms the interaction between stars is described by the Hamiltonian of a ferromagnet in the external magnetic field:

$$=\underset{i,j}{}[𝒥_zs_i^zs_j^z+𝒥_ys_i^ys_j^y]+h\underset{i}{}s_i^z+𝒞$$
(9)
where $`i,j`$ labels now nearest-neighbor pseudospins, arranged in a triangular lattice formed by the stars, $`𝐬`$ is spin-$`\frac{1}{2}`$ operator, $`𝒞=5.268J_1𝒩`$, $`𝒥_z=4a_3=0.007J_1`$, $`𝒥_y=4a_4=0.001J_1`$ and $`h=6a_2=0.092J_1`$. Here $`𝒩=N/12`$ is the number of stars in the lattice. Factor 6 appears in the expression for $`h`$ because each star interacts with 6 neighbors. One can see that the magnetic field in the effective Hamiltonian Eq. (9) is much more larger than the exchange. So in this approximation stars behave like almost free spins in the external magnetic field and the ground state of KAF has long-range singlet order which settles on the triangular star lattice and is formed by stars in $`\mathrm{\Psi }_{}`$ states.

$`V^3`$ corrections.—The field remains the largest term of the effective Hamiltonian and the KAF ground state has the same long-range order if one takes into account $`V^3`$ terms in the perturbation series. For two-stars coupling $`V^3`$ corrections have the form given by Eq. (6). Terms $`V^3`$ imply analysis of the three-stars interaction as well. Nonzero contributions from them one gets only for the configuration presented in Fig. 3(b). The secular matrix for three stars is of the size $`8\times 8`$. We have calculated $`V^3`$ corrections with a very high precision using the integral representation similar to that Eq. (7) for the second term in Eq. (6). All the operator exponents were expanded up to the power 150. As a result the low-energy properties of the KAF are described by the effective Hamiltonian:

$$=\underset{i,j}{}\left[𝒥_zs_i^zs_j^z+𝒥_xs_i^xs_j^x+𝒥_ys_i^ys_j^y\right]+h\underset{i}{}s_i^z+𝒞,$$
(10)
where all parameters are given in the Table 1. It describes two-stars coupling. We omit the three pseudospins terms in Eq. (10) which have the form $`s_i^zs_j^zs_k^z`$ and $`s_i^zs_j^ys_k^y`$ and describe three-stars interaction. The corresponding coefficients are of the order of $`10^3J_1`$ and $`10^4J_1`$, respectively, and are negligible in comparison with those of the presented terms. It should be stressed that within our precision Hamiltonian Eq. (10) is an exact mapping of the original Heisenberg model in the low-energy sector (excitation energy $`\omega \mathrm{max}\{𝒥_z,𝒥_y,𝒥_x,h\}J_1`$).

As one can conclude from the Table 1 studying $`V^3`$ corrections the common shift given by them remains much larger than the levels splitting in both cases of two- and three-stars coupling. Meanwhile the values of $`V^3`$ perturbation terms are approximately two times smaller than those of $`V^2`$ ones. So the change in the effective Hamiltonian from $`V^3`$ terms is significant and analysis of the perturbation series can not be finished at this point for correct establishing of the effective Hamiltonian. Unfortunately such a work requires great computer capacity we do not dispose now. We have to restrict ourself with this precision here.

One can judge about the applicability of the perturbation series from the following values of the ground state energy of two interacting stars shown in Fig. 4 obtained numerically and using first two orders of perturbation theory. Ground state energy of two non-interacting stars is $`9J_1`$. That of two interacting ones which we have calculated numerically by power method Parlett  is $`9.62J_1`$. Whereas the ground state energy obtained using Table 1 is $`9.42J_1`$ (contribution from $`V^2`$ and $`V^3`$ terms are $`0.27J_1`$ and $`0.15J_1`$, respectively).

Effective Hamiltonian structure.—Despite the perturbation theory works quite badly in the star model and many perturbation terms are to be taken into account, we can proof now that Eq. (10) is the most general form of the effective Hamiltonian assuming that $`n`$-pseudospins couplings with $`n>2`$ are small as it was in the case of $`n=3`$ discussed above. Let us consider possible terms of the form $`s_i^zs_j^+`$, $`s_i^zs_j^{}`$, $`s_i^+`$ and $`s_i^{}`$. In these cases numbers of functions $`\mathrm{\Psi }_+`$ and $`\mathrm{\Psi }_{}`$ to the right hand of the corresponding matrix elements differ from those to the left hand by unit. As it was also pointed out above, kagomé lattice contains lines of symmetry reflections and star Hamiltonian and inter-star interaction are invariant with respect to these reflections. As $`\mathrm{\Psi }_+`$ are invariant and $`\mathrm{\Psi }_{}`$ changes the sign under these transformations, the matrix elements are equal to themselves with the opposite sign and so must be zero. Other possible term $`s_i^xs_j^y=(i/4)(s_i^+s_j^+s_i^{}s_j^{}+s_i^{}s_j^+s_i^+s_j^{})`$ can not appear in the effective Hamiltonian because the corresponding matrix elements should be imaginary.

Ground state.—As is clear from the Table 1, $`𝒥_x`$ and $`h`$ are the largest parameters of the Hamiltonian Eq. (10) in our approximation. So KAF behaves like Ising antiferromagnet in perpendicular magnetic field. In this case the classical value of the field at which spin-flip occurs is $`h_{sf}=𝒥_x`$ which is approximately 2.6 times smaller than $`h`$. So the ground state should remain ordered with all the stars in $`\mathrm{\Psi }_{}`$ state.

The ground state energy and that of the upper edge of the singlet band calculated using Table 1 are $`(4.5J_1+\mathrm{\Delta }𝒞+h/2+3𝒥_z/4)𝒩=0.452J_1N`$ and $`(4.5J_1+\mathrm{\Delta }𝒞h/2+3𝒥_z/4)𝒩=0.437J_1N`$, respectively. Corrections from $`𝒥_x`$ to these values in the first non-zero order of perturbation theory are given by $`(3/16)𝒩𝒥_x^2/h`$ and are negligible. At the same time the ground state energy of the largest cluster with $`N=36`$ considered before numerically was $`0.438J_1N`$ Waldtmann et al. (1998). So we believe that clusters used in the previous studies were too small to reflect the Heisenberg KAF low-energy sector at $`J_2=0`$ properly.

Interaction $`J_2`$.—We show now that in spite of its smallness the next-nearest neighbor interaction can play an important role for low-energy properties. We have calculated $`J_2`$ corrections to the parameters of the effective Hamiltonian Eq. (10) for the first and for the second terms in Eq. (6) only. There are 12 intrinsic $`J_2`$ interactions in each star which splits the doubly degenerated ground state and, as is seen from the Table 1, gives a contribution to the magnetic field $`h`$ and to the constant $`𝒞`$.

As is clear from Fig. 4 the two-stars coupling is given now by the operator $`\stackrel{~}{V}=J_2(𝐒_1^{(1)}𝐒_2^{(2)}+𝐒_2^{(1)}𝐒_1^{(2)}+𝐒_2^{(1)}𝐒_3^{(2)}+𝐒_3^{(1)}𝐒_2^{(2)})`$. Corrections proportional to $`J_2`$ were calculated by the same way as above and are presented in the Table 1 too. It is seen that the contribution of next-nearest interactions to the magnetic field becomes significant if $`|J_2|0.1J_1`$. If $`J_2<0`$ (ferromagnet interaction) they can even change the sign of $`h`$.

Effect of next-nearest ferromagnetic coupling for KAF properties was previously studied in Ref. Lecheminant et al. (1997) numerically on finite clusters with $`N27`$ in the wide range of the values of $`J_2`$. It was shown there that at $`|J_2|/J_11`$ the ground state has $`\sqrt{3}\times \sqrt{3}`$ magnetic structure. At $`|J_2|/J_11`$ the ground state was found to be disorder and there is a band of singlet excitations inside the triplet gap. As we demonstrated above in our approach this band is a result of star ground state degeneracy.

$`T^2`$ specific heat behavior.—There were many speculations on the low-temperature dependence of the KAF specific heat $`CT^2`$ observed experimentally for $`S=3/2`$ (see Refs. Sindzingre et al. (2000); Kotov et al. (2001) and references therein). As we obtained above, low-$`T`$ properties are described by effective Hamiltonian Eq. (10) of a magnet which has spectrum of the form $`ϵ_𝐪=\sqrt{(cq)^2+\mathrm{\Delta }^2}`$ at $`q1`$ and can be in ordered or disordered phases depending on particular values of the Hamiltonian parameters. Small $`\mathrm{\Delta }^{}`$ here signifies the proximity to the quantum critical point in which $`CT^2`$. Such a situation arises in the singlet dynamics of the model of interactive plaquets Kotov et al. (2001) discussed above. We do not present here the corresponding analysis because parameters of the effective Hamiltonian could be changed in the further orders of the perturbation theory.

Experimental verification.—In both cases of ordered ground state and disordered one the approach presented in this paper can be checked by inelastic neutron scattering: corresponding intensity for the singlet-triplet transitions should have periodicity in the reciprocal space corresponding to the star lattice. This picture is similar to observed one in the dimerized spin-Pairls compound $`\mathrm{CuGeO}_3`$ Regnault et al. (1996). In this case inelastic magnetic scattering has a periodicity which corresponds to the dimerized lattice.

Cases of $`S>1/2`$.—It should be noted that our consideration of $`S=1/2`$ KAF can not be expanded directly on the cases of larger spins. Although functions presented in Fig. 2, where a bold line denotes the singlet state of the corresponding two spins, remain eigenfunctions of the Hamiltonian for $`S>1/2`$, we have found numerically that they are not ground states of the star with $`S=1`$ and $`S=3/2`$. All details of calculations are presented in Appendix B. So an other approach should be proposed for KAFs with $`S>1/2`$.

## IV Conclusion

In this paper we present a model for the low-energy physics of spin-$`\frac{1}{2}`$ kagomé Heisenberg antiferromagnets (KAFs). The spin lattice can be represented as a set of stars which are arranged in a triangular lattice and contain 12 spins (see Fig. 1). Each star has two degenerate singlet ground states with different symmetry which can be described in terms of pseudospin. It is shown that interaction between the stars leads to the band of singlet excitations which determines the KAF low-energy properties. The low-energy dynamics is described Hamiltonian of a spin-$`\frac{1}{2}`$ magnet in the external magnetic field given by Eq. (10). The Hamiltonian parameters are calculated in the first three orders of perturbation theory and are summarized in the Table 1. Within our precision KAF has an ordered singlet ground state with all the stars in the state given by Eq. (3). The ground state energy is lower than that calculated in previous finite clusters studies. We show that our model can not be expanded directly on KAFs with $`S>1/2`$.

The approach discussed in this paper can be verified experimentally on inelastic neutron scattering: the corresponding intensities for singlet-triplet transitions should have periodicity in the reciprocal space corresponding to the star lattice.

We are grateful to A.G. Yashenkin for interesting discussions. This work was supported by Russian State Program ”Collective and Quantum Effects in Condensed Matter”, the Russian Foundation for Basic Research (Grants No. 03-02-17340, 00-15-96814) and Russian State Program ”Quantum Macrophysics”.

## Appendix A Irreducible representations of the group $`C_{6v}`$

The symmetry group $`C_{6v}`$ contains six rotations $`C^k`$ on the angles $`2\pi k/6`$ ($`k=0,1,..,5`$) and six reflections which can be written as $`C^ku_1`$, where $`u_1`$ is operator of a reflection. One-dimensional irreducible representations can be presented as follows Petrashen and Trifonou ; Landau and Livshitz :

$`C^k`$ $`1,`$ $`u_1`$ $`1;`$ (11)
$`C^k`$ $`(1)^k,`$ $`u_1`$ $`1;`$ (12)
$`C^k`$ $`1,`$ $`u_1`$ $`1;`$ (13)
$`C^k`$ $`(1)^k,`$ $`u_1`$ $`1.`$ (14)
For two-dimensional ones we have Petrashen and Trifonou ; Landau and Livshitz

$`C^k`$ $`\left(\begin{array}{cc}e^{i\frac{2\pi l}{6}k}& 0\\ 0& e^{i\frac{2\pi l}{6}k}\end{array}\right),`$ $`u_1`$ $`\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),`$ (19)
where two inequivalent representations are given by $`l=1`$ and $`l=2`$.

## Appendix B Star with $`S=1`$ and $`S=3/2`$

In this Appendix we present details of numerical calculations shown that functions presented in Fig. 2, where a bold line denotes the singlet state of the corresponding two spins, are not ground states of the star with $`S=1`$ and $`S=3/2`$ as it was for $`S=1/2`$.

A simple numerical method for determination of the maximum by module eigenvalue of a Hermitian operator $`H`$ (power method Parlett ) was used. It is based on the following statement. Let us consider a state of the system $`f=_ic_i\psi _i`$, where the sum may not include all the $`H`$ eigenfunctions. The maximum by module eigenvalue $`E_{extr}`$ for given $`f`$ is determined by the expression

$$\underset{n\mathrm{}}{lim}\frac{f|H^{n+1}|f}{f|H^n|f}=E_{extr}.$$
(20)
It becomes evident if one notes that $`f|H^n|f=_i|c_i|^2E_i^n`$.

The Eq. (20) can be used in numerical calculations in the following way. The corresponding expression is calculated for $`n=1,2,..,n_{max}`$. So one can control the convergence comparing results with different $`n`$. Studying a full set of vectors $`f`$ and taking large enough $`n_{max}`$ to match the necessary precision one can find the largest by module eigenvalue of $`H`$.

In the case of the star the maximum eigenvalue of the Hamiltonian is $`E_{max}=18S^2J_1`$ (this energy has the state in which all the spins are along the same direction) and the energy of singlet states shown in Fig. 2 is $`E_{ss}=6S(S+1)J_1`$. As $`E_{max}>|E_{ss}|`$ for $`S>1/2`$, we have to take $`H=_0WI`$ to investigate the lower $`_0`$ levels, where $`_0`$ is the star Hamiltonian given by Eq. (1), $`I`$ is an identical matrix and $`W=(E_{max}+E_{ss})/2+J_1`$. Thus eigenvalues of $`H`$ are shifted down according to those of $`_0`$ by the same value $`W`$ so as the largest by module $`H`$ eigenvalue to be equal to the $`_0`$ ground state energy minus $`W`$.

We have not found the ground state energy for the star with $`S=1`$ and $`S=3/2`$ by this method because a full set of vectors $`f`$ should be examined for that. This operation demands much computer time. But we have obtained studying a number of vectors $`f`$ that there are states lower by the energy at least by the value $`1.8J_1`$ than those discussed above. The method gave $`E_{extr}`$ with the prescribed precision to the second decimal place at $`n_{max}=100÷300`$ depending on $`f`$ and $`S`$.