text:
# Numerical Analysis of the Double Scaling Limit in the bosonic part of the IIB Matrix Model presented by S.Horata

## 1 Introduction

Recently, several models have been proposed as non-perturbative formulations of string theory. The IIB matrix model, which is zero-volume limit of the ten-dimensional large $`N`$ supersymmetric Yang-Mills theory, has been proposed as a constructive definition of the type IIB superstring theory. It is defined by the following action,

$$S=\frac{1}{g^2}\mathrm{tr}\left(\frac{1}{4}[A_\mu ,A_\nu ]^2+\frac{1}{2}\overline{\psi }\mathrm{\Gamma }^\mu [A_\mu ,\psi ]\right),$$
(1)
where $`A_\mu `$ and $`\psi `$ are $`N\times N`$ traceless Hermitian matrices. The interesting feature of this model is that the space-time coordinates are considered as the eigenvalues of these matrices. Then, we expect that the fundamental issues including the dimensionality and the quantum gravity can be understood by studying the dynamics of the model.

To take the continuum limit ($`N\mathrm{}`$) for the IIB Matrix model, a sensible double scaling limit should be determined dynamically for the scaling property of the two important quantities of the model, the string scale ($`\alpha ^{}`$) and the string coupling constant ($`g_{str}`$).

$`\alpha ^{}`$ $``$ $`g^2N^{a+b}=g^2N^\gamma Constant,`$
$`g_{str}`$ $``$ $`N^{ab},`$ (2)
where one have a restriction, $`a=b`$ for the finite value of the string coupling constant ($`g_{str}`$).

For studying the dynamical aspects of the model, some Monte Carlo simulations have been performed. In Ref., the bosonic model has been studied with a analytical method as well as a numerical one. The numerical results support the $`1/D`$ expansion as an effective tool detecting the large $`N`$ scaling behavior of the model. The leading term of the $`1/D`$ expansion at $`D>3`$ suggests that the exponent ($`\gamma `$) takes a value of 1.

In analogy with the two-dimensional model, the Eguchi-Kawai model, we study the scaling property of the Wilson loop in ten-dimension. The area law of the Wilson loop operator has been found in the four-dimensional model in Ref.. Then we consider that the scaling property of the string tension is estimated in ten-dimensional model with investigating the scaling property of the Wilson loop in ten-dimension.

In this article, we present the data which show the existence of the double scaling limit in the model and discuss our numerical results.

## 2 Monte Carlo simulation and scaling behavior

For numerical simulation of the model, we consider a bosonic version of the IIB matrix model, namely a model which can be obtained from the IIB matrix model by dropping the fermions by hand. In the IIB matrix model, supersymmerty is expected to be very important as is discussed in many places of the literature. However, it is difficult to calculate with numerical method. In ref., numerical study of a supersymmetric version of the IIB matrix model is discussed. In this article, we study the most simplest model. The bosonic model of the IIB matrix model is given by

$$S_{bosonic}=\frac{1}{g^2}\mathrm{tr}[A_\mu ,A_\nu ][A^\mu ,A^\nu ],$$
(3)
where $`A_\mu `$ are $`N\times N`$ Hermitian matrices representing the ten-dimensional gauge fields. The coupling constant ($`g`$) is nothing but a scale parameter and is absorbed with the rescaling of the gauge field ($`A_\mu `$) as $`A_\mu \frac{1}{\sqrt{g}}A_\mu `$.

we consider the partition function of bosonic model in ten-dimension,

$$Z=𝑑Ae^{\frac{1}{g^2}\mathrm{tr}[A_\mu ,A_\nu ][A^\mu ,A^\nu ]},$$
(4)
where the measure of gauge fields is defined by

$$dA=\underset{\mu =1}{\overset{10}{}}\left[\underset{i=1}{\overset{N}{}}\underset{j=i}{\overset{N}{}}d(A_\mu ^{ij})\right].$$
(5)
The action is quadratic with respect to each component, which means that we can update each component by generating gaussian random number in the heat-bath and the Metropolis algorithm.

To investigate the scaling behavior of the string scale ($`\alpha ^{}`$), we study the Wilson loop operator in the IIB matrix model.

The Wilson loop operator and the large $`N`$ behavior have been studied with the light-cone string field theory of the type IIB superstring. The Wilson loop operator ($`w(C)`$) is defined as,

$$w(C)=\mathrm{tr}(v(C)),$$
(6)
where $`C`$ denotes the closed path and $`v(C)`$ is defined as $`v(C)=U_\mu \mathrm{}U_\mu =P_C\mathrm{exp}(i𝑑\sigma k_\mu A_\mu )`$ in the bosonic model. The matrices ($`U_\mu `$) are considered as the unitary matrices,

$$U_\mu =\mathrm{exp}(i𝑑lA_\mu ).$$
(7)
In the ordinary lattice gauge field theory, the expectation value of the Wilson loop operator which spreads a large area behaves as follows,

$$w(I,J)\mathrm{exp}(KI\times J),$$
(8)
where $`I`$ and $`J`$ are the side lengths of the rectangular loop and $`K`$ denotes the string tension. In the same analogy, we study the Wilson loop operator of the bosonic model of the IIB matrix model.

From the scaling relation of two-dimensional Eguchi-Kawai model,

$$g^2NConstant,$$
(9)
It is expected that the similar scaling relation also holds in the IIB matrix model as

$$\alpha ^{}g^2N^\gamma Constant.$$
(10)
The exponent ($`\gamma `$) should be determined dynamically from the model. For the large $`N`$ limit, the parameter ($`g^2N^\gamma `$) must be fixed in the IIB matrix model in the same manners as the parameter ($`g^2N`$) must be fixed in the Eguchi-Kawai model.

It is considered that the bosonic model is equivalent to the $`D>2`$ Eguchi-Kawai model in the weak coupling limit. Since the $`U(1)^D`$ symmetry rotates all the eigenvalues by the same angle, the following expansion is valid in the weak coupling region,

