text:
# Green’s Function Approach to B. E. Condensation in the alkali atoms

## Abstract

A new scheme has been proposed to solve the B. E. condenstates in terms of Green’s function approach. It has been shown that the radial wave function of two interacting atoms, moving in a common harmonic oscillator potential modified by an effctive interaction, satisfies an intregral equation whose kernel is separable. The solution of the integral equation can be written in terms of the harmonic oscillator wave functions. The ground state wave function of the system can be written in terms of these solutions.

Department of Physics, Presidency College, Calcutta 700 073,India.

PACS numbers: 03.75Fi,02.60Cb,32.80Pj Keywords : B. E. Condensate, Effective Interaction, Green’Function.

The recent experimental realization of Bose-Einstein condensation \[BEC\] in ultra cold atomic gases has generated intensive theoretical and practical interest regarding the properties of Bose gas. Non-linear Schrodinger equation in this context known as Gross-Pitaevskii \[ G-P \] equation ,has been used to study the properties of the interacting spinless Bosons , trapped in harmonic potential . Altough the G-P equation is widely accepeted as a workable model for description of the B E condenstates at low temperature, the actual dynamics of B E condenstates is not fully captured by the G-P equation because of the non-integrability of the G-P equation and explicit solution is not known . It is to be noted that for the positive scattering length, the G-P equation is generally solved numerically only under cylindrically symmetric systems. The Thomas Fermi approximation  has been introduced to get some analytical solution for the description of B E condenstates. On the other hand for the negative scattering length the G-P equation is yet to be properly investigated properly. Another alternative approach to study the dynamics of the B.E.condenstates is the time dependent variation technique  in which one assumes fixed profile of the condenstate and computes the evolution of some parameters such as the width by variational techniques. It should also be mentioned that the negative scattering length provides stable solution of the G-P equation only under certain conditions for the number of particles and the size of the trap . When those conditions are not fulfiled the condenstates become unstable. Therefore realization of large condenstates will provide reliable experimental data for the precise understanding of the condensation process. In this paper we shall present a new solution for the condenstate on the basis of Green’s function approach of scattering theory. Similar method was used in the treatment of the nuclear many body problems by Bethe and Goldstone; Brueckner and Gammel  for the understanding of the various aspects nuclear many body problem. We have derived an integral equation for the description of two interacting bosons in the presence of many bosons. We shall consider a system of N spinless interacting atoms all having the same hyperfine species, in equilibrium and at very low temperature and densities characteristic of alkali atoms under Bose- Einstein condensation (BEC) condition. We assume that each of the atoms is moving in the influence of a common harmonic oscillator potential which is called a trapping potential($`V_{ext}`$) and an effective two-body interaction acting between them which modifies their behaviour . The essential questions we wish to investigate are : How does the effective two body interaction modify the two-body wave function?. How does the energies of the two interacting bosons system change when the two- body effective interaction is switched on ? It is welknown that the scattering length ($`a_s`$) is very small compared to the interatomic distant($`r_{int}`$)  under BEC in the dilute atomic alkali gases. Therefore the behaviour of the atoms is extremely sensitive to the details of the trap and the effective two- body interaction. Since the effective two- body interaction is very weak the quantitative calculation based on perturbation theory in the interatomic interaction is hoped to be highly reliable for the alkali gases. Furthermore the de Broglie wave length associated with individual atom is comparable to interatomic $`r_{int}`$ distance therefore the correlation between atoms becomes important. Hence the construction of the many -body condenstate state wave function of the system in terms of two-body wave function will be interesting. The standard many-body wave function in terms of Hartree-Fock ansatz using the single particle solution of the G-P equation does not take into accoutn the effect of correlation between atoms which is very important in the case of BEC alkali gases.In this letter we have presented an analytical formalism to investigate the above aspects of BEC problem without any numerical estimate for a particular choice of the effective interaction.That will be communicated shortly.

Let us consider a system of interacting $`N`$ spinless atoms and now arbitrarily pick out any two interacting spinless atoms, say 1 and 2, denoting their relative and center-of- mass co-ordinates, by $`\stackrel{}{r}`$ and $`\stackrel{}{R}`$ respectively . We shall investigate the state ($`\varphi `$) of these two interacting atoms trapped in a common harmonic potential in the presence $`(N2)`$ interacting atoms. The most general states of the many-body system will then be described by a wave function $`\mathrm{\Psi }(\stackrel{}{r},\stackrel{}{R},\chi ,t)`$ where $`\chi `$ represents schemitically the co-ordinates of the remaining $`(N2)`$ particles. It is the dependence on $`\stackrel{}{r}`$ that is crucial but the parameter dependence of $`\mathrm{\Psi }`$ on $`\stackrel{}{R},\chi `$ and $`t`$ will be ignored . In the present context we assume that simple BEC is realised in the two-particle state $`\varphi (\stackrel{}{r})`$ whose equatiom of motion will be determined using Schrodinger equation. From the knowledge of the wave-function of two-atom system $`\varphi (\stackrel{}{r})`$, the many- body ground state of B E condenstate system of alkali atoms can be written. The hamiltonian for the two- atom system in the common harmonic oscillator potential of frequency $`\omega `$ is given by

$$H_0=\frac{p_1^2}{2m}+\frac{p_2^2}{2m}+\frac{m\omega ^2r_1^2}{2}+\frac{m\omega ^2r_2^2}{2}.$$
(1)
Introducing the cannonical co-ordinate transform

$$\stackrel{}{r}=\stackrel{}{r}_2\stackrel{}{r}_1,\stackrel{}{p}=\frac{(\stackrel{}{p}_2\stackrel{}{p}_1)}{2},$$
(2)
and

$$\stackrel{}{R}=\frac{(\stackrel{}{r}_2+\stackrel{}{r}_1)}{2},\stackrel{}{P}=(\stackrel{}{p}_2+\stackrel{}{p}_1).$$
(3)
The hamiltonian is transformed to

$$H_0=\frac{(P^2+M^2\omega ^2R^2)}{2M}+\frac{(p^2+\mu ^2\omega ^2r^2)}{2\mu },$$
(4)
where M=2m,$`\mu =\frac{m}{2}`$. The Schrodinger equation for the system with $`\mathrm{}=1`$ is

$$[\frac{(_R^2+M^2\omega ^2R^2)}{2M}+\frac{(_r^2+\mu ^2\omega ^2r^2)}{2\mu }]\varphi (\stackrel{}{r})\varphi (\stackrel{}{R})=E^0\varphi (\stackrel{}{r})\varphi (\stackrel{}{R}).$$
(5)
The Schrodinger equation for the relative motion is then

$$H_0\varphi (\stackrel{}{r})=(\frac{^2}{2\mu }+kr^2)\varphi (r)=E^0\varphi (\stackrel{}{r}),$$
(6)
where

$$k=\frac{\mu \omega ^2}{2}=\frac{m\omega ^2}{4}.$$
(7)
The solutions are given by

$$\varphi _{nl}^m(\stackrel{}{r})=\frac{R_{nl}(r)}{r}Y_{lm}(\widehat{r}),$$
(8)
where $`Y_{lm}(\widehat{r})`$ are the usual spherical harmonics. The radial wave function $`R_{nl}(r)`$ is given by

$$R_{nl}(r)=N_{nl}exp(\frac{\nu r^2}{2})r^{l+1}v_{nl}(r),$$
(9)
where $`\nu =\frac{\omega \mu }{\mathrm{}}`$ and $`v_{nl}(r)`$ is the associated Laguerre polynomial:

$$v_{nl}(r)=L_{n+l+\frac{1}{2}}^{l+\frac{1}{2}}(\nu r^2).$$
(10)
$`N_{nl}`$ is determind by the normalization condition and $`E_{nl}^0=\mathrm{}\omega (2n+l+\frac{3}{2})`$.

We shall now treat the effective spin independent two-body interaction $`v(r)`$ which is very weak between two atoms , as a peturbation and solve the Schrodinger equation using the Green’s function method to determine the perturbed wave function for the two- atom system. The unperturbed wave function is given by

$$\varphi _{nl}(\stackrel{}{r})=\frac{R_{nl}(r)}{r}Y_{lm}(\widehat{r}).$$
(11)
Similarly the perturbed wave function can be written as

$$\psi _{n_1l_1}(r)=\frac{u_{n_1l_1}}{r}Y_{l_1m_1}(\widehat{r}).$$
(12)
The Schrodinger equation for the perturbed system is

$$[H_0+v(r)]\frac{u_{n_1l_1}}{r}Y_{m_1l_1}(\widehat{r})=E_{n_1l_1}\frac{u_{n_1l_1}}{r}Y_{m_1l_1}(\widehat{r}),$$
(13)
where $`E_{n_1l_1}`$ is the perturbed energy eigen value. We now write

$$H_0\frac{R_{nl}(r)}{r}Y_{lm}(\widehat{r})=\frac{1}{r}[\frac{1}{m}\frac{^2}{r^2}+V_l(r)]R_{nl}(r)Y_{nl}(\widehat{r})=E_{nl}^0\frac{R_{nl}(r)}{r}Y_{nl}(\widehat{r}).$$
(14)
Multiplying the perturbed equation(13) by $`Y_{n_1l_1}^{}(\widehat{r})`$ on the left and integrating over the solid angle $`d\widehat{r}`$ we obtain