$$U_\mu e^{i\alpha _\mu }e^{iA_\mu },$$
(11)
where $`\alpha _\mu `$ take constant values due to the $`U(1)^D`$ symmetry and $`A_\mu `$ are small. The bosonic model action can be obtained by expanding the action of the ten-dimensional Eguchi-Kawai model in terms of $`A_\mu `$. When the higer order terms of $`A_\mu `$ can be neglected, we can obtain the area law of the Wilson loop operator as

$`w(I\times J)`$ $`=`$ $`<\mathrm{tr}(U_\mu \mathrm{}U_\mu )_{I\times J}>`$
$``$ $`<\mathrm{tr}(e^{iA_\mu }\mathrm{}e^{iA_\mu })_{I\times J}>`$
$``$ $`\mathrm{exp}(K(I\times J)).`$
In Ref., the area dependence of the Wilson loop operator has been measured and found the area law in four-dimension, $`D=4`$.

We calculate the eigenvalue of Wilson loop operator numerically in ten-dimension, $`D=10`$,

$$w(C)=<\frac{1}{N}\mathrm{tr}(e^{iA_\mu }\mathrm{}e^{iA_\nu })>.$$
(13)
We take the loop ($`C`$) as the rectangular ($`I\times J`$) where we select any two direction ($`\mu ,\nu `$) in ten-dimension. We consider that in the bosonic model the isotropy of the ten-dimensional space-time is not broken down spontaneously, since the dependence of the direction of the rectangular is not found.

From the numerical results, the Wilson loop operator closes to exponential curve with the large size area as

$$w(I,J)\mathrm{exp}(KI\times J),$$
(14)
where $`K`$ denotes string tension. We plot the results in Fig.1.

Then we can estimate the string scale ($`\alpha ^{}`$) as

$$\alpha ^21/K^2g^2N^{1.07(1)}=Constant.$$

## 3 Summary

Let us summarize the main points made in our calculation. We confirm that the Wilson loop operator in the ten-dimensional bosonic model obeys the area law. For the scaling behavior of the bosonic model of the IIB matrix model, our numerical estimation is

$$\alpha ^2g^2N^{1.07(1)}=Constant.$$
(15)
Furthermore the numerical result suggests that the large $`N`$ behavior of the square root of the string tension approximately equal to the inverse of the extent of the space-time.

$$K^{1/2}g^1N^{0.54(1)}R^1.$$
(16)
It means that the Planck scale of the theory has the same scaling property of the extent of the space time in ten-dimension. The numerical result is consistent to the suggestion from the string field theory on light-cone frame and the $`1/D`$ expansion.


format: text

text:
# The Structure of Self-Gravitating Polytropic Systems with 𝑛 around 5.

## 1. Introduction

Galaxies, star clusters, and other stellar systems are stationary assemblages of stars bound by a self-consistent, mean gravitational potential. Their structure is determined by the nonlinear Poisson equation, which, in most cases, may be solved only numerically. Therefore, any analytical study acquires special interest.

The simplest model of a gravitationally bound stellar system is a spherically symmetric polytrope with a power-law distribution function of stars. The gravitational potential of a polytropic sphere is determined by the Lane–Emden equation:

$$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}\psi (r)\right)=\psi ^n,$$
(1)
where $`r`$ is the dimensionless spherical radius, $`\psi (r)`$ is the dimensionless potential, and $`n`$ is the polytropic index (see Binney & Tremaine, 1987 for more discussion). The density is found from $`\rho =c_n\psi ^n`$, where $`c_n`$ is a constant which depends on $`n`$.

The stellar (collisionless) polytropes are related to hydrostatic equilibria of self-gravitating spheres of a polytropic gas with the equation of state $`P=K\rho ^\gamma `$, where $`\gamma =1+\frac{1}{n}`$, $`P`$ is the gas pressure, and $`K`$ is a constant. The gas polytropes with $`\gamma <6/5`$ (i.e., $`n>5`$) are usually used to model the structure of molecular clouds which has a major effect on the star formation that occurs within them. Note that the $`\gamma 1`$ ($`n\mathrm{}`$) limit corresponds to the case of an isothermal gas. The gravitationally stable polytropes111 Locally adiabatic polytropes with polytropic indexes $`\gamma >6/5`$ are stable if their adiabatic indexes $`\gamma _{\mathrm{ad}}>4/3`$; globally adiabatic polytropes are unconditionally stable for $`\gamma =\gamma _{\mathrm{ad}}>4/3`$, see discussion below and, e.g., McKee & Holliman II (1999). with $`\gamma >6/5`$ (i.e., $`n<5`$) often serve as models of stars, because they can have a zero pressure boundary at a finite radius — the star surface. Interestingly, the Lane-Emden equation appears also as an evolution equation for a scalar field in some versions of the quintessence model (Liddle & Scherrer, 1999). Hence the present study is also of practical interest.

There are two special cases for which nonsingular analytical solutions of equation (1) are known, namely $`n=1`$ and $`n=5`$. We performed a study of (1) using a powerful mathematical technique — the Lie group analysis of differential equations (Ibragimov, 1985). This cumbersome analysis, the details of which are omitted from this publication, yields that, for any $`n`$ except $`n=1`$ and $`5`$, equation (1) is non-integrable in a closed form.222 Strictly speaking, equation (1) has a zero-dimensional Lie algebra, meaning that there is no a point transformations group admitted by (1). This does not prove, however, the global non-integrability. A remarkable counter-example is a class of solitonic nonlinear differential equations which “became” integrable after the relation of these equations to the Bäcklund transformations group have been discovered. In most cases, however, the absence of a Lie algebra often means nonintegrability of a differential equation. When $`n=1`$, equation (1) becomes the linear Helmholtz equation, and when $`n=5`$, the solution has been discovered by Schuster (1883) and takes a simple form:

$$\psi (r)=\left(1+\frac{1}{3}r^2\right)^{1/2}.$$
(2)
For other values of $`n`$, the general solution may be found only numerically.

The present study is motivated by a very peculiar dependence of the radial profiles of the gravitational potential and density on the value of $`n`$ when $`n5`$. Several solutions of equation (1) with the “core” boundary conditions (i.e., $`\psi (0)=1,_r\psi (0)=0`$, where $`_r\frac{}{r}`$) are represented in Figure The Structure of Self-Gravitating Polytropic Systems with $`n`$ around 5.. A small change in $`n`$ results in a drastic modification of the profiles. Solutions with $`n>5`$ exhibit a multiple core structure, — the fact which has probably been known but never acknowledged in the literature. In this solutions, the regions where $`\psi `$ falls as a power-law are interchanged with flat cores where $`\psi \text{constant}`$. For $`n<5`$, the solutions are spatially bound, i.e., $`\psi `$ and, hence, $`\rho `$ vanish at a finite radius. The radius at which the second core begins to form in the $`n>5`$ solutions is comparable to the radius at which the density vanishes in the $`n<5`$ solutions.

In this paper we develop a functional series expansion technique which allows the derivation of solutions with the polytropic index $`n\pm |ϵ|`$, given the solution with the index $`n`$ (for any $`n`$ and $`ϵ,|ϵ|1`$), to any accuracy in powers of $`ϵ`$ and present it in §2. We apply this technique to the $`n=5`$ solution (2) and analytically derive a first-order ($`ϵ^1`$) approximate solution in §3. This solution perfectly agrees with the exact numerical $`n<5`$ solutions everywhere, and with the $`n>5`$ solutions at radii smaller than the second core radius. We also provide a useful approximation for $`\psi (r)`$ for $`n5`$ to any power in $`ϵ`$. These results are directly applicable to the polytropic gas systems with $`\gamma 6/5`$. In §4 we discuss stability of stellar polytropes. We discuss the results and make some observational predictions in §5.

## 2. Analytical study of polytropic solutions: formalism

Let us rewrite equation (1) in a form of the nonlinear Poisson equation with some boundary conditions,

$`^2\psi _n(𝐱)+\psi _n^n(𝐱)=0,`$ (3)
$`\psi _n(𝐱)|_S=f(𝐱),`$
where $`^2`$ is the three-dimensional Laplace operator, $`f(𝐱)`$ is some specified function, and $`S`$ is a two-dimensional surface. The subscript $`n`$ explicitly shows that $`\psi _n(𝐱)`$ is a function of the parameter $`n`$, as well as three-dimensional coordinates $`𝐱`$. Our goal is to find a solution $`\psi _{n+ϵ}(𝐱)`$, given the solution $`\psi _n(𝐱)`$ for positive or negative $`ϵ`$. For $`|ϵ|1`$, such a solution may be expressed as a functional expansion in powers of $`ϵ`$:

$$\psi _{n+ϵ}(𝐱)=\psi _n(𝐱)+ϵ\frac{\psi _i(𝐱)}{i}|_{i=n}+\frac{ϵ^2}{2!}\frac{^2\psi _i(𝐱)}{i^2}|_{i=n}+\mathrm{}.$$
(4)
Note, the derivative terms above are functions of coordinates $`𝐱`$ evaluated at a particular value of the polytropic index. Hereafter we use the short-hand notation for the derivatives, $`\frac{^m\psi _i(𝐱)}{i^m}|_{i=n}_n^m\psi _n(𝐱),m=1,2,\mathrm{}`$. The derivatives $`_n^m\psi _n`$ are found by solving equations (5) and (6) below. Differentiating (3) with respect to $`n`$, we obtain

$`^2(_n\psi _n)+n\psi _n^{n1}(_n\psi _n)=\psi _n^n\mathrm{ln}\psi _n,`$ (5)
$`_n\psi _n(𝐱)|_S=0,`$
where the last equation means that boundary conditions are independent of $`n`$. Equation (5) is linear in the function $`_n\psi _n`$ and, hence, is easier to solve than the original (3). Once the solution of a homogeneous equation \[i.e., equation (5) with vanishing right hand side\] is found, it is straightforward to find a solution of inhomogeneous equation (5), as shown in the Appendix. The quantity $`_n\psi _n`$ represents the first-order correction in equation (4). Higher-order terms are obtained by taking higher derivatives of (3) and solving for appropriate functions. For instance, upon differentiating (3) $`m`$ times, the equation for the $`m`$-th derivative, $`_n^m\psi _n`$, takes the form

$`^2(_n^m\psi _n)+n\psi _n^{n1}(_n^m\psi _n)`$
$`=H_{n,m}(\psi _n,_n\psi _n,\mathrm{},_n^{m1}\psi _n),`$ (6)
$`_n^m\psi _n(𝐱)|_S=0.`$
Here $`H_{n,m}(\mathrm{})`$ is a function which depends on $`n`$, $`m`$, $`\psi _n`$, and all previously found functions $`_n^k\psi _n,1km1`$. Note that the homogeneous solutions of equation (6) for all $`m`$ are the same and identical to the homogeneous solution of equation (5). Hence, the calculation of $`_n^m\psi _n`$ is staightforward (see Appendix).

## 3. The structure of stellar polytropes with $`n=5\pm |ϵ|`$

Equations (4), (6) allow the solutions $`\psi _{n+ϵ}`$ to be calculated from a solution $`\psi _n`$. As mentioned above, there are two cases in which $`\psi _n`$ is available. For $`n=1`$, equation (3) is linear and the complete three-dimensional solution, $`\psi _1(𝐱)`$, is known. For $`n=5`$, only the spherically symmetric solution, $`\psi _5(r)`$, has been discovered, which is given by equation (2). Here we consider the latter case and find an analytical expression for the first-order correction.

The function $`\psi _5(r)`$ is given by equation (2). Equation (5) for $`n=5`$ takes the form,