$$[E_{n_1l_1}(\frac{1}{m}\frac{d^2}{dr^2}+V_{l_1}(r)]u_{n_1l_1}(r)=v_{l_1}(r)u_{n_1l_1}(r),$$
(15)
where

$$v_{l_1}(r)=𝑑\stackrel{}{r}Y_{l_1m_1}^{}(\widehat{r})v(r)Y_{l_1m_1}(\widehat{r})$$
(16)
In order to solve this equation we introduce Green‘s $`G_{n_1l_1^{}}(\stackrel{}{r},\stackrel{}{r}^{})`$ corresponding to equ.(15) satisfies the equation

$$(E_{n_1l_1}H_0)G_{n_1l_1}(\stackrel{}{r},\stackrel{}{r}^{})=\delta (\stackrel{}{r}\stackrel{}{r}\mathrm{`}).$$
(17)
Multiplying through on the left by $`Y_{l_1m_1}^{}(\widehat{r})`$ and integrate over the solid angle $`dr`$ we obtain

$$[E_{n_ll_1}(\frac{1}{m}\frac{d^2}{dr^2})+V_{l_1}(r)]g_{n_1l_1}(r,r^{})=\delta (rr^{}),$$
(18)
where

$$g_{n_1l_1}(r,r^{})=\underset{n_2}{}\frac{R_{n_2l_1}(r)R_{n_2l_1}(r^{})}{E_{n_1l_1,n_2l_1}},$$
(19)
where

$$E_{n_1l_1,n_2l_2}=(E_{n_1l_1}E_{n_2l_1}^0).$$
(20)
Hence the perturbed wave function $`u_{n_1l_1}(r)`$ satisfies an integral equation with a separable kernel. Therefore we write

$$u_{n_1l_1}(r)=_0^{\mathrm{}}𝑑r^{}g_{n_1l_1}(r,r^{})v_{l_1}(r^{})u_{n_1l_1}(r^{}).$$
(21)
Finally we write

$$u_{n_1l_1}(r)=\underset{n_2=0}{\overset{\mathrm{}}{}}\frac{R_{n_2l_1}(r)}{E_{n_1l_1},n_2l_1}_0^{\mathrm{}}𝑑r^{}R_{n_2l_1}(r^{})v_{l_1}(r^{})u_{n_1l_1}(r^{}).$$
(22)
The equation for $`R_{nl_1}(r)`$ is

$$[E_{nl_1}^0(\frac{1}{m}\frac{d^2}{dr^2}+V_l{}_{1}{}^{}(r))]R_{nl_1}(r)=o.$$
(23)
Multiplying equ.(15)and equ.(23) on the left by $`R_{nl_1}(r)`$ and $`u_{n_1l_1}`$respectively subtracting and integrating over $`r^{}`$ we obtain

$$E_{n_1l_1}E_{nl_1}^0=\frac{_0^{\mathrm{}}𝑑rR_{nl_1}(r)v_l(r)u_{n_1l_1}(r)}{_0^{\mathrm{}}𝑑rR_{nl_1}(r)u_{n_1l_1}(r)}$$
(24)
and the perturbed wavefunction $`u_{n_1l_1}(r)`$ in a closed form can be written as

$$u_{n_1l_1}(r)=\underset{n_2=0}{\overset{\mathrm{}}{}}K_{n_1l_1,n_2l_1}R_{n_2l_1}(r),$$
(25)
where

$$K_{n_1l_1,n_2l_1}=\frac{1}{E_{n_1l_1,n_2,l_1}}_0^{\mathrm{}}𝑑rR_{n_2l_1}(r)v_{l_1}(r)u_{n_1l_1}(r).$$
(26)
. To this expression of $`u_{n_1l_1}`$ we should add the solution of $`(E_{n_1l_1}H_0)\chi _{n_1l_1}=0`$. Since $`E_{n_1l_1}E_{nl_1}^0`$ the equation $`(E_{n_1,l_1}H_0)\chi _{n_1l_1}=0`$ has no solution.Thus it is not necessary to add a particular solution to the $`u_{n_1l_1}`$ . From equations(25) and (26)it follows that

$$K_{n_1l_1,n_2l_1}=\frac{1}{E_{n_1l_1,n_2l_1}}\underset{n_3=0}{\overset{\mathrm{}}{}}K_{n_1l_1,n_3l_1}_0^{\mathrm{}}𝑑rR_{n_2l_1}(r)v_{l_1}(r)R_{n_3l_1}(r).$$
(27)
.

Results and Discussion

The main conclusions of the paper are the following: It should be noted that $`n_3`$ can in principle run from $`o`$ to $`\mathrm{}`$, the above equ.(27) reprsents an infinite set of simultaneous homogeneous algebraic equations in $`K_{n_1l_1,n_2l_1}`$. The consistency of the these equations imposes the condition of vanishing of the determinant formed from the coefficients of $`K_{n_1l_1,n_2l_1}`$. This condition determines the energy shift $`E_{n_1l_1,nl_1}`$of the two interacting atoms. The corresponding eigenvectors when substituted in eq.(26) determine the perturbed radial function $`u_{n_1l_1}(r)`$. Since B E condensation takes place strictly in the ground state of the system therefore a few values of $`n_3`$ in the sum will provide a very reliable estimate of the two-particle wave function. Because of the separable structure of eq.(21) we can also include various types of central and non central of two-body effective interaction and the variation of the scattering length can also be incorporated to investigate sensitivity of the two-body condenstate wave function and the stability of the Bose-Einstein condensation in alkali atoms on the effective two-body interaction which can not be so easily taken care in the formalism of G-P equation. When the effective two-body interaction becomes non-central in nature then the radial equation for $`u_{n_1l_1}(r)`$ i.e eq.(21) will satisfy a coupled integral equation. When the effective two-body interaction is a pseudo- potential type, as commonly used in the calculation of G-P equation, the determination of the two-body wave function $`u_{n_1l_1}(r)`$ becomes extermely simple in this formalism. We can also investigate the importance two- particle correlation in BEC. Finally the construction of the condenstate wave function describing the BEC of the in terms of the calculated two- particle wave function to study various aspects of BEC becomes possible. These are the main conclusions of the paper.


format: text

text:
# V2051 Oph’s disc evolution on decline from superoutburst

## 1 Introduction

Cataclysmic variables (CVs) are close, interacting binaries in which a Roche-lobe filling main sequence star (the secondary) looses matter to its somewhat more massive white dwarf companion. In the absence of strong magnetic fields the angular momentum conservations leads to the formation of an accretion disc around the primary component through which the matter from the secondary is accreted onto the compact object. A detailed description of such systems is given in Warner (1995), an introduction in Hellier (2001).

If the mass transfer rate is below a certain limit, these discs undergo outbursts, during which the system brightens by a few magnitudes. The cause of such outbursts is the accumulation of material in the disc leading to a temperature increase and subsequently ionization of hydrogen. The changed opacity leads to a sudden further increase in temperature. The system is in outburst. However, with such a low mass accretion rate as is typical for a dwarf nova this state cannot be maintained. Subsequently, due to a decrease in surface density and temperature the hydrogen eventually recombines and the disc returns to its pre-outburst, quiescent state. Osaki (1996) gives an overview on mechanisms and outburst types.

Certain CVs additionally show so-called superoutbursts. These outburst last longer than normal outbursts, but show a similar maximum brightness. Additionally, they are associated with so-called superhumps, short brightenings with a period slightly dissimilar to the orbital period. The superoutbursts are usually associated with the accretion disc’s apsidal precession (see e.g. Patterson et al. 2002).

The eclipsing cataclysmic variable V2051 Oph was discovered by Sanduleak (1972). Only few outbursts of the system have been reported (Warner & O’Donoghue 1987, Warner & Cropper 1983, Bateson 1980) with a maximal magnitude in B of about 13 in September 1980, possibly in April 1982 and in May 1984. It is possible that outbursts were missed due to the faintness of the system (Wenzel 1984). In May 1998 the first superoutburst was detected (Kiyota & Kato 1998).

We observed V2051 Oph during this superoutburst between maximum and near quiescence. This paper presents an analysis of these data using the Eclipse Mapping technique invented by Horne (1985).

## 2 The data

Our observing set consists of 10 eclipses taken in four nights between 23 and 27 May 1998 (interrupted by a cloudy night) during the superoutburst first reported by R. Stubbings through VSNET. The data were gathered by W. Offutt with the 60cm Ritchey-Chretien telescope at the W&B observatory in New Mexico/USA using a #21 Orange filter which covers R band and the long wavelength region of the V band, see Fig. 1. For reasons described in Section 2.1 the eclipse light curves were not averaged or phase binned, but left at the original phase resolution of 10 to 15 sec (about 0.0023 in phase). A log of the superoutburst data is shown in Table 1.

According to VSNET observers the superoutburst started on May 17 when the system rose to a magnitude of 13.5. It reached the maximum on May 18 with a magnitude of 11.7 and ended after May 27. Since this system in quiescence is usually too faint for amateur astronomers it is not possible to pin down the exact length of the outburst or even get a complete superoutburst light curve.

Our series therefore started during decline of the superoutburst and ended when the system got too faint for the telescope with the given filter. This means we did not cover the full decline, but only down to near the quiescent state.

For the phasing of the data we used the ephemeris of Echevarría & Alvarez (1993). Since the white dwarf is most certainly outshown by the accretion disc and due to the irregular eclipse profile (see following Section) we did not attempt to determine a new ephemeris from these superoutburst data.

### 2.1 The light curves

Fig. 2 shows the observed eclipse light curves. A prelimimary analysis of the superoutburst data was shown in Vrielmann & Offutt 1999. But in that case we averaged the light curves of each night, losing information about the short term variation of the disc.

As can be seen in Fig. 2, the out-of-eclipse levels and the eclipse profile change dramatically throughout the four nights. The out-of-eclipse level drops from about 35 mJy down to 10 mJy and the eclipse minimum appears at positive or negative phases. The shifts in eclipse minimum indicate that the disc brightness distribution is asymmetric during superoutburst, with occasionally the leading or the following lune of the disc being brighter.

The eclipse shape is sometimes symmetric (as in eclipses S01, S49, S50, S51 and S65) and in other cases very asymmetric with out-of-eclipse levels before eclipse being higher (S02,S03,S66) or lower (S33,S34). These general differences will appear in the reconstructions as distinctive features.

## 3 The analysis

### 3.1 Eclipse Mapping

The eclipse mapping technique was invented by Horne (1985) as a tomographic tool to reconstruct the intensity distribution in the accretion disc by fitting the eclipse profile. Thus, intensity maps of the accretion disc are derived using a maximum entropy algorithm (MEM, Skilling & Bryan 1984). The MEM algorithm serves to avoid the ambiguities of reconstructing a two dimensional map from a one-dimensional data set (the light curve) and introduces a certain degree of symmetry in the map111For the current mapping we used a solid arc smearing with smearing parameters 0.1$`R_{\mathrm{SS}L1}`$ in radius and 20 in azimuth at the disc edge.. Otherwise, a minimum of assumptions about the system are made as to minimally influence the results by premonitions.

The only assumptions are that the geometry of the system is known (e.g. the secondary star fills its Roche-lobe, the accretion disc is flat or of pre-determined 3-dimensional shape) and that the accretion disc is not changing during the course of the eclipse as any short term variation will lead to smeared out features in the reconstructions. While the first assumption is generally relatively well fullfilled, the latter assumption is usually not, since accretion discs show short-term fluctuations in brightness called flickering. In order to minimize this effect, usually a number of eclipse light curves are averaged.

In the present case we did not average the eclipse light curves, because we were interested in possible changes in the disc between eclipses, e.g. due to a precession of the disc. The precession may cause a change of the orientation of the elliptical disc between two consecutive eclipses. Therefore, we do not want to lose this information by averaging over the eclipse light curves of a specific night. This means, on the other hand, that we have to be careful not to over-interpret structures in the light curve which are due to flickering.

### 3.2 The system parameters

The inclination angle $`i=83.3^{}`$, the mass ratio $`q=0.19`$, the white dwarf mass $`_{wd}=0.78_{}`$ and radius $`R_{wd}=0.0244R_{}`$ as well as the scale parameter $`R_{\mathrm{SS}L1}/R_{}=0.42`$ were taken from Baptista et al. (1998). For the distance we used the Physical Parameter Eclipse Mapping estimate by Vrielmann, Stiening & Offutt (2002) of 146 pc.

### 3.3 The disc geometry

In a first attempt to reconstruct the outburst disc we used a geometrically infinitesimally thin disc in the orbital plane. However, the disc appeared to be bright in the back part (disc region opposite the seconary star) – showing a so-called back-to-front asymmetry –, indicating that the disc has a certain, non-negligible thickness. This can also be seen in Z Cha during outburst (Robinson, Wood & Wade 1999). Therefore, for the study of the outburst disc we use a 3-dimensional disc geometry.

We constructed a circular disc geometry with a fixed opening angle. At the edge we allowed emission to be reconstructed on a vertical ribbon around the disc. This geometrical model introduces two new parameters, the opening angle and the disc radius (or disc size). However, it turned out that the exact value of the opening angle did not influence the results dramatically. Similar results were found by Robinson et al. (1995) and Robinson (1999, private communication).

Since we used an inclination angle of $`i=83.3^{}`$, we arbitrarily chose a half opening angle of $`\alpha =6.7^{}`$ (which allows us to see just about all parts of the disc and the white dwarf) and a somewhat smaller angle of $`\alpha =4^{}`$. A comparison of the reconstructions showed that in almost all cases the entropy corresponding to the maps with the larger angle is somewhat larger than the entropies for the thinner disc. In the remaining cases (eclipses S01, S34, S51) the entropies of the reconstructions of both angles were almost identical. Therefore we chose the larger angle of $`\alpha =6.7^{}`$.

Using a fixed opening angle means introducing an edge of the disc in the form of a ribbon around the disc perimeter with 126 pixels onto which intensity could be reconstructed as well as onto the circular disc surface. This ribbon is set vertical to the orbital plane around the disc at the disc edge $`R_e`$, therefore we can identify the ribbon pixel uniquely by their azimuthal angle. We see at most half of this ribbon, which half depending on the orbital phase. During eclipse, a further part of this ribbon is occulted by the secondary. A similarly flared disc geometry was used by Bobinger et al. (1997).

Since the 3:1 resonance radius in V2051 Oph’s accretion disc lies at $`R_{3:1}=0.68R_{\mathrm{SS}L1}`$, the superoutburst disc is expected to be larger than this. On the other hand, the Roche-lobe of the primary confines the disc to within about $`R_R=0.8R_{\mathrm{SS}L1}`$. We chose a disc edge at $`R_e=0.7R_{\mathrm{SS}L1}`$.

If the radius $`R_e`$ is chosen too small we should see this as enhanced emission from the outer regions of the disc; if it is too large, the reconstructions simply should not show much intensity in the outer regions of the disc. Therefore, since the range between $`R_{3:1}`$ and $`R_R`$ is small, the exact choice will not influence our results dramatically. However, because of the MEM algorithm used in Eclipse Mapping, the reconstructed brightness distribution in the disc – especially close to the disc edges – is severely smeared out (see also Appendix A) and will not be affected by a small difference in true to assumed radius.

The introduction of a ribbon allowed us to use the light curves as they were, without levelling the out-of-eclipse light curves. A flat disc geometry cannot cope with different out-of-eclipse levels, but in the flared disc approach intensity can be reconstructed onto the ribbon. For example an orbital hump can be reconstructed on the leading side of the ribbon, very well visible before eclipse and more or less invisible after. However, in the present case we are dealing with superoutburst data. We can therefore not distinguish between an orbital hump and a superhump – especially since our observations only cover the immediate eclipses. Nevertheless, we did not attempt to level out the out-of-eclipse light curves.

The use of a circular disc geometry seems to contradict the expectation for a precession of the disc with an elongated, elliptical disc shape. However, in order to avoid any assumption about the direction of the semi major axis in the disc, we use a circular disc geometry. In the case that the disc is truely elongated, we will face similar scenarios as described above for a too small or too large $`R_e`$, this time only applying to parts of the disc, i.e. we would either see a lack of emission or enhancement close to the edge of two regions $`180^{}`$ apart from each other.

To identify the location of structures in the disc we use the radius and azimuth. Azimuth $`0^{}`$ indicates the direction from the white dwarf to the secondary, negative azimuths correspond to the following lune, positive angles to the leading lune of the accretion disc (the side where the accretion stream hits the disc).

### 3.4 The maps

The Eclipse Mapping fits to the light curves (Fig. 2) are good, we reached a $`\chi ^2`$ of 1 for all light curves except for light curves S50, S65 and S66 where we aimed only at a $`\chi ^2`$ of 1.5 instead of 1. Aiming at a lower $`\chi ^2`$ would have led to artificial structures and noise in the maps.

Since the outburst disc is usually optically thick, the reconstructed brightness distribution can directly be translated into a map of the brightness temperature and treated as good approximations of the true disc temperatures.

The disc reconstructions are displayed in Fig. 3. The inner disc appears to be very hot on the first night, with temperatures up to 100 000 K in the disc centre and around 50 000 K in the intermediate disc ($`0.2R_{\mathrm{SS}L1}`$). The temperature drops in the following nights to values around 30 000 K on May 25 and 26 and 20 000 K in the last night within a radius of about $`0.2R_{\mathrm{SS}L1}`$. These values are compatible with theoretical considerations (Smak 2000, Menou, Hameury & Stehle 1999).

In the intermediate disc the temperature either follows a steady state distribution (e.g. S02, S33, S49, S66) or is flat (S34, S51) or inbetween. The mass accretion rate in the former cases is about $`10^{18}`$ gs-1 in S02, $`3\times 10^{17}`$ gs-1 in S33, $`5\times 10^{17}`$ gs-1 in S49 and $`10^{17}`$ gs-1 in S66. The last value is about a magnitude above the quiescence value of $`10^{16}`$ gs-1 (Vrielmann et al. 2002).

Fig. 4 shows the temperature distribution along the ribbon. A comparison with Fig. 3 shows that the reconstructed temperature in this ribbon is significantly higher than in the disc edge and in some cases comparable to the temperature at a radius of about $`0.2R_{\mathrm{SS}L1}`$. As such feature does not appear in the test cases, it is unlikely to be a systematic effect caused by the method. It is possible that this indicates that the opening angle of the disc is chosen too small, thus concentrating the emission on a too small surface at the disc edge. Furthermore, this would mean that parts of the disc surface are eclipsed by the high disc edge. As the jump is strongest in those eclipses which show the strongest asymmetry (i.e. S02, S03, S33, S34, with the exception of S51) this might indicate an uneven height of the disc edge and a variablility of the flaring angle during the course of the decline.

While for each night the distribution is similar, large differences appear for different nights, reflecting a slow but significant change of the disc. In Section 4.2 we describe what we can learn from these plots.

The grey-scale plots in Fig. 5 give some indication of the emissivity and possibly shape of the disc. However, since the MEM algorithm introduces a certain amount of radial and azimuthal smearing (and therefore a certain amount of symmetry) and due to a fixed circular geometry we will not necessarily be able to directly reconstruct an elliptical disc. A comparison to the test case in Appendix A, however, shows that in all maps we can see an asymmetry possibly indicating an elliptical – definately non-circular – intensity distribution. The fact that in the maps the “front” instead of the “back” part is bright might indicate that the opening angle we chose is too large, however, the bright region is not always symmetric around azimuth $`0^{}`$.

An interesting feature is the arc like structure in S01 and S03 and possibly S49 (merged with the central part of the disc). Similar structures have been observed in IP Peg during decline from outburst (Vrielmann 1997).

### 3.5 The uneclipsed component

Table 2 and Fig. 6 give the reconstructed uneclipsed component in the #21 Orange filter for the observed outburst nights. The scatter of the uneclipsed flux values in each night is an indication for the error of these values of roughly half a mJy. The flux level drops dramatically during the four nights. Since the filters R and #21 Orange are similar, one can compare this flux with the uneclipsed quiescent flux in the R band of 0.66 mJy (see Vrielmann et al. 2002). This indicates that even in the fourth night the uneclipsed light was not yet at its quiescent level.

## 4 Discussion

### 4.1 The disc evolution during superoutburst

Assuming that in superoutburst the disc emits as an optically thick source, we derive brightness temperatures from observations in an orange filter (part of V and up to I). The reconstructed maps show that the disc cools down considerably during the four nights.

Only in the intermediate eclipses (i.e. S34, S49, S50, S51, S66) do we see indications for a cooling front as expected by theoretical calculations (e.g. Cannizzo 1993). Measuring the location of the cooling front in S34 and S49 we derive a speed of $`0.12`$ km s-1, measuring the speed of the cooling front between S51 and S66 we also derive a speed of about $`0.12`$ km s-1, whereas these two cooling fronts are not identical.

These speeds are comparable to the value Menou et al. (1999) found with numerical simulations and who predict a fraction of a km s-1. Furthermore, it is comparable with the speeds derived with Eclipse Mapping of $`0.43`$ to $`0.06`$ km s-1 for EX Dra (decreasing during outburst, Baptista & Catalàn 2001), of $`0.8`$ km s-1 for IP Peg (Bobinger et al. 1997) and $`0.7`$ to $`0.14`$ km s-1 for OY Car in superoutburst (Bruch, Beele & Baptista 1996).

Now, why do we see two distinctly different cooling fronts? The explanation can be found in numerical calculations by e.g. Menou et al. (2000) who find that an inward moving cooling front can reflect an outward moving heating front which again reflects into an inward moving cooling front. If we assume this scenario we derive a speed of $`+1.8`$ km s-1 for the reflected heating front. A speed of a few km s-1 is typical for heating fronts (Menou et al. 1999). Furthermore, the reflection of the heating front occurs at a disc radius $`R_{front}`$ of about $`4\times 10^9`$ cm, compatible with Menou et al.’s (2000) finding that they occur at $`R_{front}\genfrac{}{}{0pt}{}{<}{}10^{10}`$ cm in their simulations.

If such multiple reflections of cooling and heating fronts do occur in V2051 Oph, we should observe re-flares as predicted by numerical simulations (Menou et al. 2000). Unfortunately, the coverage of this object by amateur astronomers (VSNET, AAVSO) during the so far observed superoutbursts is too sparse to confirm or reject the appearance of re-flares. Detailed observations during the next super outburst are highly desirable.

### 4.2 The disc precession

During superoutburst, we see a partially strong, variable peak in the intensity distribution on the edge of the disc. It is probably caused by an elliptical, precessing disc: Whenever one of the vertices of the semimajor axis of the disc is close to the secondary, the disc will experience tidal forces leading to precession and local heating of the disc. The disc will react on the temperature increase by flaring up near the secondary. So, the location of the peak gives a rough estimate for the phase of the apsidal precession. In Appendix A we show that the location of a bright area (not to be confused with the bright spot) on the disc edge can be relatively well determined.

Looking at the distribution of light along the rim of the disc, we identify three structures present in most reconstructions: (a) a central peak close to azimuth $`0^{}`$ (facing the secondary); (b) a secondary peak in the profile of the central peak; and (c) a broad underlying brightness distribution with a maximum at 90 and/or -90. The latter peaks at the same side of the disc as the secondary peak is located. Similar features appear in the test (Appendix A).

The central peak is also present in the test cases and therefore most likely an artefact, although it could also be associated with the bright spot where the ballistic stream hits the accretion disc. Such a spot was not considered in the test case. The location of the bright spot is only marginally influenced by a apsidal motion of an elliptical, precessing disc. As apparent from the tests, the large scale variation on the disc rim is most likely an artefact.

If we measure the position of the second peak, we find a 10 prograde shift per eclipse or a 10-20 backwards shift per day. If this can be considered an indication for the apsidal motion of the disc, we arrive at a value for the precession period $`P_{prec}`$ of 51h to 54h, slightly larger than two days (therefore the apparent backward shift during consecutive nights). This would mean the disc is slightly larger than the 3:1 resonance radius $`R_{3:1}`$ of 0.68$`R_{\mathrm{SS}L1}`$ and the superhump period ($`P_{sh}=P_{prec}P_{orb}/(P_{prec}P_{orb}`$), Warner 1995) should be 2.5 min or 3% larger than the orbital period $`P_{orb}=1.4982`$ hours, i.e. $`P_{sh}`$(predicted) = 1.5422$`{}_{}{}^{h}\pm `$ 0.0013h. This is in very good agreement with the superhump period of 0.06423d = 1.54152h independently determined by Kiyota & Kato (1998) during the same super outburst.

This value for the precession period is very typical for SU UMa type dwarf novae, e.g. DV UMa (67h, Patterson et al. 2000), HT Cas (52h), Z Cha (51h), OY Car (65h), IR Gem (47h), AQ Eri (52h) to name but a few or even the nova V1974 Cyg with 44h (Retter, Leibowitz & Ofek 1997). Numerical calculations of super outbursts of dwarf novae with a similar mass ratio as V2051 Oph give on average a slightly larger ratio of $`P_{sh}/P_{orb}`$ of about 1.06 (Murray 1998) compared to our 1.03.

### 4.3 The mass accretion rate

The mass accretion rate $`\dot{}`$ in superoutburst appears very high, around $`10^{18}`$ gs-1 in our first night and varying throughout the disc. This is about two magnitudes larger than in quiescence (about $`10^{16}`$ gs-1, Vrielmann et al. 2002). Towards the end of our observing run (somewhat before the end of the outburst) the disc cools down and shows a mass accretion rate of about $`10^{17}`$ gs-1, still a magnitude higher than the quiescent level.

Only occasionally the disc appears to follow the steady state models and then only in intermediate disc regions between 0.1$`R_{\mathrm{SS}L1}`$ and 0.2$`R_{\mathrm{SS}L1}`$. At other times the distribution is rather flat. During the last eclipse, when the disc is near the quiescent level, indication for a bright spot can be seen at the expected location where the accretion stream hits the disc.

Since the disc undergoes dramatic changes during the course of an outburst with cooling fronts travelling through the material, we would not expect it to settle down completely into a steady state. However, the fact that parts of the disc do seem to appear in steady state suggests that the time scale for returning into the steady state is relatively short, of the order of a day or less.

According to Tylenda’s (1981) calculations, the mass accretion rate in superoutburst is too high and the disc is too small (outer radius at $`2.310^{10}`$ cm) to show optically thin areas – the disc should be everywhere optically thick. This agrees very well with our assumption to consider the brightness temperatures as good approximations of the gas temperatures.

### 4.4 The disc wind

The Eclipse Mapping shows that the uneclipsed component is decreasing in strength on decline from outburst, but never disappearing completely, as the analysis of the quiescence data show (Vrielmann et al. 2002). A similar variability has been found for EX Dra by Baptista & Catalán (1999) where the out-of-eclipse flux in R dropped during decline from outburst and practically disappeared at minimum light. At that time EX Dra’s disc is reduced to the immediate vicinity of the white dwarf. Baptista & Catalán suggest this variable part of the uneclipsed component is due to a disc wind and a chromosphere ejected by the inner part of the disc. The variation in these components is presumably due to a variable mass accretion rate.

Our geometrical model of the superoutburst disc is very simple; the disc may in fact be thicker, tilted, warped or only flare close to the secondary. The latter geometrical model would introduce too many free parameters (extent and location of flared part) which we cannot determine from our data. However, the uneclipsed component is a flux independent of orbital phase and therefore unlikely to originate in the (eclipsed parts of the) disc itself, but rather in a region above the disc. The easiest explanation is the same as for EX Dra: a variable disc wind triggered by the outburst.

## acknowledgments

We thank all amateur astronomers who are dedicating much of their spare time to continously monitor dwarf novae in order to detect outbursts. SV is funded by the South African Claude Harris Leon Foundation through a postdoctoral fellowship.

## Appendix A Test for reconstructing elliptical discs

In order to test how well an elliptical accretion disc can be reconstructed, we performed a test in which we constructed an artificial elliptical disc with an azimuthal angle of the semimajor axis of $`20^{}`$ and a bright area at the vertex of the semimajor axis closest to the secondary (see Fig. 7 (left)). This bright area is a simple simulation of the flaring of the disc near the secondary. The geometry of the disc surface is the same as that described in Section 3.3.

For this artificial disc we calculated the light curve (see Fig. 8) that is to be fitted with the Eclipse Mapping algorithm. Artificial noise with a signal-to-noise level of 100 has been added to the light curve. This so constructed light curve shows a clear asymmetry with different slopes for the ingress and egress, different out-of-eclipse levels before and after eclipse as well as a minimum at a phase $`0`$. These features are also apparent in the superoutburst light curves of V2051 Oph (Fig.2).

The application of the Eclipse Mapping algorithm to these light curves lead to the reconstructed maps in Fig. 7 (right). On the surface of the disc only the “back” part of the disc (around an azimuth of $`160^{}`$) shows some brightening, the front part (around $`20^{}`$) seems to be unchanged. This is probably due to a forshortening of the front part of the disc because of the relatively large opening angle. In any case, the poor reconstruction of the elliptical shape in this test case means that we cannot expect to see an elliptical shape of the disc emission in Eclipse maps of observed light curves.

However, on the edge of the disc (the “ribbon”) a brightening can be seen – even if much fainter than in the original map. In Fig. 9 we have plotted the light distribution along the disc edge for reconstructions using eight different orientation angles between $`0^{}`$ and $`35^{}`$. For larger angles the flaring will be negligible.

Fig. 9 shows that in the reconstructions there is a large scale brightness variation, often a peak near azimuth $`0^{}`$ and for disc orientation angles $`>10^{}`$ a peak that coincides well with the original position of the bright area on the disc edge. Although the reconstructed brightness distribution for a disc orientation angle of $`0^{}`$ is not similar to the original one, it is at least symmetric around azimuth $`0^{}`$. For a disc orientation angle of $`5^{}`$ the reconstruction near azimuth $`0^{}`$ is rather noisy and should be recognizable in real data as a disc with an orientation angle close to $`0^{}`$.

This test shows that a peak in the brightness distribution along the disc edge can be determined with a positional accurancy of $`5^{}`$.


text:
# Spectral variability of the nucleus of M33 in a Chandra /ACIS observation.

## 1 INTRODUCTION

M33 is a nearby (795 kpc, van den Bergh, 1991) Sc galaxy. Its nucleus is the brightest X-ray source in the Local Group: with a luminosity of $`10^{39}`$ erg/sec, it accounts for $`70\%`$ of the total luminosity of the galaxy in the 0.15-4.5 keV band (Trinchieri, Fabbiano & Peres, 1991). The nature of this source (M33-X8 in Trinchieri, Fabbiano & Peres, 1991) is still unclear. Because of its position, coincident with the optical center of the galaxy, it has been suggested that X8 could be a low-luminosity AGN (see, e.g. Markert & Rallis, 1983), but this interpretation has two main drawbacks: no activity has been found in other energy bands and an upper limit of only $`1500M_{}`$ has been put on the mass of the central black hole (Gebhardt et al., 2001). Dubus et al. (1997) have found an apparent periodicity of $`106`$ days in a set of ROSAT data with a six years coverage, that they ascribe to a “superorbital” modulation in a binary system with a giant star orbiting around a black hole. However, this finding has not been confirmed by the subsequent observations (see e.g. Parmar et al., 2001), that failed to match the phase prediction. Today, the most favoured hypothesis is that M33-X8 is an intermediate (or even stellar) mass black hole, likely in a binary system. Recent studies with ASCA and BeppoSAX data (Takano et al., 1994; Parmar et al., 2001) agree in describing its X-ray energy spectrum with a power-law with a high energy exponential cutoff with $`2`$ keV $`e`$-folding energy, or with a multicolor disk black-body (MCD), both consistent with the spectra of Galactic black-hole candidates and of ultra luminous sources (ULXs, see. e.g., Makishima et al., 2000).

In this work we examine a Chandra/ACIS-S observation of the nucleus of M33, and we discuss its variability and spectral characteristics.

## 2 THE OBSERVATION AND DATA REDUCTION

M33 was observed with Chandra/ACIS on July 6th 2001 (Obs Id: 2023), with an exposure time of 91.8 ksec. The observation was pointed to the star formation region NGC604, with the aimpoint on ACIS-I, therefore the nucleus of the galaxy is $`12^{}`$ off-axis and falls on ACIS-S3. Though the image of the nucleus is highly degraded, this proves to be an advantage in terms of spectral and temporal analysis, as the photons from the nuclear source are not affected by pile-up effects, a major problem if the source were on-axis.

Due to the degradation of the point spread function (PSF) with increasing off-axis angle, the source appears distorted and elongated with an elliptical shape. The source counts were thus extracted from an elliptical region centered on the position of M33-X8 as determined from ROSAT/HRI data by Schulman & Bregman (1995), i.e. RA=$`01^h33^m50.92^s`$, Dec=$`30^{}39^{}36.7^{\prime \prime }`$ (J2000). The axes of the ellipse were determined examining the radial profiles along the two main direction of the PSF and chosing the radii that enclose 99% of the emission in each direction, namely 40 and 25 arcsec, with a major axis position angle of $`53^{}`$. The background was extracted from a concentric elliptical annular region with 110 and 70 arcsec major and minor semi-axes respectively, excluding the portions of the annulus falling outside the chip.

The data analysis was performed using different tools in the CIAO 2.1 software package, the imaging application DS9 v.2.1, the timing analysis software package Xronos v.5.18 and the spectral analysis software XSPEC v11.2.

## 3 TEMPORAL ANALYSIS

The observation is continuous with a time coverage of 91.8 ksec, an average count rate of $`1.9`$ counts/sec in the 0.2-8.0 keV band and a time resolution of 3.24 s, that allows for a detailed temporal analysis with moderate resolution. The light-curve of M33-X8 has been binned in 40 intervals, each 2250 s long. The source is evidently variable: the light curves show a significant variability of more than 10% from maximum to minimum (Figure 1\- Panel A). Panels B, C, D in Figure 1 show the light-curve in three energy sub-ranges (0.2-1.5 keV, 1.5-4.5 keV, 4.5-8.0 keV), chosen to search for energy-dependent variability: the soft band samples the spectral region that hosts most of the absorption edges and spectral lines, the medium and hard bands cover the continuum and their ratio measures the slope of the continuum itself. We applied to each band (including the full band) various tests: the Kolmogorov-Smirnof test (see, e.g., Conover, 1971), the $`\chi ^2`$ test, the Collura test (a phase averaged $`\chi ^2`$ test plus an interpretation with a pulse model, Collura et al., 1987). These tests (see Table 1) confirm the presence of variability with high statistical significance in the full band and in the soft and medium bands. The Collura test indicates a variability timescale of $`5000`$ sec for the total, soft and medium bands.

The hardness ratio plotted in Figure 1 (Panel E) is calculated as $`HR=\frac{HardSoft}{Hard+Soft}`$ where Hard (H) is the number of counts in the 1.5-8.0 keV band and Soft (S) is the number of counts in the 0.2-1.5 keV band. This curve strongly suggests that some spectral variability occurred during the observation. A $`\chi ^2`$ test against the hypothesis of constant value (equal to the average value) gives $`\chi ^2/\nu =57.5/37`$, with a 2% probability that the hardness ratio is constant. As suggested by Figure 1 and, more strongly, by Figure 2, there is a correlation between the medium + hard band flux and the hardness ratio. This suggests that the variability is stronger in these bands than in the soft one. We find no evidence of correlations between total flux and hardness ratio (Figure 2, bottom panel), and X-ray colors with changing flux (Figure 3). Colors are defined as R1= Medium/Soft and R2=Hard/Medium, with Soft, Medium and Hard energy bands as in Table 1 (see also Figure 1).

A power spectrum analysis performed on the photon counts with a 3.24 s resolution in the frequency range $`10^50.3`$ Hz yields a flat frequency distribution that can be fit with a single power law with a slope of $`4\times 10^2`$ with $`\chi ^2/\nu =64.7/75`$, and is consistent within $`2\sigma `$ with the Poisson noise in the frequency range sampled by the observation.

## 4 SPECTRAL ANALYSIS

### 4.1 Integrated spectrum

The source spectrum appears to be affected by a PHA to PI conversion problem111see http://asc.harvard.edu/cal/Links/Acis/acis/Cal\_projects/index.html (note added by Alexey Vikhlinin on 3/1/2002) that causes spurious emission and absorption features between $`1.5`$ and $`2.5`$ keV. Pending resolution of this problem, the spectral analysis has been conducted excluding the energy range affected by this problem. This left us with two data ranges: 0.3-1.5 keV and 2.5-8.0 keV, that were analyzed simultaneously. Moreover, the spectrum has been corrected for the quantum efficiency degradation of ACIS at low energies, using the CIAO tool corrarf. We find that the total spectrum can be described with a three components model: a powerlaw with $`\mathrm{\Gamma }2.0`$, a MCD with kT$`{}_{int}{}^{}1.0`$ and a Raymond-Smith model (RS) with kT$`0.18`$ and solar abundance. The residuals, however, still show an excess at $`1`$ keV, that can be fitted with a gaussian at $`0.9`$ keV with $`\sigma 0.12`$ keV. After this new component, the reduced $`\chi ^2`$ becomes 1.27, with an improvement of 82 over 3 additional d.o.f., yielding an F-test probability lower than 0.1%. This bump could be the signature of an overabundance of Fe and/or Ne in an environment of warm matter. We thus tried to model it by allowing the Fe and Ne abundance to vary in the RS model, but the model was not able to reproduce it: the best result has been obtained having the Ne abundance four times higher than the solar value, but the residuals still show a narrow gaussian-like feature at 1 keV that cannot be accounted for by varying the abundance of other elements. Thus, the feature cannot be associated unambiguosly to any of the continuum components. Figure 4 shows the best fit model residuals, including also the range that has been excluded from the fit. It also shows that there is no evidence of any systematic excess in the residuals, thus the value of the reduced $`\chi ^2`$ can be entirely ascribed to the statistic scattering of the data. The full set of best fit parameter values is reported in Table 2. The total flux is $`F_{0.38.0keV}=2.35\times 10^{11}`$ erg sec-1cm-2, that, at the distance of M33, gives a luminosity of $`1.5\times 10^{39}`$ erg s-1.

The flux of the RS component is roughly consistent with the value of the 20” extended emission reported by Schulman & Bregman (1995) with ROSAT/HRI. To investigate this possibility further, we extracted the spectrum of the circumnuclear region of M33 X-8 from one of the publicly archived on-axis observations of this source (Obs Id: 786, Exposure time: 50 ksec), using the same elliptical region described in Section 2 but excluding a circle with radius 4.92” (10pixels) centered on the nuclear source (heavily piled up in this observation). The resulting spectrum ($`3200`$ counts) does not match the RS in either shape or intensity. It can be well described by a power-law with $`\mathrm{\Gamma }=1.4`$, with a flux of $`4\times 10^{13}`$ erg sec-1cm-2in the 0.3-8.0 keV band ($`L_{0.38.0}=2.5\times 10^{37}`$ erg/s), $``$1/4 of the RS flux. So, most of the RS component originates from a region with a radius smaller than $`5`$”, corresponding to $`17`$ pc at the distance of M33.

### 4.2 Spectral variability

Figure 1-E suggests spectral variability. In order to investigate further this variability, we have divided the data into three phases, according to the value of the hardness ratio. We define the “HR-low” state as the one where the hardness ratio (HR) is lower than -0.22, the “HR-medium” state as the one with -0.22$`<`$HR$`<`$-0.18 and the “HR-high” state as the one where HR$`>`$-0.18. We have selected three spectra according to this criterion and we have analyzed them separately, using the same spectral component as in Section 4.1. The results are in Table 2. We find that the Gaussian component is not needed for the HR-low and HR-high states. This is not surprising, as each of these spectra includes only $`10\%`$ of the total counts of the source, and the gaussian component is relatively faint. Examining the other two spectral components we find that: a) the power-law photon index $`\mathrm{\Gamma }`$ does not change significantly in the three spectra, while the flux decreases from the HR-low to the HR-medium and HR-high states; b) the errors of the MCD components in the HR-low state are rather large, however there is evidence of the disk kT increasing from the medium to the hard state; c) we do not detect any variability in the RS, consistent with possibly extended emission from hot ISM. Figure 5 shows the three spectra, deconvolved from the instrumental response, with the best fit model components. From this figure we see a progression of larger MCD contribution (relatively to the power-law) going from the HR-low to the HR-high spectral states.

## 5 DISCUSSION

The light curve of M33-X8 shows a significant low-amplitude variability, with a $`10\%`$ variation on a timescale of a few thousands of seconds. This timescale is particularly evident both in the full band light curve, where the source count rate drops and suddenly goes back to the original value in $`3\times 10^4s`$, and in the medium band light curve, where the count rate oscillates a few times around its average value of $`0.71`$ c/s. This result strongly confirms the variability on a 3000 s timescale found by Peres et al. (1989) in the Einstein data of the nucleus of M33. However, the power spectrum is flat, and does not evidence any periodical or quasi-periodical oscillation. The hardness ratio is also variable, showing evidence of spectral variability, that appears to be correlated mostly with the emission at energy $`1`$ keV.

The spectral analysis reveals the presence of four components. Two of these, a Raymond-Smith component and a broad “Gaussian” line feature contribute relatively little to the total spectrum, and do not present any compelling evidence of variability. We will not discuss them any further. The two main components are a strong $`\mathrm{\Gamma }2`$ power-law and a kT$`=0.91.3`$ keV multicolor disk component (MCD). The parameters of the power-law and of the MCD are in good agreement with the results found by Parmar et al. (2001), Takano et al. (1994), Dubus & Rutledge (2002). However, their contribution to the total spectrum is variable, as showed by the HR-selected spectra (Table 2): the HR-low phase is characterized by a brighter power-low and a fainter disk. While we do not find any further variation of the power-law flux from the HR-medium to the HR-high state, the MCD flux seems to increase with increasing kT. This result may suggest that the variability is linked with fluctuations in the accretion rate.

For the MCD component, the lower limit of 37 km to the inner radius of the disk (Table 2; the limit is due to the presence of the $`\sqrt{\mathrm{cos}\theta }`$ factor accounting for the inclination angle of the disk) implies a lower limit to the mass of the compact object of $`4.2M_{}`$. The Eddington limit for this mass is L$`{}_{Edd}{}^{}=6.3\times 10^{38}`$ erg/sec, that is consistent with the power-law + disk emission of the source. This is in agreement with the conclusion of Parmar et al. (2001), Takano et al. (1994) and Dubus & Rutledge (2002) that M33 X-8 is most likely a normal binary system containing a stellar mass black hole that happens to be at the center of the galaxy, rather than an AGN. Although not extreme, the spectral variability of M33 X-8 is reminiscent of what has been observed in black hole binaries. In the Galactic black hole binaries we usually observe a “high/soft state”, dominated by a hot disk, sometimes associated with a faint power-law tail, and a “low/hard state” dominated by a power-law spectrum, with a colder disk component (see, e.g., Esin, Mc Clintock and Narayan, 1997). The soft state is explained with high accretion rate, with the accreting disk going down toward the last stable orbit, while the hard state is associated to a lower accretion rate that causes the disruption of the inner part of the disk and the presence of a comptonized corona wrapping the central source. In the case of M33 X-8 both disk and power-law component are always present, but the power-law dominated HR-low state corresponds with a lower flux in the medium energy band (1.5-4.5 keV). With increasing medium band flux the spectrum evolves to be dominated by the disk component.

The source that resembles more M33 X-8 is probably LMC X-3 (Nowak et al., 2001), a black hole candidate with a probable mass of $`9M_{}`$ and a luminosity of up to $`4\times 10^{38}`$ erg/s. M33 X-8 presents many similarities with LMC X-3: the thermal spectrum, that can be modelled with a MCD (with kT$`0.8`$ in LMC X-3 and kT$`1.2`$ in M33 X-8) and a steep power-law; the presence of variability on timescales of some ten thousands of seconds; a long timescale periodicity (116 d for M33 X-8 and 99 or 198 d for LMC X-3, see Dubus et al., 1997 and Cowley et al., 1991) that appears to vary or even disappear in a few years (Nowak et al., 2001; Parmar et al., 2001). Also, LMC X-3 is characterized by a lack of fast variability (on timescales $`1`$ s), but this feature cannot be investigated in M33 X-8, due to the lack of data with high enough temporal resolution. However the power-law component is relatively more prominent in M33 X-8.

Concluding, the Chandra spectral variability observations of the nucleus of M33 (X-8) strongly agree with the suggestion that this source is an accreting $`5`$ M black hole with a variable accretion disk, comparable to Galactic black hole binaries.

We thank M. Krauss and F. Nicastro for their help in solving the calibration problems, J. MacDowell for interesting discussions and A. Kong, D. Schwartz and J. Drake for useful comments on the manuscript. This research has been supported by CNAA and by NASA contract NAS 8-39073 (CXC).


text:
# Elliptic Selberg Integrals and Conformal Blocks

## 1. Introduction

The Selberg integral is the integral

$$B_p(\alpha ,\beta ,\gamma )=_{\mathrm{\Delta }_p}\underset{j=1}{\overset{p}{}}t_j^{\alpha 1}(1t_j)^{\beta 1}\underset{0j<k1}{}(t_jt_k)^{2\gamma }dt_1\mathrm{}dt_p,$$

where $`\mathrm{\Delta }_p=\{t^p|\mathrm{\hspace{0.17em}0}t_p\mathrm{}t_11\}.`$ The Selberg integral is a generalization of the beta function. It can be calculated explicitly,

$$B_p(\alpha ,\beta ,\gamma )=\frac{1}{p!}\underset{j=0}{\overset{p1}{}}\frac{\mathrm{\Gamma }(1+\gamma +j\gamma )\mathrm{\Gamma }(\alpha +j\gamma )\mathrm{\Gamma }(\beta +j\gamma )}{\mathrm{\Gamma }(1+\gamma )\mathrm{\Gamma }(\alpha +\beta +(p+j1)\gamma )}.$$

The Selberg integral has many applications, see \[A1, A2, As, D, DF1, DF2, M, S\]. In this paper, we present elliptic versions of the Selberg integral.

## 2. Conformal Blocks on the Torus

Let $`\tau `$ be such that $`\mathrm{Im}\tau >0`$. Let $`\kappa `$ and $`p`$ be non-negative integers satisfying $`\kappa 2p+2`$. The KZB-heat equation is the partial differential equation

(1)
$$2\pi i\kappa \frac{u}{\tau }(\lambda ,\tau )=\frac{^2u}{\lambda ^2}(\lambda ,\tau )+p(p+1)\rho ^{}(\lambda ,\tau )u(\lambda ,\tau ).$$

Here, the prime denotes the derivative with respect to the first argument, and $`\rho `$ is defined in terms of the first Jacobi theta function,

$$\vartheta _1(\lambda ,\tau )=2q^{\frac{1}{8}}\mathrm{sin}(\pi \lambda )\underset{j=1}{\overset{\mathrm{}}{}}(1q^je^{2\pi i\lambda })(1q^je^{2\pi i\lambda })(1q^j),\rho (\lambda ,\tau )=\frac{\vartheta _1^{}(\lambda ,\tau )}{\vartheta _1(\lambda ,\tau )},$$

where $`q=e^{2\pi i\tau }`$. Holomorphic solutions of the KZB-heat equation with the properties,

1. $`u(\lambda +2,\tau )=u(\lambda ,\tau )`$,
2. $`u(\lambda +2\tau ,\tau )=e^{2\pi i\kappa (\lambda +\tau )}u(\lambda ,\tau )`$,
3. $`u(\lambda ,\tau )=(1)^{p+1}u(\lambda ,\tau )`$,
4. $`u(\lambda ,\tau )=𝒪((\lambda mn\tau )^{p+1})`$ as $`\lambda m+n\tau `$ for any $`m,n`$

are called conformal blocks (or elliptic hypergeometric functions) associated with the family of elliptic curves $`/+\tau `$ with the marked point $`z=0`$ and the irreducible $`sl_2`$ representation of dimension $`2p+1`$. It is known that the space of conformal blocks has dimension $`\kappa 2p1`$.

## 3. Integral Representations of Conformal Blocks

Introduce special functions

$$\sigma _\lambda (t,\tau )=\frac{\vartheta _1(\lambda t,\tau )\vartheta _1^{}(0,\tau )}{\vartheta _1(\lambda ,\tau )\vartheta _1(t,\tau )},E(t,\tau )=\frac{\vartheta _1(t,\tau )}{\vartheta _1^{}(0,\tau )}.$$

Consider the theta functions

$$\theta _{\kappa ,n}(\lambda ,\tau )=\underset{j}{}e^{2\pi i\kappa (j+\frac{n}{2\kappa })^2\tau +2\pi i\kappa (j+\frac{n}{2\kappa })\lambda },n/2\kappa .$$

They form a basis of the space of theta functions of level $`\kappa `$. They satisfy the equations

$$\theta _{\kappa ,n}(\lambda +1,\tau )=(1)^n\theta _{\kappa ,n}(\lambda ,\tau ),\theta _{\kappa ,n}(\lambda +\tau ,\tau )=e^{\pi i\kappa (\lambda +\frac{\tau }{2})}\theta _{\kappa ,n+\kappa }(\lambda ,\tau )$$

and have the modular properties
$$\theta _{\kappa ,n}(\lambda ,\tau +1)=e^{\pi i\frac{n^2}{2\kappa }}\theta _{\kappa ,n}(\lambda ,\tau ),\theta _{\kappa ,n}(\frac{\lambda }{\tau },\frac{1}{\tau })=\sqrt{\frac{i\tau }{2\kappa }}e^{\pi i\kappa \frac{\lambda ^2}{2\tau }}\underset{m=0}{\overset{2\kappa 1}{}}e^{\pi i\frac{mn}{\kappa }}\theta _{\kappa ,m}(\lambda ,\tau ),$$

where $`|\text{arg}(i\tau )|<\pi /2`$. Let $`\theta _{\kappa ,n}^s`$ denote the symmetrization of $`\theta _{\kappa ,n}`$ with respect to $`\lambda `$, $`\theta _{\kappa ,n}^s(\lambda ,\tau )=\theta _{\kappa ,n}(\lambda ,\tau )+\theta _{\kappa ,n}(\lambda ,\tau )`$.

Define $`u_{\kappa ,n}`$ by

$$u_{\kappa ,n}(\lambda ,\tau )=u_{p,\kappa ,n}(\lambda ,\tau )=J_{p,\kappa ,n}(\lambda ,\tau )+(1)^{p+1}J_{p,\kappa ,n}(\lambda ,\tau ),$$

where

$$\begin{array}{c}J_{p,\kappa ,n}(\lambda ,\tau )=_{\mathrm{\Delta }_p}\underset{j=1}{\overset{p}{}}E(t_j,\tau )^{\frac{2p}{\kappa }}\underset{1j<kp}{}E(t_jt_k,\tau )^{\frac{2}{\kappa }}\times \hfill \\ \hfill \underset{j=1}{\overset{p}{}}\sigma _\lambda (t_j,\tau )\theta _{\kappa ,n}(\lambda +\frac{2}{\kappa }\underset{j=1}{\overset{p}{}}t_j,\tau )dt_1\mathrm{}dt_p.\end{array}$$

The branch of the logarithm is chosen in such a way that arg $`(E(t,\tau ))0`$ as $`t0^+`$, and the integral is understood as the analytic continuation from the region where all of the exponents in the integrand have positive real parts.

###### Theorem 3.1.

\[FV1\] For all $`n`$, the integrals $`u_{\kappa ,n}(\lambda ,\tau )`$ are solutions of the KZB-heat equation having the properties (i)-(iv).

###### Theorem 3.2.

\[FSV2\] We have

1. $`u_{\kappa ,n}=u_{\kappa ,n+2\kappa }`$ and $`u_{\kappa ,n}=e^{2\pi ipn/\kappa }u_{\kappa ,n}`$.
2. The set $`\{u_{\kappa ,n}(\lambda ,\tau )|n=p+1,\mathrm{},\kappa p1\}`$ is a basis for the space of conformal blocks. The integrals $`u_{\kappa ,n}`$ are identically zero for all other values of $`n`$ in the interval from $`0`$ to $`\kappa `$.

## 4. Transformations acting on the space of conformal blocks

Introduce four transformations $`A`$, $`B`$, $`T`$, and $`S`$ defined by

$$Au(\lambda ,\tau )=u(\lambda +1,\tau ),Bu(\lambda ,\tau )=e^{\pi i\kappa (\lambda +\frac{\tau }{2})}u(\lambda +\tau ,\tau ),$$

$$Tu(\lambda ,\tau )=u(\lambda ,\tau +1),Su(\lambda ,\tau )=e^{\pi i\kappa \frac{\lambda ^2}{2\tau }}\tau ^{\frac{1}{2}\frac{p(p+1)}{\kappa }}u(\frac{\lambda }{\tau },\frac{1}{\tau }),$$

where we fix $`\mathrm{arg}\tau (0,\pi )`$.

###### Proposition 4.1.

If $`u(\lambda ,\tau )`$ is a solution of the KZB-heat equation, then $`Au(\lambda ,\tau )`$, $`Bu(\lambda ,\tau )`$, $`Tu(\lambda ,\tau )`$, and $`Su(\lambda ,\tau )`$ are solutions too. Moreover, the transformations $`A`$, $`B`$, $`T`$ and $`S`$ preserve the properties (i)-(iv).

The proofs that $`T`$ and $`S`$ preserve the space of conformal blocks are given in \[EK\]. The proofs that $`A`$ and $`B`$ also preserve this space are straightforward and follow from the equations

$$\vartheta _1(\lambda +1,\tau )=\vartheta _1(\lambda ,\tau ),\vartheta _1(\lambda +\tau ,\tau )=e^{\pi i(2\lambda +\tau )}\vartheta _1(\lambda ,\tau ).$$

###### Lemma 4.2.

Restricted to the space of conformal blocks, the transformations $`A`$, $`B`$, $`T`$, and $`S`$ satisfy the relations

$$A^2=I,B^2=I,S^2=(1)^pie^{\pi i\frac{p(p+1)}{\kappa }}I,(ST)^3=(1)^pie^{\pi i\frac{p(p+1)}{\kappa }}I,$$

$$SAS^1=B,AB=(1)^\kappa BA,TB=i^\kappa BAT,$$

where $`I`$ denotes the identity transformation.∎

###### Lemma 4.3.

We have

$$Au_{\kappa ,n}(\lambda ,\tau )=(1)^nu_{\kappa ,n}(\lambda ,\tau ),Bu_{\kappa ,n}(\lambda ,\tau )=e^{2\pi i\frac{pn}{\kappa }}u_{\kappa ,\kappa n}(\lambda ,\tau ).\mathit{}$$

Let $`(t_{m,n})`$ and $`(s_{m,n})`$ be the matrices of the transformations $`T`$ and $`S`$, respectively, with respect to the basis $`\{u_{\kappa ,n}(\lambda ,\tau )|p+1n\kappa p1\}`$, namely, $`Tu_{\kappa ,n}=_{m=p+1}^{\kappa p1}t_{m,n}u_{\kappa ,m},Su_{\kappa ,n}=_{m=p+1}^{\kappa p1}s_{m,n}u_{\kappa ,m}`$. In Theorem 4.4, we give formulas for the matrices of $`T`$ and $`S`$ in terms of Macdonald polynomials of type $`A_1`$.

The Macdonald polynomials \[Ma\] of type $`A_1`$ are $`x`$-even polynomials in terms of $`ϵ^{mx}`$, where $`m`$. They depend on two parameters $`k`$ and $`n`$, where $`k`$ and $`n`$ are non-negative integers. They are defined by the conditions:

1. $`P_n^{(k)}(x)=ϵ^{nx}+ϵ^{nx}+`$ lower order terms, except for $`P_0^{(k)}(x)=1`$,
2. $`P_m^{(k)},P_n^{(k)}=0`$ for $`mn`$, where

$$f,g=\frac{1}{2}\text{Const Term}\left(fg\underset{j=0}{\overset{k1}{}}\left(1ϵ^{2(j+x)}\right)\left(1ϵ^{2(jx)}\right)\right).$$

###### Theorem 4.4.

\[FSV2\] Let $`ϵ=e^{\pi i/\kappa }`$. For $`p+1m,n\kappa p1`$, we have

$$t_{m,n}=ϵ^{\frac{n^2}{2}}\delta _{mn},$$

$$s_{m,n}=\frac{e^{\frac{\pi i}{4}}}{\sqrt{2\kappa }}ϵ^{p(nm)\frac{p(p+1)}{2}}(ϵ^mϵ^m)\left(\underset{j=1}{\overset{p}{}}(ϵ^{n+j}ϵ^{nj})\right)P_{np1}^{(p+1)}(m),$$

where $`\delta _{mn}=1`$ for $`m=n`$ and $`0`$ otherwise.

## 5. Integral identities

To formulate our main result, we need functions $`\eta `$, $`\varphi _1`$, $`\varphi _2`$, and $`\varphi _3`$. The Dedekind $`\eta `$-function is the function $`\eta (\tau )=q^{1/24}_{j=1}^{\mathrm{}}(1q^j)`$. We have $`\vartheta _1^{}(0,\tau )=2\pi \eta ^3(\tau )`$. Consider the functions \[W\]

$$\varphi _1(\tau )=\frac{\eta (\tau )^2}{\eta (\frac{\tau }{2})\eta (2\tau )}=q^{\frac{1}{48}}\underset{j=1}{\overset{\mathrm{}}{}}\left(1+q^{j\frac{1}{2}}\right),\varphi _2(\tau )=\frac{\eta (\frac{\tau }{2})}{\eta (\tau )}=q^{\frac{1}{48}}\underset{j=1}{\overset{\mathrm{}}{}}\left(1q^{j\frac{1}{2}}\right),$$

$$\varphi _3(\tau )=\sqrt{2}\frac{\eta (2\tau )}{\eta (\tau )}=\sqrt{2}q^{\frac{1}{24}}\underset{j=1}{\overset{\mathrm{}}{}}(1+q^j).$$

We have

$$\varphi _1\left(\frac{1}{\tau }\right)=\varphi _1(\tau ),\varphi _2\left(\frac{1}{\tau }\right)=\varphi _3(\tau ),\varphi _3\left(\frac{1}{\tau }\right)=\varphi _2(\tau ),$$

$$\varphi _1(\tau +1)=e^{\frac{\pi i}{24}}\varphi _2(\tau ),\varphi _2(\tau +1)=e^{\frac{\pi i}{24}}\varphi _1(\tau ),\varphi _3(\tau +1)=e^{\frac{\pi i}{12}}\varphi _3(\tau ).$$

Let

$$c_{\kappa ,n}=c_{p,\kappa ,n}=(2\pi )^{\frac{p(p+1)}{\kappa }}e^{\pi i\frac{p(3p1)}{2\kappa }}e^{\pi i\frac{p+1}{2}}B_p(\frac{n+1}{\kappa },\frac{2p}{\kappa },\frac{1}{\kappa })\underset{j=1}{\overset{p}{}}\left(1e^{2\pi i\frac{n+j}{\kappa }}\right).$$

Here, $`B_p(\alpha ,\beta ,\gamma )`$ is the Selberg integral.

###### Theorem 5.1.

We have ten series of identities,

(2)
$$u_{2p+2,p+1}(\lambda ,\tau )=c_{2p+2,p+1}\vartheta _1(\lambda ,\tau )^{p+1},$$

(3)
$$u_{2p+3,p+1}(\lambda ,\tau )=c_{2p+3,p+1}\eta (\tau )^{\frac{3(p+1)}{2p+3}}\vartheta _1^{p+1}(\lambda ,\tau )\theta _{1,0}(\lambda ,\tau ),$$

(4)
$$u_{2p+3,p+2}(\lambda ,\tau )=c_{2p+3,p+2}\eta (\tau )^{\frac{3(p+1)}{2p+3}}\vartheta _1^{p+1}(\lambda ,\tau )\theta _{1,1}(\lambda ,\tau ),$$

(5)
$$u_{2p+4,p+2}(\lambda ,\tau )=2^{\frac{2(p+1)}{2p+4}}c_{2p+4,p+2}\left(\varphi _3(\tau )\eta (\tau )^1\right)^{\frac{4(p+1)}{2p+4}}\vartheta _1^{p+1}(\lambda ,\tau )\theta _{2,1}^s(\lambda ,\tau ),$$

(6)
$$\begin{array}{c}u_{2p+4,p+1}(\lambda ,\tau )+(1)^{p+1}e^{2\pi i\frac{p(p+1)}{2p+4}}u_{2p+4,p+3}(\lambda ,\tau )=\hfill \\ \hfill c_{2p+4,p+1}\left(\varphi _2(\tau )\eta (\tau )^1\right)^{\frac{4(p+1)}{2p+4}}\vartheta _1^{p+1}(\lambda ,\tau )(\theta _{2,0}(\lambda ,\tau )\theta _{2,2}(\lambda ,\tau )),\end{array}$$

(7)
$$\begin{array}{c}u_{2p+4,p+1}(\lambda ,\tau )+(1)^pe^{2\pi i\frac{p(p+1)}{2p+4}}u_{2p+4,p+3}(\lambda ,\tau )=\hfill \\ \hfill c_{2p+4,p+1}\left(\varphi _1(\tau )\eta (\tau )^1\right)^{\frac{4(p+1)}{2p+4}}\vartheta _1^{p+1}(\lambda ,\tau )(\theta _{2,0}(\lambda ,\tau )+\theta _{2,2}(\lambda ,\tau )),\end{array}$$

(8)
$$\begin{array}{c}u_{2p+6,p+1}(\lambda ,\tau )+(1)^{p+1}e^{2\pi i\frac{p(p+1)}{2p+6}}u_{2p+6,p+5}(\lambda ,\tau )=\hfill \\ \hfill 2^{\frac{3(p+1)}{2p+6}}c_{2p+6,p+1}\left(\varphi _3(\tau )\eta (\tau )\right)^{\frac{6(p+1)}{2p+6}}\vartheta _1^{p+1}(\lambda ,\tau )(\theta _{4,0}(\lambda ,\tau )\theta _{4,4}(\lambda ,\tau )),\end{array}$$

(9)
$$\begin{array}{c}u_{2p+6,p+2}(\lambda ,\tau )+(1)^pe^{2\pi i\frac{p(p+2)}{2p+6}}u_{2p+6,p+4}(\lambda ,\tau )=\hfill \\ \hfill c_{2p+6,p+2}\left(\varphi _2(\tau )\eta (\tau )\right)^{\frac{6(p+1)}{2p+6}}\vartheta _1^{p+1}(\lambda ,\tau )(\theta _{4,1}^s(\lambda ,\tau )+\theta _{4,3}^s(\lambda ,\tau )),\end{array}$$

(10)
$$\begin{array}{c}u_{2p+6,p+2}(\lambda ,\tau )+(1)^{p+1}e^{2\pi i\frac{p(p+2)}{2p+6}}u_{2p+6,p+4}(\lambda ,\tau )=\hfill \\ \hfill c_{2p+6,p+2}\left(\varphi _1(\tau )\eta (\tau )\right)^{\frac{6(p+1)}{2p+6}}\vartheta _1^{p+1}(\lambda ,\tau )(\theta _{4,1}^s(\lambda ,\tau )\theta _{4,3}^s(\lambda ,\tau )),\end{array}$$

(11)
$$\begin{array}{c}u_{2p+8,p+2}(\lambda ,\tau )+(1)^{p+1}e^{2\pi i\frac{p(p+2)}{2p+8}}u_{2p+8,p+6}(\lambda ,\tau )=\hfill \\ \hfill c_{2p+8,p+2}\eta (\tau )^{\frac{8(p+1)}{2p+8}}\vartheta _1^{p+1}(\lambda ,\tau )(\theta _{6,1}^s(\lambda ,\tau )\theta _{6,5}^s(\lambda ,\tau )).\end{array}$$

The integrals in Theorem 5.1 are appropriately called the elliptic Selberg integrals. The identity (2) appears in \[FSV1\] and in \[FV1\] for $`p=1`$.

## 6. Differential equations

In Lemmas 6.16.3, let  denote the derivative with respect to $`\lambda `$, let $`\dot{}`$ denote the derivative with respect to $`\tau `$, and let $`v(\lambda ,\tau )=\vartheta _1^{p+1}(\lambda ,\tau )_{j=0}^{\kappa 2p2}c_j(\tau )\theta _{\kappa 2p2,j}^s(\lambda ,\tau )`$.

###### Lemma 6.1.

The function $`v(\lambda ,\tau )`$ is a solution of the KZB-heat equation if and only if the differential equation

(12)
$$\begin{array}{c}\frac{\kappa }{p+1}\underset{j=0}{\overset{\kappa 2p2}{}}\left(\frac{d}{d\tau }c_j\right)\theta _{\kappa 2p2,j}^s=(2p+2\kappa )\frac{\dot{\vartheta _1}}{\vartheta _1}\underset{j=0}{\overset{\kappa 2p2}{}}c_j\theta _{\kappa 2p2,j}^s\hfill \\ \hfill 2\underset{j=0}{\overset{\kappa 2p2}{}}c_j\dot{\theta }_{\kappa 2p2,j}^s+\frac{1}{\pi i}\frac{\vartheta _1^{}}{\vartheta _1}\underset{j=0}{\overset{\kappa 2p2}{}}c_j(\theta _{\kappa 2p2,j}^s)^{}\end{array}$$

holds.

The proof of Lemma 6.1 uses the identities

(13)
$$2\pi i(2p+2)\dot{(\vartheta _1^{p+1})}(\lambda ,\tau )=(\vartheta _1^{p+1})^{\prime \prime }(\lambda ,\tau )+p(p+1)\rho ^{}(\lambda ,\tau )\vartheta _1^{p+1}(\lambda ,\tau ),$$

(14)
$$2\pi i\kappa \dot{\theta }_{\kappa ,m}^s(\lambda ,\tau )=(\theta _{\kappa ,m}^s)^{\prime \prime }(\lambda ,\tau ).$$

###### Proof of Lemma 6.1.

Applying the differential operator $`2\pi i\kappa /\tau `$ to $`v(\lambda ,\tau )`$ gives

$$\begin{array}{c}2\pi i\kappa [(p+1)\dot{\vartheta _1}\vartheta _1^p\underset{j=0}{\overset{\kappa 2p2}{}}c_j\theta _{\kappa 2p2,j}^s+\vartheta _1^{p+1}\underset{j=0}{\overset{\kappa 2p2}{}}\left(\frac{d}{d\tau }c_j\right)\theta _{\kappa 2p2,j}^s\hfill \\ \hfill +\vartheta _1^{p+1}\underset{j=0}{\overset{\kappa 2p2}{}}c_j\dot{\theta }_{\kappa 2p2,j}^s].\end{array}$$

Applying the differential operator $`^2/\lambda ^2+p(p+1)\rho ^{}(\lambda ,\tau )`$ to $`v(\lambda ,\tau )`$ gives

$$\begin{array}{c}(\vartheta _1^{p+1})^{\prime \prime }\underset{j=0}{\overset{\kappa 2p2}{}}c_j\theta _{\kappa 2p2,j}^s+2(p+1)\vartheta _1^{}\vartheta _1^p\underset{j=0}{\overset{\kappa 2p2}{}}c_j(\theta _{\kappa 2p2,j}^s)^{}\hfill \\ \hfill +\vartheta _1^{p+1}\underset{j=0}{\overset{\kappa 2p2}{}}c_j(\theta _{\kappa 2p2,j}^s)^{\prime \prime }+p(p+1)\rho ^{}\vartheta _1^{p+1}\underset{j=0}{\overset{\kappa 2p2}{}}c_j\theta _{\kappa 2p2,j}^s.\end{array}$$

Applying (13) and (14), we obtain the result. ∎

###### Lemma 6.2.

If $`v(\lambda ,\tau )`$ is a solution of the KZB-heat equation, then the functions $`c_j(\tau )`$ satisfy the differential equation

$$\begin{array}{c}\frac{\kappa }{p+1}\underset{j=0}{\overset{\kappa 2p2}{}}\left(\frac{d}{d\tau }c_j(\tau )\right)\theta _{\kappa 2p2,j}(0,\tau )=\hfill \\ \hfill (2p+2\kappa )\left(\frac{d}{d\tau }\mathrm{ln}\vartheta _1^{}(0,\tau )\right)\underset{j=0}{\overset{\kappa 2p2}{}}c_j(\tau )\theta _{\kappa 2p2,j}(0,\tau )\\ \hfill +2(\kappa 2p3)\underset{j=0}{\overset{\kappa 2p2}{}}c_j(\tau )\frac{d}{d\tau }\theta _{\kappa 2p2,j}(0,\tau ).\end{array}$$

###### Proof.

For any fixed $`\lambda `$, equation (12) gives a differential equation for the functions $`c_j(\tau )`$. We take the limit of that equation as $`\lambda 0`$. In the ratio $`\dot{\vartheta _1}/\vartheta _1`$, both the numerator and denominator tend to zero, so the limit of this term as $`\lambda 0`$ is equal to the limit of the ratio of the derivatives of the numerator and the denominator. The limit of the ratio $`(_{j=0}^{\kappa 2p2}c_j(\theta _{\kappa 2p2,j}^s)^{})/\vartheta _1`$ is calculated in the same way, since each $`\theta _{\kappa 2p2,j}^s`$ is a symmetric function and therefore $`(\theta _{\kappa 2p2,j}^s)^{}(0,\tau )=0`$. Then the result follows from (14) . ∎

###### Lemma 6.3.

If $`v(\lambda ,\tau )`$ is a solution of the KZB-heat equation, then the functions $`c_j(\tau )`$ satisfy the differential equation

$$\begin{array}{c}\frac{\kappa }{p+1}\underset{j=0}{\overset{\kappa 2p2}{}}\left(\frac{d}{d\tau }c_j(\tau )\right)\theta _{\kappa 2p2,\kappa 2p2j}(0,\tau )=\hfill \\ \hfill (2p+2\kappa )\left(\frac{d}{d\tau }\mathrm{ln}\vartheta _1^{}(0,\tau )\right)\underset{j=0}{\overset{\kappa 2p2}{}}c_j(\tau )\theta _{\kappa 2p2,\kappa 2p2j}(0,\tau )\\ \hfill +2(\kappa 2p3)\underset{j=0}{\overset{\kappa 2p2}{}}c_j(\tau )\frac{d}{d\tau }\theta _{\kappa 2p2,\kappa 2p2j}(0,\tau ).\end{array}$$

###### Proof.

We take the limit of (12) as $`\lambda \tau `$. This limit is calculated in terms of the limit $`\lambda 0`$ using the formulas

$$\frac{}{z}\theta _1(z,\tau )|_{z=\lambda +\tau }=e^{2\pi i\lambda \pi i\tau }\left(2\pi i\theta _1(\lambda ,\tau )\theta _1^{}(\lambda ,\tau )\right),$$

$$\frac{}{z}\theta _1(\lambda +\tau ,z)|_{z=\tau }=e^{2\pi i\lambda \pi i\tau }\left(\pi i\theta _1(\lambda ,\tau )+\theta _1^{}(\lambda ,\tau )\dot{\theta }_1(\lambda ,\tau )\right),$$

$$\frac{}{z}\theta _{\kappa ,m}^s(z,\tau )|_{z=\lambda +\tau }=e^{\pi i\kappa \lambda \pi i\kappa \frac{\tau }{2}}\left(\pi i\kappa \theta _{\kappa ,\kappa m}^s(\lambda ,\tau )+(\theta _{\kappa ,\kappa m}^s)^{}(\lambda ,\tau )\right),$$

$$\frac{}{z}\theta _{\kappa ,m}^s(\lambda +\tau ,z)|_{z=\tau }=e^{\pi i\kappa \lambda \pi i\kappa \frac{\tau }{2}}\left(\pi i\frac{\kappa }{2}\theta _{\kappa ,\kappa m}^s(\lambda ,\tau )(\theta _{\kappa ,\kappa m}^s)^{}(\lambda ,\tau )+\dot{\theta }_{\kappa ,\kappa m}^s(\lambda ,\tau )\right).$$

It is straightforward to calculate the limit of the left hand side. Using the above formulas, the limit as $`\lambda \tau `$ of the right hand side is equal to the limit as $`\lambda 0`$ of the expression

$$\begin{array}{c}e^{\pi i(2p+2\kappa )(\lambda +\frac{\tau }{2})}((2p+2\kappa )\frac{\dot{\vartheta _1}}{\vartheta _1}\underset{j=0}{\overset{\kappa 2p2}{}}c_j\theta _{\kappa 2p2,\kappa 2p2j}^s\hfill \\ \hfill 2\underset{j=0}{\overset{\kappa 2p2}{}}c_j\dot{\theta }_{\kappa 2p2,\kappa 2p2j}^s+\frac{1}{\pi i}\frac{\vartheta _1^{}}{\vartheta _1}\underset{j=0}{\overset{\kappa 2p2}{}}c_j(\theta _{\kappa 2p2,\kappa 2p2j}^s)^{}).\end{array}$$

This limit is calculated using L’Hôpital’s rule. ∎

## 7. Identities for theta functions

In the next section, we give the proofs of the integral identities in Theorem 5.1. We will use the following results.

###### Lemma 7.1.

We have $`\theta _{2,1}^s(\lambda )=\vartheta _1(\lambda +1/2)`$.

Lemma 7.1 is proved by comparing the Fourier series expansions of the functions.

###### Corollary 7.2.

We have $`2\theta _{2,1}(0)=\eta (\tau )\varphi _3(\tau )^2.`$

###### Lemma 7.3.

Let

$$f_1(\tau )=\frac{\theta _{4,1}(0)\theta _{4,3}(0)}{\eta (\tau )},f_2(\tau )=\frac{\theta _{4,1}(0)+\theta _{4,3}(0)}{\eta (\tau )},f_3(\tau )=\frac{\theta _{4,0}(0)\theta _{4,4}(0)}{\sqrt{2}\eta (\tau )}.$$

Then $`f_1(\tau )=\varphi _1(\tau )^1,f_2(\tau )=\varphi _2(\tau )^1,f_3(\tau )=\varphi _3(\tau )^1.`$

The proof of Lemma 7.3 is based on the following result.

###### Lemma 7.4.

\[W\] Suppose $`g_1(\tau )`$, $`g_2(\tau )`$, and $`g_3(\tau )`$ are holomorphic functions on the upper half plane $`_+`$ satisfying the following conditions.

1. The functions $`g_1`$, $`g_2`$, and $`g_3`$ can be written in the forms

$$g_1(\tau )=q^{\frac{a}{48}}\underset{j=0}{\overset{\mathrm{}}{}}a_jq^{\frac{j}{2}},g_2(\tau )=q^{\frac{a}{48}}\underset{j=0}{\overset{\mathrm{}}{}}(1)^ja_jq^{\frac{j}{2}},g_3(\tau )=q^{\frac{a}{24}}\underset{j=0}{\overset{\mathrm{}}{}}b_jq^j,$$

where $`a`$ is an integer, $`a_j,b_j`$, and $`a_0=1`$.
2. The functions $`g_1`$, $`g_2`$, and $`g_3`$ have the modular properties

$$g_1\left(\frac{1}{\tau }\right)=g_1(\tau ),g_2\left(\frac{1}{\tau }\right)=g_3(\tau ),g_3\left(\frac{1}{\tau }\right)=g_2(\tau ).$$

Then $`g_i(\tau )=\varphi _i(\tau )^a`$, for $`i=1,2,3`$.

###### Proof of Proposition 7.3.

We have

$$\theta _{4,1}(0)=\underset{j}{}q^{4(j+\frac{1}{8})^2}=q^{\frac{1}{16}}\underset{j}{}q^{\frac{1}{2}(8j^2+2j)},\theta _{4,3}(0)=\underset{j}{}q^{4(j+\frac{3}{8})^2}=q^{\frac{1}{16}}\underset{j}{}q^{\frac{1}{2}(8j^2+6j+1)},$$

$$\theta _{4,0}(0)=\underset{j}{}q^{4j^2}=\underset{j}{}q^{(2j)^2},\theta _{4,4}(0)=\underset{j}{}q^{4(j+\frac{1}{2})^2}=\underset{j}{}q^{(2j+1)^2}.$$

Hence, the functions

$$\stackrel{~}{f}_1(\tau )=q^{\frac{1}{24}}(\theta _{4,1}(0)\theta _{4,3}(0)),\stackrel{~}{f}_2(\tau )=q^{\frac{1}{24}}(\theta _{4,1}(0)+\theta _{4,3}(0)),$$

$$\stackrel{~}{f}_3(\tau )=2^{\frac{1}{2}}q^{\frac{1}{24}}(\theta _{4,0}(0)\theta _{4,4}(0))$$

are holomorphic functions on $`_+`$ with the property P1 for $`a=1`$ and $`g_i=\stackrel{~}{f}_i`$. Let $`y(q)=_{j=1}^{\mathrm{}}c_jq^j`$ be defined by the condition $`1+y(q)=q^{1/24}\eta (\tau ).`$ Then for $`i=1,2,3`$,

$$f_i(\tau )=\frac{\stackrel{~}{f}_i(\tau )}{1+y(q)}=\stackrel{~}{f}_i(\tau )(1y(q)+y(q)^2+\mathrm{})$$

are holomorphic functions on $`_+`$ with the property P1 for $`a=1`$ and $`g_i=f_i`$. One checks that $`f_1`$, $`f_2`$ and $`f_3`$ have the property P2 using the modular properties of $`\theta _{4,n}(0)`$ and $`\eta (\tau )`$. ∎

###### Lemma 7.5.

We have $`\theta _{6,1}(0)\theta _{6,5}(0)=\eta (\tau ).`$

Lemma 7.5 is proved by comparing the infinite series expansions of the functions.

## 8. The proof of Theorem 5.1

### 8.1. Proof of (2)

For $`\kappa =2p+2`$, the space of conformal blocks is one-dimensional. The right hand side of (2) is a solution of (1) with the properties (i)-(iv) \[FV1\]. According to Theorem 3.1, the left hand side also has these properties. Thus the two functions are proportional. The coefficient of proportionality is calculated by comparing the leading terms of $`\vartheta _1^{p+1}`$ and $`u_{\kappa ,p+1}`$ in the limit as $`\tau i\mathrm{}`$. The leading term of $`\vartheta _1^{p+1}`$ is $`(i)^{p+1}q^{(p+1)/8}(e^{\pi i\lambda }e^{\pi i\lambda })^{p+1}`$. Let $`dt=dt_1\mathrm{}dt_p`$. The leading term of $`u_{\kappa ,p+1}`$ is

$$\begin{array}{c}_{\mathrm{\Delta }_p}\underset{j=1}{\overset{p}{}}\left(\frac{e^{\pi it_j}e^{\pi it_j}}{2\pi e^{\frac{\pi i}{2}}}\right)^{\frac{2p}{2p+2}1}\underset{1j<kp}{}\left(\frac{e^{\pi i(t_jt_k)}e^{\pi i(t_jt_k)}}{2\pi e^{\frac{\pi i}{2}}}\right)^{\frac{2}{2p+2}}\hfill \\ \hfill (\left(\underset{j=1}{\overset{p}{}}\frac{e^{\pi i(\lambda t_j)}e^{\pi i(\lambda t_j)}}{e^{\pi i\lambda }e^{\pi i\lambda }}\right)q^{\frac{(p+1)^2}{4(2p+2)}}e^{\pi i(p+1)(\lambda +\frac{2}{2p+2}_{j=1}^pt_j)}\\ \hfill +(1)^{p+1}\left(\underset{j=1}{\overset{p}{}}\frac{e^{\pi i(\lambda +t_j)}e^{\pi i(\lambda +t_j)}}{e^{\pi i\lambda }e^{\pi i\lambda }}\right)q^{\frac{(p+1)^2}{4(2p+2)}}e^{\pi i(p+1)(\lambda +\frac{2}{2p+2}_{j=1}^pt_j)})dt.\end{array}$$

The above expression is equal to

$$(2\pi e^{\frac{\pi i}{2}})^{\frac{p(p+1)}{2p+2}+p}e^{\pi i(\frac{2p^2}{2p+2}+p)}q^{\frac{p+1}{8}}(e^{\pi i\lambda }e^{\pi i\lambda })^p\underset{l=0}{\overset{p}{}}I_l(e^{\pi i(2l+1)\lambda }e^{\pi i(2l+1)\lambda }),$$

where
$$I_l=_{\mathrm{\Delta }_p}f_l(t_1,\mathrm{},t_p)\underset{j=1}{\overset{p}{}}e^{2\pi it_j(\frac{p+2}{2p+2}+\frac{1}{2})}(1e^{2\pi it_j})^{\frac{2p}{2p+2}1}\underset{1j<kp}{}(e^{2\pi it_j}e^{2\pi it_k})^{\frac{2}{2p+2}}dt,$$

for some function $`f_l`$, symmetric in the variables $`t_1,\mathrm{},t_p`$. Comparing the coefficients of $`e^{\pi i(p+1)\lambda }`$ in the leading terms, we find that $`u_{\kappa ,p+1}=i^{p+1}(2\pi e^{\frac{\pi i}{2}})^{\frac{p(p+1)}{2p+2}+p}e^{\pi i(\frac{2p^2}{2p+2}+p)}I_p\vartheta _1^{p+1}`$. To complete the proof, it remains to compute $`I_p`$. It is not difficult to show that $`f_p(t_1,\mathrm{},t_p)=_{j=1}^pe^{\pi it_j}`$. Let $`x_j=e^{2\pi it_j}`$. Let $`\stackrel{~}{\mathrm{\Delta }}_p`$ be the image of $`\mathrm{\Delta }_p`$ under the map $`t_jx_j`$. We have

$$I_p=(2\pi i)^p_{\stackrel{~}{\mathrm{\Delta }}_p}\underset{j=1}{\overset{p}{}}x_j^{\frac{p+2}{2p+2}1}(1x_j)^{\frac{2p}{2p+2}1}\underset{1j<kp}{}(x_jx_k)^{\frac{2}{2p+2}}dx.$$

Applying the Stokes theorem, we deform the contour $`\stackrel{~}{\mathrm{\Delta }}_p`$ to get

$$I_p=(2\pi i)^p\underset{j=1}{\overset{p}{}}(e^{2\pi i\frac{j+p+1}{2p+2}}1)_{\mathrm{\Delta }_p}\underset{j=1}{\overset{p}{}}x_j^{\frac{p+2}{2p+2}1}(1x_j)^{\frac{2p}{2p+2}1}\underset{1j<kp}{}(x_jx_k)^{\frac{2}{2p+2}}dx.$$

Observe that

$$_{\mathrm{\Delta }_p}\underset{j=1}{\overset{p}{}}x_j^{\frac{p+2}{2p+2}1}(1x_j)^{\frac{2p}{2p+2}1}\underset{1j<kp}{}(x_jx_k)^{\frac{2}{2p+2}}dx$$

is the Selberg integral $`B_p((p+2)/(2p+2),2p/(2p+2),1/(2p+2))`$. This completes the proof.

### 8.2. Proof of (3) and (4)

Let $`\kappa =2p+3`$. Then any solution of (1) with the properties (i)-(iv) has the form $`v(\lambda ,\tau )=\vartheta _1^{p+1}(\lambda ,\tau )(c_0(\tau )\theta _{1,0}(\lambda ,\tau )+c_1(\tau )\theta _{1,1}(\lambda ,\tau ))`$. Let $`A`$ be the transformation introduced in section $`4`$. By Proposition 4.1, $`Av(\lambda ,\tau )=(1)^{p+1}\vartheta _1^{p+1}(\lambda ,\tau )(c_0(\tau )\theta _{1,0}(\lambda ,\tau )c_1(\tau )\theta _{1,1}(\lambda ,\tau ))`$ is also a solution. Hence, for $`j=0`$ or $`1`$, the function $`v_j(\lambda ,\tau )=c_j(\tau )\vartheta _1^{p+1}(\lambda ,\tau )\theta _{1,j}(\lambda ,\tau )`$ gives a solution too. Moreover, $`Av_j=(1)^{p+1+j}v_j`$. By Theorem 3.2, the integrals $`u_{\kappa ,p+1}`$ and $`u_{\kappa ,p+2}`$ span the space of conformal blocks. By Lemma 4.3, $`Au_{\kappa ,p+1}=(1)^{p+1}u_{\kappa ,p+1}`$ and $`Au_{\kappa ,p+2}=(1)^pu_{\kappa ,p+2}`$. So for $`j=0`$ or $`1`$, the integral $`u_{\kappa ,p+1+j}`$ is proportional to $`v_j`$. By Lemma 6.2, $`c_j(\tau )`$ must satisfy the differential equation

$$\frac{\kappa }{p+1}\left(\frac{d}{d\tau }c_j(\tau )\right)\theta _{1,j}(0,\tau )=\left(\frac{d}{d\tau }\mathrm{ln}\vartheta _1^{}(0,\tau )\right)c_j(\tau )\theta _{1,j}(0,\tau ).$$

The function $`c_j(\tau )=\left((2\pi )^1\vartheta _1^{}(0,\tau )\right)^{(p+1)/\kappa }=\eta (\tau )^{3(p+1)/\kappa }`$ is a solution of this equation for $`j=0`$ and $`j=1`$. The coefficients of proportionality are computed in the limit as $`\tau i\mathrm{}`$, cf. the proof of (2). This completes the proof.

### 8.3. Proof of (5), (6), and (7)

Let $`\kappa =2p+4`$. Let $`v_j(\lambda ,\tau )=\vartheta _1^{p+1}(\lambda ,\tau )\theta _{2,j}^s(\lambda ,\tau ),\mathrm{\hspace{0.17em}0}j2`$. Any solution of (1) with the properties (i)-(iv) has the form $`v(\lambda ,\tau )=_{j=0}^2c_j(\tau )v_j(\lambda ,\tau )`$. Let $`A`$ be the transformation in section $`4`$. By Proposition 4.1, $`Av=_{j=0}^2(1)^{p+j+1}c_jv_j`$ is also a solution. Hence the function $`c_1v_1`$ gives a solution too. Moreover, it is an eigenvector of $`A`$ with eigenvalue $`(1)^p`$. By Theorem 3.2, the space of conformal blocks is three-dimensional with spanning set $`\{u_{\kappa ,n}|p+1np+3\}`$. According to Lemma 4.3, the eigenspace of $`A`$ corresponding to the eigenvalue $`(1)^p`$ is one-dimensional and is spanned by $`u_{\kappa ,p+2}`$. It follows that $`u_{\kappa ,p+2}`$ is proportional to $`c_1v_1`$. By Lemma 6.2, $`c_1(\tau )`$ must satisfy the differential equation

$$\frac{\kappa }{p+1}\left(\frac{d}{d\tau }c_1(\tau )\right)\theta _{2,1}(0,\tau )=2\left(\frac{d}{d\tau }\mathrm{ln}\vartheta _1^{}(0,\tau )\right)c_1(\tau )\theta _{2,1}(0,\tau )+2c_1(\tau )\frac{d}{d\tau }\theta _{2,1}(0,\tau ).$$

The function $`c_1(\tau )=\left(4\pi \theta _{2,1}(0,\tau )\vartheta _1^{}(0,\tau )^1\right)^{2(p+1)/\kappa }`$ is a solution of the above equation. By Corollary 7.2, $`c_1(\tau )=\left(\varphi _3(\tau )\eta (\tau )^1\right)^{4(p+1)/\kappa }`$. The coefficient of proportionality is computed in the limit $`\tau i\mathrm{}`$, cf. the proof of (2). This proves (5). To prove (6), we apply the transformation $`S`$ to both sides of (5). To prove (7), we apply the transformation $`T`$ to both sides of (6). This completes the proof.

### 8.4. Proof of (8), (9), and (10)

Let $`\kappa =2p+6`$. Let $`v_j(\lambda ,\tau )=\vartheta _1^{p+1}(\lambda ,\tau )\theta _{4,j}^s(\lambda ,\tau ),\mathrm{\hspace{0.17em}0}j4`$. Any solution of (1) with the properties (i)-(iv) has the form $`v(\lambda ,\tau )=_{j=0}^4c_j(\tau )v_j(\lambda ,\tau )`$. Let $`A`$ and $`B`$ be the transformations in section $`4`$. By Proposition 4.1, $`Av=_{j=0}^4(1)^{p+j+1}c_jv_j`$ is also a solution. Hence the function $`c_0v_0+c_2v_2+c_4v_4`$ gives a solution too. Moreover, $`B(c_0v_0+c_2v_2+c_4v_4)=(1)^{p+1}(c_4v_0+c_2v_2+c_0v_4)`$ is also a solution. So there exists a solution of the form $`c(\tau )(v_0v_4)`$. It is an eigenvector of $`A`$ with eigenvalue $`(1)^{p+1}`$ and an eigenvector of B with eigenvalue $`(1)^p`$. We show that the subspace of conformal blocks with this property is one-dimensional. By Theorem 3.2, the space of conformal blocks is five-dimensional with spanning set $`\{u_{\kappa ,n}|p+1np+5\}`$. By Lemma 4.3, the eigenspace of $`A`$ corresponding to the eigenvalue $`(1)^{p+1}`$ is three-dimensional and is spanned by $`u_{\kappa ,p+1}`$, $`u_{\kappa ,p+3}`$, and $`u_{\kappa ,p+5}`$. By Lemma 4.3, the transformation $`B`$ preserves the subspace $`u_{\kappa ,p+1},u_{\kappa ,p+3},u_{\kappa ,p+5}`$. The matrix of $`B`$ restricted to this subspace is

$$\left(\begin{array}{ccc}0& 0& e^{2\pi i\frac{p(p+5)}{\kappa }}\\ 0& (1)^{p+1}& 0\\ e^{2\pi i\frac{p(p+1)}{\kappa }}& 0& 0\end{array}\right).$$

Thus the restriction of $`B`$ to $`u_{\kappa ,p+1},u_{\kappa ,p+3},u_{\kappa ,p+5}`$ has eigenvalues $`(1)^p`$ and $`(1)^{p+1}`$ of multiplicities $`1`$ and $`2`$, respectively. The eigenspace corresponding to the eigenvalue $`(1)^p`$ is spanned by the vector $`u_{\kappa ,p+1}+(1)^{p+1}e^{2\pi ip(p+1)/\kappa }u_{\kappa ,p+5}.`$ It follows that this vector is proportional to $`c(\tau )(v_0v_4)`$. By Lemma 6.2, $`c(\tau )`$ must satisfy the differential equation

$$\begin{array}{c}\frac{\kappa }{p+1}\left(\frac{d}{d\tau }c(\tau )\right)(\theta _{4,0}(0,\tau )\theta _{4,4}(0,\tau ))=\hfill \\ \hfill 4\left(\frac{d}{d\tau }\mathrm{ln}\vartheta _1^{}(0,\tau )\right)c(\tau )(\theta _{4,0}(0,\tau )\theta _{4,4}(0,\tau ))+6c(\tau )\frac{d}{d\tau }(\theta _{4,0}(0,\tau )\theta _{4,4}(0,\tau )).\end{array}$$

The function $`c(\tau )=\left((2\pi )^2(\theta _{4,0}(0,\tau )\theta _{4,4}(0,\tau ))^3\vartheta _1^{}(0,\tau )^2\right)^{2(p+1)/\kappa }`$ is a solution of the above equation. By Lemma 7.3, we have $`c(\tau )=2^{3(p+1)/\kappa }(\varphi _3(\tau )\eta (\tau ))^{6(p+1)/\kappa }`$. The coefficient of proportionality is computed as in the proof of (2). This proves (8). To prove (9) and (10), apply the transformations $`S`$ and $`TS`$, respectively, to both sides of (8). This completes the proof.

### 8.5. Proof of (11)

Let $`\kappa =2p+8`$. Let $`v_j(\lambda ,\tau )=\vartheta _1^{p+1}(\lambda ,\tau )\theta _{6,j}^s(\lambda ,\tau ),\mathrm{\hspace{0.17em}0}j6.`$ Any solution of (1) with the properties (i)-(iv) has the form $`v(\lambda ,\tau )=_{j=0}^6c_j(\tau )v_j(\lambda ,\tau )`$. Let $`A`$ and $`B`$ be as in section $`4`$. By Proposition 4.1, $`Av=_{j=0}^6(1)^{p+j+1}c_jv_j`$ is also a solution. Hence the function $`c_1v_1+c_3v_3+c_5v_5`$ gives a solution too. The function $`B(c_1v_1+c_3v_3+c_5v_5)=(1)^{p+1}(c_5v_1+c_3v_3+c_1v_5)`$ also gives a solution. So there is a solution of the form $`c(\tau )(v_1v_5)`$ which is an eigenvector of $`A`$ and $`B`$ with eigenvalue $`(1)^p`$ under both transformations. We show that the subspace of conformal blocks with this property is one-dimensional. By Theorem 3.2, the space of conformal blocks is seven-dimensional with spanning set $`\{u_{\kappa ,n}|p+1np+7\}`$. By Lemma 4.3, the three-dimensional eigenspace of $`A`$ corresponding to the eigenvalue $`(1)^p`$ is spanned by $`u_{\kappa ,p+2}`$, $`u_{\kappa ,p+4}`$, and $`u_{\kappa ,p+6}`$. The matrix of $`B`$ restricted to $`u_{\kappa ,p+2},u_{\kappa ,p+4},u_{\kappa ,p+6}`$ is

$$\left(\begin{array}{ccc}0& 0& e^{2\pi i\frac{p(p+6)}{\kappa }}\\ 0& (1)^{p+1}& 0\\ e^{2\pi i\frac{p(p+2)}{\kappa }}& 0& 0\end{array}\right).$$

Thus, $`B`$ has eigenvalues $`(1)^p`$ and $`(1)^{p+1}`$ of multiplicities $`1`$ and $`2`$, respectively. The eigenspace corresponding to the eigenvalue $`(1)^p`$ is spanned by the vector $`u_{\kappa ,p+2}+(1)^{p+1}e^{2\pi ip(p+2)/\kappa }u_{\kappa ,p+6}.`$ So this vector is proportional to $`c(\tau )(v_1v_5)`$. By Lemma 6.2, $`c(\tau )`$ must be a solution of the differential equation

$$\begin{array}{c}\frac{\kappa }{p+1}\left(\frac{d}{d\tau }c(\tau )\right)(\theta _{6,1}(0,\tau )\theta _{6,5}(0,\tau ))=\hfill \\ \hfill 6\left(\frac{d}{d\tau }\mathrm{ln}\vartheta _1^{}(0,\tau )\right)c(\tau )(\theta _{6,1}(0,\tau )\theta _{6,5}(0,\tau ))+10c(\tau )\frac{d}{d\tau }(\theta _{6,1}(0,\tau )\theta _{6,5}(0,\tau )).\end{array}$$

The function $`c(\tau )=\left((2\pi )^3(\theta _{6,1}(0,\tau )\theta _{6,5}(0,\tau ))^5\vartheta _1^{}(0,\tau )^3\right)^{2(p+1)/\kappa }`$ is a solution of the preceding equation. By Lemma 7.5, $`c(\tau )=(\eta (\tau ))^{8(p+1)/\kappa }`$. The coefficient of proportionality is computed as in the proof of (2). Notice that the solution in (11) is invariant with respect to the action of the modular group.∎