$$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}\varphi (r)\right)+\frac{5}{(1+r^2/3)^2}\varphi (r)=\frac{1}{2}\frac{\mathrm{ln}(1+r^2/3)}{(1+r^2/3)^{5/2}},$$
(7)
where we again introduced a short-hand notation $`\varphi (r)=_n\psi _5=\frac{\psi _i(r)}{i}|_{i=5}`$. The general solution of this equation may be obtained if a special solution of the homogeneous equation \[i.e., equation (7) with the zero right hand side\] is known, as shown in the Appendix. Since $`\psi _5`$ is a power-law of a rational function, we guess that a special solution of the homogeneous equation is a function of the same kind. We have the first homogeneous solution,

$$\varphi _1(r)=\frac{r^23}{(r^2+3)^{3/2}}$$
(8)
The second homogeneous solution is obtained from equation (A4) where $`\mathrm{\Delta }=1/r^2`$, as follows from equation (A3). It reads

$$\varphi _2(r)=\frac{r^418r^2+9}{r(r^2+3)^{3/2}}.$$
(9)
The general solution is found from equation (A2) by straightforward integration. The constants of integration are found from boundary conditions. First, we require the solution to be non-diverging everywhere; hence $`C_2=0`$. Second, we consider the solutions with a core, i.e., $`\psi _5(0)=\psi _{5+ϵ}(0)=1`$; hence, $`\varphi (0)=_n\psi _5(0)=0`$. The last condition yields $`C_1=\sqrt{3}/12`$. Finally, we have

$$\psi _{5+ϵ}(r)\psi _5(r)+ϵ_n\psi _5(r),$$
(10)
where

$`\psi _5(r)`$ $`=`$ $`\left(1+{\displaystyle \frac{r^2}{3}}\right)^{1/2},`$ (11)
$`_n\psi _5(r)`$ $`=`$ $`{\displaystyle \frac{(r^418r^2+9)}{16r(r^2+3)^{3/2}}}\text{arc\hspace{0.17em}tan}\left({\displaystyle \frac{r}{\sqrt{3}}}\right)`$ (12)
$`+{\displaystyle \frac{3^{3/2}/4}{(r^2+3)^{3/2}}}\left[{\displaystyle \frac{7r^2}{36}}{\displaystyle \frac{1}{4}}+\mathrm{ln}\left(1+{\displaystyle \frac{r^2}{3}}\right)\right].`$
The approximate solution (10) is shown in Figure The Structure of Self-Gravitating Polytropic Systems with $`n`$ around 5. by dotted curves for several values of $`ϵ`$. The function $`_n\psi _5`$ is also shown by the dashed line. It is seen that solutions with $`n<5`$ are very well approximated by the first-order expansion (10) at all radii. In contrast, for solutions with $`n>5`$, this approximation is valid for $`r`$ less than or comparable to the second core radius, i.e., where the second “plateau” ends. To obtain the profiles at larger radii, higher-order terms in $`ϵ`$ must be kept. The next, $`ϵ^2`$-term can also be calculated analytically. It is, however, cumbersome and involves generalized hypergeometric functions and, therefore, is not presented here. As for the higher-order corrections, the complexity of the resulting integrals seems to preclude further analytical treatment.

We can now estimate the inner radius of the second core as follows. The asymptotic value of $`_n\psi _5`$ at $`r1`$ is $`_n\psi _5(\mathrm{})=\pi /32`$. The transition from a power-law decay to a flat core occurs when the first and second terms are of comparable magnitude, $`\psi _5(r_{})\pi |ϵ|/32`$. Thus, the transition radius is

$$r_{}\frac{32\sqrt{3}}{\pi |ϵ|}.$$
(13)
Note that for $`ϵ<0`$ this radius determines the size of the system, i.e., the radius at which the density vanishes.

Without proof we suggest that the inner and outer core radii of other cores are also inversely proportional to the powers of $`ϵ`$. Using this conjecture, we present a heuristic expression which approximates the $`n=5+|ϵ|`$ solution quite well. The expression is best written as a continuous fraction,

$$\psi _{5+|ϵ|}(r)\frac{1}{\sqrt{1+{\displaystyle \frac{r^2/r_1^2}{1+ϵ^2{\displaystyle \frac{r^2/r_2^2}{1+ϵ^4{\displaystyle \frac{r^2/r_3^2}{1+\mathrm{}}}}}}}}},$$
(14)
where $`r_1,r_2,r_3,\mathrm{}`$ are constants independent of $`ϵ`$ which are to be found by fitting to the exact numerical solution with $`ϵ0`$. For the first two radii, we use the analytical expressions: $`r_1=\sqrt{3}`$ and $`r_2=r_{}|ϵ|`$. Truncating the fraction at the $`r_7`$ term and performing the nonlinear fit to the exact solution for $`ϵ=.001`$, we have

$`r_1`$ $`=`$ $`\sqrt{3}1.73,`$
$`r_2`$ $`=`$ $`32\sqrt{3}/\pi 17.6,`$
$`r_3`$ $``$ $`1.95\times 10^2,`$
$`r_4`$ $``$ $`9.29\times 10^2,`$ (15)
$`r_5`$ $``$ $`4.75\times 10^3,`$
$`r_6`$ $``$ $`1.60\times 10^4,`$
$`r_7`$ $``$ $`5.00\times 10^5.`$
The polytropes with $`n=5|ϵ|`$ are well described by the following approximate expression

$$\psi _{5|ϵ|}(r)\frac{1}{\sqrt{1+r^2/3}}|ϵ|\frac{\pi }{32},$$
(16)
which is easily obtained from equation (10). No higher order corrections are needed in this case.

How do the $`n>5`$ stellar polytropic spheres look to an observer? Assuming that the luminosity density is proportional to the number density of stars, the surface brightness profile is calculated as follows (Binney & Tremaine (1987)),

$$I(R)=I_0_R^{\mathrm{}}\frac{\rho (r)rdr}{\sqrt{r^2R^2}},$$
(17)
where $`R`$ is the radial distance from the center of an object on the plane of the sky and $`I_0`$ is a normalization. Note, if $`\rho (r)`$ falls too slowly and the integral fails to converge, an upper limit must be replaced with an effective tidal radius to truncate the system. Using approximations (14)–(15) and that $`\rho _n(r)\psi _n^n(r)`$, we numerically determine the surface brightness profiles, as presented in Figure The Structure of Self-Gravitating Polytropic Systems with $`n`$ around 5.. As one can see, photometric observations will reveal a bright core surrounded by a much larger and considerably fainter shell or “halo” throughout which the surface brightness is almost constant. Other halos which are located at larger radii are perhaps too faint to be observable.

We can also estimate the profile of the velocity dispersion of stars $`\sigma ^2=v^2v^2`$ which is analogous to the sound speed in a gas, $`\sigma ^2(r)P/\rho \psi _n(r)`$. Therefore, the velocity dispersion decreases with radius and exhibits the same multi-core structure. In particular, the faint halo is cooler than the central core and looks nearly isothermal because $`\sigma ^2\text{ constant}`$ thoughout.

## 4. Stability of the Polytropes

The stability of self-gravitating polytropes depends on how the system responses to an adiabatic perturbation, i.e., whether the internal heat transfer through the system is long or short compared to the dynamical scale (McKee & Holliman II, 1999). In the first case — a locally adiabatic polytrope — the polytropic index, which determines the global structure, may differ from the adiabatic index $`\gamma _{\mathrm{ad}}`$, which determines thermodynamical properties. The second case represents a globally adiabatic polytrope in which $`\gamma _{\mathrm{ad}}=\gamma `$.

The stability analysis of nonsingular polytropic spheres usually requires numerical calculations, even for the simplest cases (e.g., $`n=5`$). Here is a brief summary of the results (see McKee & Holliman II, 1999 and references therein). First, polytropes with $`\gamma >4/3`$ (i.e., $`n<3`$) are unconditionally stable for any value of $`\gamma _{\mathrm{ad}}`$. Second, locally adiabatic systems with $`6/5<\gamma 4/3`$ are stable if $`\gamma _{\mathrm{ad}}>4/3`$. A system with smaller $`\gamma _{\mathrm{ad}}`$, e.g., a globally adiabatic polytrope, is stable only if it is truncated at a certain radius and confined by an appropriate external pressure. Third, for the systems with $`\gamma 6/5`$, the truncation radius depends on the value of $`\gamma _{\mathrm{ad}}`$. This radius approaches infinity for

$$\gamma _{\mathrm{ad},\mathrm{}}=\frac{32\gamma (2\gamma )}{(6\gamma )^2}<\gamma .$$
(18)
As one can see, globally adiabatic polytropes with $`\gamma <4/3`$ are stable only if they are pressure confined (in an analogy with Bonnor-Ebert isothermal sphere). The truncation radius follows from the solution of equation (54) of McKee & Holliman II (1999)333 There is a typo: in the first equality the power of $`c_s`$ is 3, not 3/2; cf., their equation (10). which in our dimensionless units reads as

$$m_{\mathrm{cr}}_0^{r_{\mathrm{cr}}}4\pi r^2\rho (r)𝑑r=\left(\frac{4\pi }{n+1}\right)^{\frac{3}{2}}\mu _{\mathrm{cr}}\left[\rho (r_{\mathrm{cr}})\right]^{\frac{3}{2n}\frac{1}{2}}.$$
(19)
Here $`\mu _{\mathrm{cr}}`$ is a factor determined numerically; it is close to unity for large $`n`$ and monotonically approaches $`4.6`$ as $`n3`$. For $`n=5`$ and $`r_{\mathrm{cr}}r_1`$ this equation yields

$$r_{\mathrm{cr}}=\frac{18}{\mu _{\mathrm{cr}}}\sqrt{\frac{3}{2\pi }}.$$
(20)
Comparing this to equation (13), we see that $`r_{\mathrm{cr}}r_{}`$ corresponds to

$$ϵ\frac{16}{9}\sqrt{\frac{2}{\pi }}\mu _{\mathrm{cr}}1,$$
(21)
i.e., $`n6`$. Therefore, the multiple core structure in a globally adiabatic system may be seen if its polytropic index is $`n6`$. The stability analysis of these systems requires the use of the density profiles given by equations (10) or (14). The accuracy of our approximate solution degrades for $`ϵ>1`$. The qualitative behavior, however, remains correct (note that in the limit of an isothermal sphere $`ϵ\mathrm{}`$, some oscillations in the profile persist). On the other hand, the multiple core structure may be seen in pressure confined, locally adiabatic ($`\gamma _{\mathrm{ad}}\gamma `$) polytropic systems with $`n>5`$ and $`\gamma _{\mathrm{ad}}`$ being sufficiently close to unity, so that $`r_{\mathrm{cr}}r_{}`$.

## 5. Discussion

In this paper we have demonstrated that the radial structure of polytropic self-gravitating spheres is very sensitive to the value of the polytropic index $`n`$ when $`n5`$ and have investigated it analytically. We developed a method which allows one to obtain approximate solutions, $`\psi _{n+ϵ}`$, of the Lane-Emden equation, provided the solution $`\psi _n`$ is known and $`ϵ`$ is a positive or negative constant, $`|ϵ|1`$. We applied this method to $`n=5`$ solution and obtained the first-order (in $`ϵ`$) correction analytically. The obtained solution describes finite-size systems with $`n<5`$ and the formation of the second core in systems with $`n>5`$. Both effects occur at the radius $`r_{}32\sqrt{3}/\pi |ϵ|`$. We also provided a continuous fraction approximation for the polytropes with $`n>5`$, equation (14), and a simple approximate formula for those with $`n<5`$, equation (16).

We make the observational predictions by calculating numerically the surface brightness of a stellar polytropic system. A galaxy or a star cluster which is described by a polytrope with $`n5`$ will be observed as a bright core surrounded by a very extended faint halo or shell. This halo is cooler than the central core and looks nearly isothermal. Because of the large spatial size and low surface brightness this halo may easily be confused with a background. A way to disentangle the two is to probe a correlation between the core radii, $`r_1,r_2`$, and (if detected) $`r_3`$ and the surface brightness distribution.

Due to the duality between the stellar collisionless polytropes and the gas polytropes, our results are immediately applicable to the latter with the adiabatic index $`\gamma =1+\frac{1}{n}\frac{6}{5}`$. Such polytropic gas spheres are used in studies of molecular clouds and star formation. The solutions with $`\gamma <6/5`$ describe extended molecular clouds bounded by external pressure. Equation (14) describes the multi-core structure of a cloud with $`\gamma 6/5`$. The density and, hence, pressure is almost constant in core/halo regions. Therefore, if an external pressure, $`p_{ext}`$, is comparable to the pressure in one of the halos, then even a small variation in $`p_{ext}`$ will drastically change the radius of the cloud. Thus, a large dispersion of masses and sizes of the clouds embedded into nearly the same environments may indicate that the gas in the cloud has $`\gamma 6/5`$, although other explanations are also possible. The gravitationally stable polytropes with $`\gamma >6/5`$ are often used as models of stars because the pressure vanishes at a finite radius, — the star’s surface. If $`\gamma 6/5`$, the radial profile is well described by equation (16) and the radius of a star is given by equation (13).

We thank Ramesh Narayan for interesting discussions and an anonymous referee for very careful reading of the manuscript and many helpful suggestions. This work was supported in part by NASA grant NAG 5-2837 and NSF grant AST 9820686.

## Appendix A General solution of equation (7)

Let us consider the following differential equation,

$$y^{\prime \prime }+f(x)y^{}+g(x)y=h(x),$$
(A1)
where prime denotes the $`x`$-derivative and $`f(x),g(x)`$, and $`h(x)`$ are some functions. The general solution of this inhomogeneous equation reads,

$$y=C_1y_1+C_2y_2+y_2\frac{h}{\mathrm{\Delta }}y_1𝑑xy_1\frac{h}{\mathrm{\Delta }}y_2𝑑x,$$
(A2)
where $`C_1`$ and $`C_2`$ are constants of integration, $`y_1=y_1(x)`$ and $`y_2=y_2(x)`$ are two linearly independent solutions of the homogeneous equation $`y^{\prime \prime }+f(x)y^{}+g(x)y=0`$. Here also

$$\mathrm{\Delta }=e^{{\scriptscriptstyle f(x)𝑑x}}=y_1y_2^{}y_2y_1^{}.$$
(A3)
If one nontrivial solution of the homogeneous equation, $`y_1`$, is known, then the second may be found from the equation

$$y_2=y_1\frac{\mathrm{\Delta }}{y_1^2}𝑑x.$$
(A4)
Figs. The Structure of Self-Gravitating Polytropic Systems with $`n`$ around 5., The Structure of Self-Gravitating Polytropic Systems with $`n`$ around 5.


text:
# Detection of an X-ray periodicity in the Narrow-line Seyfert 1 Galaxy Mrk 766 with XMM-Newton Based on observations with XMM-Newton, an ESA Science Mission with instruments and contributions directly funded by ESA Member States and the USA (NASA)

## 1 Introduction

Mrk 766 is a nearby ($`z`$ = 0.013) and X-ray bright ($`10^{11}\mathrm{erg}\mathrm{s}^1`$) Narrow-Line Seyfert 1 galaxy (Walter & Fink 1993, Boller et al. 1996, Leighly et al. 1996, Leighly 1999). Due to the extreme spectral and X-ray properties found in Narrow-Line Seyfert 1 galaxies and due to its extreme X-ray brightness, Mrk 766 was proposed for the PV phase for observations with the X-ray satellite XMM-Newton. The high throughput of XMM-Newton compared to previous X-ray missions allows precise studies of the spectral properties, e.g. the shape of the soft excess and of the transition region between the soft excess and the power-law component. In addition, RGS observations allow to perform X-ray spectroscopy of accretion disk lines (e.g. Ross et al. 1999) as well as emission and absorption features originating from the broad- and narrow-line region. The XMM-Newton spectral properties of Mrk 766 are discussed in detail by Page et al. (2001) and Pounds et al. (2001). The long-period orbit of XMM allows to measure the variability power spectra of ultrasoft Narrow-Line Seyfert 1 galaxies far better than has yet been possible before and allow to constrain nonlinear variability, to search for spectral variability and to search for any quasi-periodic oscillations (e.g. Papadakis & Lawrence 1995 and references therein). In this paper we concentrate on the timing properties of Mrk 766 and present the discovery of an X-ray periodicity in Mrk 766. The object resembles IRAS 18325$``$5926 (Iwasawa et al. 1998), where a 58000 second periodicity in the 0.5$``$10 keV ASCA energy band was found.

## 2 The XMM-Newton PV observations on Mrk 766 and data reduction

The Narrow-line Seyfert 1 galaxy Mrk 766 was observed during the performance and calibration phase in revolution 0082 at May 20, 2000. The EPIC pn camera was operated in the small window mode implying a frame read out time of 30 msec so that pile-up is not significant. The full frame mode was used for the EPIC MOS cameras and the RGS cameras were operated in the standard spectroscopy+Q mode. In total 13 exposures were taken with the OM using the UVW1, UVW2 and UVM2 filter. At the beginning of the observations the EPIC cameras were blocked for 20 ksec due to high particle background. At the end of the observations all instruments suffered from data gaps. The exposure time used for the spectral and timing analysis presented in this paper is 29000 seconds. For the data reduction and data analysis the standard SAS software packages were used. The spectral fitting and light curve analysis was performed using XSPEC version 10.00 and XRONOS version 4.0, respectively.

### 2.1 Detection of a periodicity in the X-ray light curve of Mrk 766

#### 2.1.1 EPIC pn data analysis

The EPIC pn 0.2–2 keV and 5–10 keV light curves of Mrk 766 are shown in Figure 1. The most striking feature is the presence of a periodic signal in the 0.2–2 keV energy band superimposed on a longer-term count rate increase. Starting at 5000 seconds, the most significant intensity peaks are separated by 4050, 4210, 4210, 3500, 4050 and 4380 seconds. In addition, lower amplitude variations separated by about 2100 seconds are visible (cf. the peak emission at about 27500 and 31000 seconds). The 5–10 keV light curve shows no significant count rate oscillations. We have folded the light curve with different periods, ranging from 1000 and 5000 seconds, and have determined the corresponding $`\chi _\mathrm{r}^2`$ value. The best fitting period is found at 4200 seconds (cf. upper panel of Figure 2). The corresponding folded light curve is given in the middle panel of Figure 2. The power spectrum (lower panel) confirms the presence of a strong periodic signal with $`2.410^4\mathrm{Hz}`$. The second highest value is found at twice that frequency.

#### 2.1.2 Comparison with EPIC MOS and RGS data

The MOS and RGS data are basically consistent with the EPIC pn results showing a peak near 4200 seconds. In Figure 3 we show the folded light curves for the MOS-1 and RGS-1 detectors. Although the statistical errors are larger in the MOS and RGS data, the period of 4200 seconds is also confirmed by the MOS and RGS data.

#### 2.1.3 Monte-Carlo Simulations on red noise light curves

We made simulations of light curves with the observed time sequence and phase randomized for a red noise $`f^1`$ power spectrum (cf. Figure 4). The periodic signals detected in the X-ray light curve of Mrk 766 at 4200 and 2100 seconds (cf. lower panel of Figure 2) far exceed the power seen in the simulated power spectrum. Therefore we conclude that the periodicity peaks are not caused by red noise and are intrinsic to the Narrow-Line Seyfert 1 galaxy Mrk 766.

### 2.2 EPIC pn spectral analysis

The 0.2–10 keV EPIC pn spectrum of Mrk 766 is well described by a power law model with a photon index of $`\mathrm{\Gamma }=2.11\pm 0.10`$ dominating the 1 to 10 keV flux, in combination with a strong soft excess component which dominates the X-ray flux below about 1 keV. In Figure 5 we show the EPIC pn spectrum where the power-law photon index is extrapolated down to 0.2 keV.

The EPIC pn spectrum clearly indicates the presence of a strong soft excess component superimposed on an underlying power-law component. Between 0.2–0.7 keV the soft excess component contains about 40 per cent more flux than the power law component. The soft excess component is well described by a black body spectrum. Signatures for a strong warm absorber component are not detected in the EPIC pn spectrum. Although spectral residua between about 0.7–0.8 keV seem to be present in the RGS spectrum (Page et al. 2001), these residua cannot explain the spectral kink detected at about 1 keV in the EPIC pn spectrum (cf. Figure 7). We conclude that absorption due to a warm absorber, if present, is small compared to the soft excess emission and is unable to produce the spectral kink at around 1 keV. Therefore a combination of a black body component and an underlying power law is used for the spectral fitting described in the following. In Figure 6 we show EPIC pn count rate spectra selected according to 4 different phase intervals (middle panel of Figure 2). The count rate spectra were obtained from the phase intervals ranging from 0$``$0.2 (phase A), 0.2$``$0.55 (phase B), 0.55$``$0.85 (phase C), and 0.85$``$1.0 (phase D), respectively. These phase intervals correspond to different count rate states during the X-ray flux oscillations. The spectral fitting results indicate that the black body temperature remains constant within the errors for the 4 phase intervals. The black body temperature is 76 $`\pm `$ 1 eV (phase A), 76 $`\pm `$ 1 eV (phase B), 75 $`\pm `$ 1 eV (phase C) and 76 $`\pm `$ 1 eV (phase D). The absorption by neutral hydrogen is consistent with the Galactic value of $`N_{\mathrm{H},\mathrm{Gal}}=1.810^{20}\mathrm{cm}^2`$ (Dickey and Lockman 1990) for the four phase intervals. In addition, we show in Figure 7 the spectral fitting results for three different time intervals of the X-ray light curve shown in Figure 1. The spectra were selected for the high and low flux states for 3 periodicity peaks centered at 9000, 13200 and 29700 seconds, respectively. We find that the black body temperature remains constant within the errors during the X-ray oscillations for each periodicity peak.

To study the long-term variability we have divided the light curve shown in Figure 1 into 4 different time intervals, ranging from 7000$``$17100, 17100$``$24500, 24500$``$28500, and 28500$``$31000 seconds, respectively. For each time interval we have performed the combined power law plus black body fit and have determined the temperature and the bolometric luminosity of the black body component. With Figure 8 we show that the black body temperature increases from 74.6 to 80.1 eV and that bolometric luminosity increases from about $`3.5710^{44}\mathrm{erg}\mathrm{s}^1`$ to $`4.7210^{44}\mathrm{erg}\mathrm{s}^1`$ (for $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ and $`q_0=\frac{1}{2}`$).

## 3 Discussion

We first discuss the ’long-term‘ increase observed in the light curve (cf. Figure 1) which is characterized by a correlation between flux and temperature of the soft component. To our knowledge such a behavior has not been found in AGN spectra before. The correlation between the bolometric luminosity and the spectral temperature is close to the $`L_{\mathrm{bb}}T^4`$ law, consistent with a radiating surface remaining constant during the variations. The projected size of the black body radiating disk surface can be estimated using the luminosities and temperatures given in Figure 8. We find that this projected disk size is $`1.310^{25}\mathrm{cm}^2`$ which is rather a lower limit to the physical size of the emission region, because we have neglected a possible inclination of the disk, gravitational redshift effects and deviations from a black body due to scattering effects. We note that this area estimate requires a black hole mass of at least a few $`10^5`$ solar masses.

The periodic signal may be a quasi-periodic oscillation (QPO). Investigations of QPOs in galactic binary systems show a wide range of frequencies between about 0.1 and 1200 Hz (Psaltis et al. 1999). However, we note that for the galactic black hole microquasars GRS 1915+105 and GRO 1655-40 narrow frequencies of roughly 50 – 300 Hz have been detected (Morgan et al. 1997, Remillard et al. 1999). The ratio between these QPO frequencies and the strongest peak frequency found in the data of Mrk 766 lies between $`10^{56}`$ roughly consistent with the expected mass ratio of these objects. Specifically, QPOs could be produced by instabilities in the inner part of the accretion disk or by pulsating accretion if the rate is close to the Eddington limit. In these cases the period is expected to be of the order of the radial drift time scale (‘instability time scale’)

$$\tau ϵ\alpha ^1\left(\frac{r}{R_\mathrm{S}}\right)^{\frac{3}{2}}\frac{R_\mathrm{S}}{c}$$
$`(1)`$
where $`ϵ`$ 5 to 10 (Sunyaev 2000). We note that $`\tau `$ is of the order of a few thousand seconds for $`\alpha 0.1`$, $`r6R_\mathrm{S}`$ and $`M10^6`$ solar masses.

If the periodic signal arises from X-ray hot spots orbiting the black hole (e.g. Sunyaev 1973; Guilbert, Fabian & Rees 1983) we expect a variation in temperature. Using the radius dependent orbital velocity as given by Abramowicz et al. (1996), the boost model of Rybicki & Lightman (1979) and assuming an inclination angle of 45 degrees we find a temperature variation by a factor of about 1.6 and 1.4 at 3 and 4 Schwarzschild radii, respectively. As shown in Figure 7 the black body temperature remains constant within a few percent during the periodic variations requiring an inclination angle of the accretion disk of more than 85 degrees. For 85 degrees the temperature change according to the moving black body field is about a factor of 1.1 at 3 Schwarzschild radii, still inconsistent with the observations. Even for a nearly face-on geometry with an inclination angle of 89 degrees temperature changes by a factor of 1.06 are expected from the hot spot model.

Another possibility is that the observed variability is due to disk precession according to Bardeen and Petterson (1975). This causes the inner disk to settle into the plane of the black hole equator up to the so called transition radius, whereas the outer parts remain in their original orientation. The disk precession takes place in a ring which is in between. Nelson & Papaloizou (2000) found that the transition radius is of the order of about 10 to 30 gravitational radii. Since the precession of the disk is a geometrical effect it does not lead to spectral changes which is consistent with the observations. The absence of significant absorption by neutral gas above the galactic value shows that shadowing by the outer parts of the accretion disk is not responsible for the observed periodicity. The amplitude of the variations of about 7 per cent (see Figure 2) requires a precession angle dependent on the disk inclination according to $`\mathrm{\Delta }i=\frac{0.07}{\mathrm{cot}i}`$. For instance, for a disk inclination of 90 (45) degrees the precession angle is 32 (4) degrees. The angular shape of the light curve may reflect the expected disk warp. The second peak at 2100 seconds in the power spectrum and the possible third one at about 3000 seconds may be related to the warp as well. The precessing frequency produced by the Bardeen-Petterson effect depends on the disk radius $`r`$ as

$$\omega =\frac{2G}{c^2}Jr^3$$
$`(2)`$
where $`J`$ is the angular momentum of the black hole. The maximum angular momentum of a Kerr black hole is

$$J_{\mathrm{max}}=GM^2c^1$$
$`(3)`$
For the ratio between $`J`$ and $`J_{\mathrm{max}}`$ we get

$$\frac{J}{J_{\mathrm{max}}}=\frac{\omega c^3}{2G^2}\frac{r^3}{M^2}$$
$`(4)`$
Equ. 4 can be used to estimate the lower limit of the black hole mass in Mrk 766: Taking $`\omega =1.510^3`$ and the black body disk radius of $`210^{12}\mathrm{cm}`$ as an rough estimate for $`r`$, we find that a black hole mass of $`M>310^6`$ solar masses is required not to violate the angular momentum limit.

Presently, we are unable to differentiate between a strictly periodic behavior and quasi-periodic oscillations with high precision. We also need better data for the light curve and the power spectrum in order to constrain the models. In this context the forthcoming 150 ksec XMM-Newton observation will be most useful.

## 4 Summary

During the long-term variations of Mrk 766 the fitted black body temperature and luminosity of the soft component correlates as expected for a black body ($`LT^4`$) requiring a lower limit of the black body disk size of about $`10^{25}\mathrm{cm}^2`$. In addition, we have detected a strong periodic signal in the black body component with 4200 seconds. This periodicity may be of the QPO type resulting from instabilities in the inner accretion disk around the black hole. X-ray periodicity due to hot spots orbiting the black hole does not provide a plausible explanation, except for the case the disk is observed nearly face on ($`i`$ $``$ 85 degrees). As no significant periodic changes of the black body temperature are detected precession may provide another explanation for the X-ray periodicity of Mrk 766. If the Bardeen-Petterson effect is responsible for the observed periodicity, the required size of the emission region calls for a black hole mass of $`M>310^6`$ solar masses.

###### Acknowledgements.

It is a pleasure to acknowledge the efforts of the SOC and SSC teams in making the observations possible and for developing the SAS software package used to reduce the data. The XMM - Newton project is supported by the Bundesministerium für Bildung und Forschung / Deutsches Zentrum für Luft- und Raumfahrt (BMBF/DLR), the Max-Planck Society and the Heidenhain-Stiftung. We are grateful to the referee, K. Iwasawa, for very useful comments which helped us to improve the paper substantially. We are indebted to B. Aschenbach for fruitful discussions on the scientific content of this paper. TB thanks A.C. Fabian for using his software code on the red noise simulations.