text:
# 1 The temperature dependence of cos(𝛾) for 𝑛=3, 𝛼=0.001 and different values of the parameter 𝛽.

THE EXISTENCE OF A STABLE NONCOLLINEAR PHASE IN A HEISENBERG MODEL WITH COMPLEX STRUCTURE

Diana V. Shopova<sup>†,</sup>, and Todor L. Boyadjiev

CPCM Laboratory, G. Nadjakov Institute of Solid State Physics,
Bulgarian Academy of Sciences, BG-1784 Sofia, Bulgaria.

Faculty of Mathematics and Computer Science, University of Sofia,
BG-1164 Sofia, Bulgaria.

Corresponding author: email: sho@issp.bas.bg

Key words: Heisenberg model, magnetism, order parameter, phase transition.

PACS number(s):75.10.Hk, 75.30.Kz, 75.40.Cx

## Abstract

We have analyzed the properties of a noncollinear magnetic phase obtained in the mean-field analysis of the model of two coupled Heisenberg subsystems. The domain of its existence and stability is narrow and depends on the ratio between the averaged over nearest neighbours microscopic exchange parameters.

The study of Heisenberg magnets with complex structure is important for the understanding of magnetic properties of real substances . In this letter we shall present some results of the mean-field analysis based on the model of two bilinearly coupled classical Heisenberg subsystems. The proper identification of the possible magnetic phases, their domain of existence and stability can be described by this fully isotropic model for a big class of magnetic materials, in which the different types of anisotropies are several orders smaller than the exchange interaction. The importance of purely exchange models for the description of phases and phase transitions in magnetic substances is well elucidated, for example, in .

Our analysis is done with the help of the following microscopic hamiltonian:

$$H=\frac{1}{2}\underset{ij}{\overset{2N}{}}\left[J_{ij}^{(1)}𝐒_i^{(1)}𝐒_j^{(1)}+J_{ij}^{(2)}𝐒_i^{(2)}𝐒_j^{(2)}+2K_{ij}𝐒_i^{(1)}𝐒_j^{(2)}\right].$$
(1)
The variables $`𝐒_i^\alpha `$ are $`n`$-component classical spin vectors which obey the condition $`\left|𝐒_i^{(\alpha )}\right|=1`$; ($`\alpha =1,2`$). The dimensionless parametres,

$`J_{ij}^{(\alpha )}`$ $`=`$ $`{\displaystyle \frac{𝒥_{ij}^{(\alpha )}}{T}},`$ (2)
$`K_{ij}`$ $`=`$ $`{\displaystyle \frac{𝒦_{ij}}{T}},`$
where $`T`$ is the temperature, stand for the in-subsystems and intersubsystem exchange interactions, respectively. In the general case they are elements of $`(\text{}\times \text{ N})`$ matrix with $`N`$ \- the number of lattice sites in the subsystems.

The Ginzburg-Landau functional for the hamiltonian (1) is obtained in our previous paper  with the help of the Hubbard-Stratonovich transformation (HST). Somewhat different way of application of HST to the model (1) can be found in the book , where a RG analysis of the same model is done.

In the mean field calculations we neglect the spatial dependence of the order parameters and obtain the Landau free energy density in the form,

$`g{\displaystyle \frac{G}{NT}}`$ $`=`$ $`{\displaystyle \frac{\tau _1}{2}}𝝋_1^2+{\displaystyle \frac{\tau _2}{2}}𝝋_2^2+{\displaystyle \frac{g_1}{4}}\left(𝝋_1^2\right)^2+{\displaystyle \frac{g_2}{4}}\left(𝝋_2^2\right)^2+{\displaystyle \frac{v}{2}}𝝋_1^2𝝋_2^2`$
$`+v\left(𝝋_1𝝋_2\right)^2+w_1𝝋_1^2\left(𝝋_1𝝋_2\right)+w_2𝝋_2^2\left(𝝋_1𝝋_2\right).`$
The $`n`$-component real-space vectors $`𝝋_i,(i=1,2)`$, play the role of order parameters of the system. We shall consider only a positive in-subsystem exchange and the weak-coupling limit which means that the interaction between the subsystems is smaller than the in-subsystem interactions. As far as we do not take into account any magnetic anisotropies, we shall limit our considerations to magnetic substances with a cubic crystal structure as the bcc lattice with two different magnetic ions per unit cell. Under these conditions, the model (1) will describe two interpenetrating ferromagnetically ordered sublattices which interact either ferromagnetically or antiferromagnetically.

For a ferromagnetic in-subsystem exchange the coefficients in Eq. (3) will be expressed by the quantities $`S_0(𝒌),S_1(𝒌)`$ and $`\lambda _{1,2}(𝒌)`$ for $`𝒌=0`$, where $`S_0(𝒌)\text{and}S_1(𝒌)`$, are the elements of the unitary matrix, $`\widehat{S}`$, which diagonalizes the Ginzburg-Landau functional in $`𝒌`$space and $`\lambda _{1,2}(𝒌)`$ are the respective eigenvalues. The matrix elements and eigenvalues are functions of the averaged over the nearest neighbours microscopic exchange parameters (2) denoted by $`J_1,J_2,K`$, respectively, for details, see . The Landau free energy depends on the number of order parameter components $`n`$ through the parameter $`u=1/n^2(n+2)`$ which enters linearly in all coefficients in front of the fourth order terms. The coefficients in front of the second order terms $`\tau _{1,2}`$ are also functions of $`n`$; $`\tau _{1,2}=1\lambda _{1,2}/n`$.

The initial hamiltonian, Eq. (1), is degenerate with respect to the rotations of all spins by one and the same angle and this degeneracy is preserved in the Landau free energy density, Eq. (3). That is why from the mean-field equations:

$$\frac{g}{\phi _{1i}}=0,\frac{g}{\phi _{2i}}=0,(i=1,\mathrm{},n),$$

we can find for $`n>1`$, only the magnitudes of the vector order parameters $`𝝋_i,(i=1,2)`$ and their mutual orientation. Our previous analysis shows that concerning the mutual orientation of vector order parameters, there are two possibilities, namely:

$`\mathrm{cos}(𝝋_\mathrm{𝟏},𝝋_\mathrm{𝟐})`$ $`=`$ $`\pm 1\text{collinear phases}`$ (4)
$`\mathrm{cos}(𝝋_\mathrm{𝟏},𝝋_\mathrm{𝟐})`$ $``$ $`\pm 1\text{noncollinear phases}.`$
In the collinear case the system of equations for the magnitudes of vector order parameters $`𝝋_i,(i=1,2)`$ can be solved only numerically and the results are described in details in .

The respective equations for the noncollinear phase can be treated analytically. Because of the degeneracy of Landau free energy density in a purely exchange approximation for $`n>1`$, it is more appropriate to present $`𝝋_1`$ and $`𝝋_2`$ by the direction cosines $`\gamma _i`$ and $`\alpha _i`$ and their magnitudes $`\phi _1=\left(_{i=1}^n\phi _{1i}^2\right)^{1/2},\phi _2=\left(_{i=1}^n\phi _{2i}^2\right)^{1/2}`$, in the following way,

$`\phi _{1i}`$ $`=`$ $`\phi _1\gamma _i,`$ (5)
$`\phi _{2i}`$ $`=`$ $`\phi _2\alpha _i.`$
Then the Landau free energy density, Eq. (3) is rewritten with the help of Eq. (5) and afterwards minimized with respect to $`\alpha _i\gamma _i,\phi _1\text{and}\phi _2`$ under the conditions

$$\underset{i=1}{\overset{n}{}}\alpha _i^2=1,\underset{i=1}{\overset{n}{}}\gamma _i^2=1.$$

The resulting mean field equation for $`\mathrm{cos}(𝝋_1,𝝋_2)=_{i=1}^n\alpha _i\gamma _i`$ will be,

$$2v\phi _1\phi _2\underset{i=1}{\overset{n}{}}\alpha _i\gamma _i+w_1\phi _1^2+w_2\phi _2^2=0.$$
(6)
It is obvious from the above expression that the cosine between the vector order parameters is well defined for $`\phi _10`$ and $`\phi _20`$. The vector order parameters will be mutually perpendicular when $`w_1=0,w_2=0`$. This is possible for equivalent sublattices because $`w_{1,2}(J_1J_2)`$, but in this case the system of equations for $`\phi _1`$ ans $`\phi _2`$ has a solution only when $`K0`$, which means fully decoupled sublattices.

When $`J_1J_2`$ the system of equations for $`\phi _1,\phi _2`$ is simple and can be solved directly,

$`\tau _1+(g_1{\displaystyle \frac{w_1^2}{v}})\phi _1^2+(v{\displaystyle \frac{w_1w_2}{v}})\phi _2^2`$ $`=`$ $`0`$ (7)
$`\tau _2+(g_2{\displaystyle \frac{w_2^2}{v}})\phi _2^2+(v{\displaystyle \frac{w_1w_2}{v}})\phi _1^2`$ $`=`$ $`0.`$
To understand the properties of the obtained magnetic phase we need the connection between the order parameter vectors $`𝝋_i,(i=1,2)`$, and the sublattice magnetizations, $`𝒎_i,(i=1,2)`$

$`𝒎_1`$ $`=`$ $`{\displaystyle \frac{S_0}{\lambda _1}}𝝋_1{\displaystyle \frac{S_1}{\lambda _2}}𝝋_2,`$ (8)
$`𝒎_2`$ $`=`$ $`{\displaystyle \frac{S_1}{\lambda _1}}𝝋_1+{\displaystyle \frac{S_0}{\lambda _2}}𝝋_2.`$
From the above expression we can find the magnitudes of the sublattice magnetizations and the angle $`\gamma `$ between them. The requirements $`|\mathrm{cos}(\gamma )|1\text{and}|_{i=1}^n\alpha _i\gamma _i|1`$ restrict the values of parameters $`J_1,J_2\text{and}K`$, for which the noncolliner phase can exist.

The analytical expressions for magnitudes of the sublattice magnetizations and $`\mathrm{cos}(\gamma )`$ can be obtained straightforwardly after some algebra but they are very cumbersome, moreover, the temperature dependence of $`m_1,m_2`$ and $`\mathrm{cos}(\gamma )`$ cannot be comprehended directly from formulae. For this reason we prefer to give a graphical presentation of our results. In order to compare the properties of the noncollinear and collinear phases, obtained previously, we shall rewrite the solutions for the noncollinear phase with the help of the dimensionless parameters, introduced in :

$$x=\frac{T}{𝒥_1+𝒥_2}.$$

$$\alpha =\frac{𝒥_1𝒥_2}{𝒥_1+𝒥_2},\beta =\frac{2𝒦}{𝒥_1+𝒥_2}.$$

The quantity $`x`$ is called the reduced temperature; $`\alpha `$ can be considered as a measure of the difference between the two sublattices. We assume that, $`𝒥_1>𝒥_2`$, because the sublattices are symmetric with respect to the interchange $`12`$; then, $`0<\alpha <1`$. The dimensionless intersublattice interaction, $`|\beta |<1`$, which is a requirement of the weak-coupling limit $`𝒥_1𝒥_2𝒦^2>0`$.

In Fig. 1 the variation of $`\mathrm{cos}(\gamma )`$ with the reduced temperature $`x`$ is depicted for a fixed value of $`\alpha `$ and three different values of $`\beta `$. The behaviour of $`\mathrm{cos}(\gamma )`$ is essentially different for different values of the relation $`\alpha /\beta `$. When $`\alpha \beta `$, the noncollinear phase is metastable in the whole domain of its existence and $`cos(\gamma )<1`$, including the phase transition point. In the opposite case, i.e, $`\alpha <\beta `$, the noncolliner phase becomes the most stable phase in the temperature interval $`x_n<xx^{},\text{where}x^{}`$ is the transition temperature of the noncollinear phase and $`x_n`$ is the temperature, below which the noncollinear phase does not exist. For $`\beta =0.003`$, the angle between the sublattice magnetizations vary from zero at $`x^{}`$ to $`68^{}\text{at}x_n`$. When $`\beta `$ grows (see the curve for $`\beta =0.005`$) the domain of the existence for the noncollinear ferrimagnetic phase becomes more narrow and for $`\beta 0.0056`$ it does not exist.

The calculations show that in the region of stability of the noncollinear phase the magnetization $`𝑴=𝒎_\mathrm{𝟏}+𝒎_\mathrm{𝟐}`$ is perpendicular to the staggered magnetization $`𝑳=𝒎_\mathrm{𝟏}𝒎_\mathrm{𝟐}`$. This is the condition for the existence and stability of a noncollinear ferrimagnetic phase, derived on the basis of a symmetry group analysis by Andreev and Marchenko .

In Figs. 2 and 3 we compare the Landau free energy density as a function of the reduced temperature $`x`$ for all magnetic phases we have found (collinear and noncollinear).

The phase with the highest transition temperature, denoted by $`S`$ in Figs. 2, 3, has a collinear ferrimagnetic structure with a compensation point at lower temperatures. It is stable and occurs through a second order phase transition. The symbols for the high- and low-temperature metastable collinear phases are $`M_i,i=1,2`$, respectively. $`M_3`$ is a collinear metastable phase, which occurs by a second order phase transition, not given in . The Landau energy density of the noncollinear phase is denoted by $`New`$.

For $`\alpha >\beta `$ the noncollinear phase has an energy slightly higher than the stable collinear phase - Fig. 2. When the difference between the interaction in the sublattices becomes smaller than the intersublattice interaction, the noncollinear phase has lower energy than the $`S`$-phase in the whole domain of its existence - Fig. 3.

We can argue that with the decrease of the temperature the system orders by a second order phase transition in the stable collinear ferrimagnetic phase and with the further lowering of temperature a second phase transition occurs and a new ordered noncollinear phase appeares. At lower temperatures the noncollinear phase ceases to exist. It is obvious from the above figures that the domain of existence and stability of the noncollinear phase is very narrow and it can occur only for magnetic substances, in which the following relation between the averaged microscopic exchange parameters holds true: $`(𝒥_1𝒥_2)<𝒦<\sqrt{𝒥_1𝒥_2}`$.

Up to now we have considered only the ferromagnetic exchange between the sublattices ($`\beta >0`$). As the calculations show the respective noncollinear phase for $`\beta <0`$, i.e., for an antiferromagnetic coupling, exhibits a similar temperature behaviour as the ferrimagnetic noncollinear phase described above.

We did not study the noncollinear phase in the general case of $`n`$-component order parameters. The above calculations are made for $`n=3`$ which is the most interesting case up to our opinion. The analysis shows that for $`n=2`$ the domain of existence of the noncollinear phase is very tiny and we did not find values of the parameters $`\alpha \text{and}\beta `$, for which it can be stable.

The difference between the energies of the stable collinear ferrimagnetic and the noncollinear ferrimagnetic phases is very small. In the real magnetic substances even a small anisotropic interaction may change the above described picture, but this problem needs a separate investigation.

Acknowledgements

The authors thank Dimo Uzunov for the valuable discussions.


format: text

text:
# Three-dimensional calculations of high and low-mass planets embedded in protoplanetary discs

## 1 Introduction

Young planets interact with the surrounding discs from which they form by accreting mass and exerting torques. In the core-accretion model for planet formation, the initial growth proceeds as solids accumulate to form a planetary core (Lissauer 1995, Wuchterl, Guillot & Lissauer 2000). The core may become a terrestrial planet or may further develop to become a gas giant through accretion of gas. The torques resulting from planet-disc interactions cause a planet to migrate typically inwards (Lin et al. 2000; Ward and Hahn 2000).

Numerical hydrodynamical calculations have been carried out in two dimensions in order to understand better the dynamics of a circular orbit planet embedded in a gaseous disc (e.g. Bryden et al. 1999; Kley 1999; Lubow, Seibert & Artymowicz 1999; Nelson et al. 2000). These studies concentrated on the interaction between a 1 MJ planet and a gas disc that orbit a 1 M star. The tidal forces caused by the planet create a gap in the disc. In spite of the presence of the gap, disc mass flow on to the planet continues through the gap with high efficiency. Nearly all the flow through the gap is accreted by the planet at a rate comparable to the rate at which mass accretion would occur in the disc in the absence of a planet.

The flow within the Roche lobe of the planet is highly non-axisymmetric and involves shocks produced by colliding gas streams (Lubow et al. 1999; D’Angelo, Henning & Kley 2002). Torques on the planet are exerted by circumstellar disc material that lies outside the gap. In addition, torques are exerted locally by the material that flows close to the planet. The net torque results in inward migration. Gas accretion can continue to planet masses of order 10 MJ, at which point tidal forces are sufficiently strong to prevent flow into the gap.

These earlier studies were limited by their neglect of effects in the vertical direction (perpendicular to the orbit plane). The Roche lobe radius of a 1 MJ planet orbiting a 1 M star is comparable to the local disc thickness. The effects in the vertical direction are even more important for lower-mass planets. In addition, sufficiently low-mass planets do not open a gap in the disc. Analytic models of planet migration in the non-gap case, sometimes called Type I migration, were carried out in two dimensions (Ward 1997). The results implied that planetary migration timescales are much shorter than the disc lifetimes. Planets would then be accreted by the central star. Alternatively, the planets could reside at the circumstellar disc inner edge, if there is a central hole in the disc. This situation poses problems for planet survival and for gas giant planet formation. Recent two and three dimensional analytic calculations by Tanaka, Takeuchi, and Ward (2002) of the migration rates of low-mass planets obtained smaller values (by about an order of magnitude), making survival more plausible.

The two-dimensional numerical results also indicated that accretion within the Roche lobe of a 1 MJ planet is driven by shocks (Lubow et al. 1999; D’Angelo et al. 2002). It is not clear whether the same flow pattern would persist in three dimensions. Recently, Kley, D’Angelo, and Henning (2001) computed the three-dimensional flow for 0.5 and 1 MJ planets and found minor differences in the accretion and migration rates compared to the two-dimensional case. However, the flow within the Roche lobe was unresolved.

In this paper, we analyze the non-linear interactions between a planet and a gaseous disc by using global three-dimensional numerical simulations that resolve the flow within the Roche lobe of a 1 MJ planet. We are interested in determining the flow patterns, the accretion rates, and migration rates for planets whose masses range between 1 M (Earth mass) and 1 MJ (Jupiter mass). Each planet is assumed to be in a circular orbit about the central star. The mass and orbital radius of each planet are fixed during the simulation; consequently, the results do not include the effects on the flow of planetary migration.

The outline of the paper is as follows. Section 2 describes the computational procedure. Section 3 provides the results. In Section 4, we discuss the implications of the results for giant planet formation. Section 5 contains our conclusions.

## 2 Computational method

### 2.1 Basic equations

We use a computational method that is similar to that of Lubow et al. , except that they performed two-dimensional vertically averaged calculations whereas we solve the problem in three dimensions.

We assume a viscous model for the disc turbulence, with the usual $`\alpha `$ prescription. The origin of the coordinate system is taken to be the centre of the star. (We ignore the slight centre-of-mass shift caused by the planet.) The disc self-gravity is ignored. The flow is modelled in the orbital frame that rotates with the angular speed of the planet $`\mathrm{\Omega }_\mathrm{p}=\sqrt{GM_{}/r_\mathrm{p}^3}`$, where $`M_{}`$ is the mass of the star, $`r_\mathrm{p}`$ is the orbital radius of the planet, $`G`$ is the gravitational constant, and we have neglected the mass of the planet. In this frame, the flow achieves a near steady state. We adopt spherical coordinates $`(r,\theta ,\varphi )`$ with associated flow velocities in the rotating frame $`𝐮=(u_\mathrm{r},u_\theta ,u_\varphi )`$. The equations of motion for the disc are

$$\frac{\rho }{t}+(\rho 𝐮)=0,$$
(1)
$`{\displaystyle \frac{s_\mathrm{r}}{t}}+(s_\mathrm{r}𝐮)`$ $`=`$ $`\rho r\mathrm{sin}^2\theta \left({\displaystyle \frac{u_\varphi }{r\mathrm{sin}\theta }}+\mathrm{\Omega }_p\right)^2+{\displaystyle \frac{\rho u_\theta ^2}{r}}`$ (2)
$`{\displaystyle \frac{p}{r}}\rho {\displaystyle \frac{\mathrm{\Phi }}{r}}+f_r,`$
$`{\displaystyle \frac{s_\theta }{t}}+(s_\theta 𝐮)`$ $`=`$ $`\rho r^2\mathrm{sin}\theta \mathrm{cos}\theta \left({\displaystyle \frac{u_\varphi }{r\mathrm{sin}\theta }}+\mathrm{\Omega }_p\right)^2`$ (3)
$`{\displaystyle \frac{p}{\theta }}\rho {\displaystyle \frac{\mathrm{\Phi }}{\theta }}+rf_\theta ,`$
and

$$\frac{s_\varphi }{t}+(s_\varphi 𝐮)=\frac{p}{\varphi }\rho \frac{\mathrm{\Phi }}{\varphi }+r\mathrm{sin}\theta f_\varphi ,$$
(4)
where $`\rho `$ is the gas density, $`p`$ is the gas pressure, $`s_\mathrm{r}=\rho u_\mathrm{r}`$ is the radial momentum per unit volume, $`s_\theta =\rho ru_\theta `$ is the meridional momentum per unit volume, $`s_\varphi =\rho r\mathrm{sin}\theta (u_\varphi +\mathrm{\Omega }_pr\mathrm{sin}\theta )`$ is the azimuthal angular momentum per unit volume, $`\mathrm{\Phi }`$ is the gravitational potential due to the central star and the planet, and $`𝐟=(f_\mathrm{r},f_\theta ,f_\varphi )`$ is the viscous force per unit volume that describes the effects of disc turbulence. We use an unsoftened gravitational potential

$$\mathrm{\Phi }(𝐫)=\frac{GM_{}}{r}\frac{GM_\mathrm{p}}{|𝐫𝐫_\mathrm{p}|},$$
(5)
where $`M_\mathrm{p}`$ is the planet mass. This is possible because the location of the planet is such that all gravitational force evaluations are made at a finite distance from the planet. We also performed some test calculations with gravity softened on the length scale of twice the grid resolution near the planet, but found no significant difference between the softened and unsoftened results.

Equations (1), (2), (3), and (4) express conservation of mass, radial momentum, meridional momentum, and azimuthal angular momentum, respectively. Equation (4) is written in terms of the total azimuthal angular momentum $`s_\varphi `$, rather than that in the rotating frame. The reason is that the $`s_\varphi `$ equation provides better numerical stability .

The equation of state is taken to be locally isothermal, $`p\rho T`$, with the temperature expressed as a specified function of radius, $`T(r)`$. This equation of state is appropriate for a gas that radiates internal energy gained by shocks with high efficiency. The viscosity force $`𝐟`$ is assumed to be the standard Navier-Stokes force (see eq. \[15.3\] of Landau & Lifshitz 1975; Klahr, Henning & Kley 1999). The coefficient of shear viscoity $`\mu `$ represents the effects of disc turbulence, while the bulk viscosity coefficient $`\zeta `$ is set to zero. The value for the kinematic turbulent viscosity $`\nu =\mu /\rho `$ is assumed to be constant in space and time. It can be expressed in terms of the usual $`\alpha `$ prescription of Shakura & Sunyaev . Namely, for a disc with local isothermal sound speed $`c_\mathrm{s}`$ and vertical scaleheight $`H`$, dimensionless parameter $`\alpha `$ is defined through

$$\alpha (r)=\frac{\nu }{c_\mathrm{s}H}.$$
(6)
The above equations are non-dimensionalized so that the unit of time is the inverse of the planetary orbital frequency $`\mathrm{\Omega }_\mathrm{p}`$, the unit of distance is the orbital radius of the planet $`r_\mathrm{p}`$, and $`G=1`$.

### 2.2 Numerical method

The equations are solved using a three-dimensional spherical coordinate version of the ZEUS-2D code . The code was written by K. A. Miller and J. M. Stone. It was modified to include a standard three-dimensional Navier-Stokes viscous force term. The code allows for variably-spaced gridding, which permits us to obtain higher resolution in the vicinity of the planet. The timesteps satisfy the usual Courant condition, for which we have adopted Courant number 0.3. The code provides an artificial viscosity term, but since we introduce a Navier-Stokes viscous force, we suppress this artificial viscosity. Of course, there is some intrinsic numerical viscosity because of the finite gridding. For a uniformly spaced mesh, the code is formally second-order accurate in space and first-order accurate in time. For variably spaced meshes, it is formally first-order accurate in space. However, a high level of accuracy and resolution can be attained by limiting the fractional change in mesh spacing between adjacent cells to be small, of order 1%. The code uses van Leer interpolation.

### 2.3 Numerical grid and initial conditions

The planet was fixed at location $`(r,\theta ,\varphi )=(1,\pi /2,\pi )`$. We modelled the disc in the region $`r[0.3,4.0]`$, $`\theta [\pi /24H/r,\pi /2]`$, and $`\varphi [0,2\pi ]`$. We imposed reflective boundary conditions at the radial and $`\theta `$ grid boundaries and periodic boundary conditions at the azimuthal boundary. The radial boundaries were sufficiently far from the planet that the reflected waves were not noticeable. The grid was uniform in $`\theta `$, but non-uniform in $`r`$ and $`\varphi `$. Both the $`r`$ and $`\varphi `$ grids were uniform in the vicinity of the planet, but the grid spacing increased logarithmically away from this uniform region (i.e. each zone was a certain percentage larger or smaller than the zone preceding it). In addition, the last two zones inside $`\varphi =0`$ and $`\varphi =2\pi `$ were uniform to allow periodic boundaries to be implemented easily. Each region of the radial grid had $`[0.3,0.8965]=28`$, $`[0.8965,1.1035]=72`$, $`[1.1035,4.0]=80`$ zones. The $`\varphi `$ grid had $`[0,0.19]=[2\pi 0.19,2\pi ]=2`$, $`[0.19,\pi 0.1035]=[\pi +0.1035,2\pi 0.19]=106`$, $`[\pi 0.1035,\pi +0.1035]=72`$ zones in each region. The $`\theta `$ grid modelled 4 scale heights above the disc midplane using 36 zones. Thus, the zones in the vicinity of the planet had dimensions of 0.002875 in $`r`$ and $`\varphi `$, and twice this value in $`\theta `$. This grid was decided upon after extensive testing using two and three-dimensional calculations. In particular, we found that if the change in size between two neighbouring zones is too large or if the $`r`$ and $`\varphi `$ dimensions of zones are too different, numerical instabilities occur. The instabilities are worst in two-dimensional calculations, but also appear in three-dimensional calculations that have the same gridding in $`r`$ and $`\varphi `$. We also performed two-dimensional calculations with increasing resolution in $`r`$ and $`\varphi `$ to test for convergence. Lowering the resolution in either $`r`$ or $`\varphi `$ tends to increase the accretion rates. Convergence testing of three-dimensional calculations was not practical due to the factor of 16 increase in computational time that would have been required to double the above resolution in all three directions. The above resolution was chosen as the minimum required both to avoid numerical instabilities and to give quasi-steady-state accretion rates that had converged to within a few percent in the two-dimensional calculations.

The Roche lobe of the planet had a radius

$$r_\mathrm{R}=\left(\frac{M_\mathrm{p}}{3M_{}}\right)^{1/3}r_\mathrm{p}.$$
(7)
We modelled planets from 1 Earth mass ($`3\times 10^6`$ M) to 1 Jupiter mass ($`1\times 10^3`$ M). With the mass of the star equal to 1 M, the planets’ Roche radii ranged from $`r_\mathrm{R}/r_\mathrm{p}=0.010`$ to 0.069. Thus, the Roche radius of the Jupiter-mass planet was resolved by 24 zones (i.e. the Roche lobe contains about $`3\times 10^4`$ zones), while that of the Earth-mass planet was marginally resolved by 3.5 zones (i.e. the Roche lobe contains about $`80`$ zones).

The planet was assumed to be able to accrete material without substantial expansion of its radius on the scale of its Roche lobe. Some models of planet evolution suggest that this assumption is valid only for planets of mass greater than about 50 M (Pollack et al 1996). However, these models are subject to number of major simplifying assumptions, as discussed in Wuchterl et al. (2000). We return to this point in Section 4. In any case, the accretion assumption provides a simple prescription for handling the accretion flow that allows us to compute torques on the planet due to its interaction with the disc.

To simulate the accretion on to the planet, we nearly fully removed material in the four grid zones that surrounded the location of the planet in each timestep. A residual density was retained in these zones to avoid numerical divergences. The removed material was assumed to be accreted on to the planet, although the mass of the planet was not increased. There was a pressure force directed toward the planet at the edge of the evacuated region. This force was small compared with other dynamical forces for massive planets, but was noticable for the lowest mass planets that we modelled (see Section 3.2).

A temperature profile $`T(r)r^1`$ was used at all times. It was normalised such that $`H/r_\mathrm{p}=c_\mathrm{s}/(\mathrm{\Omega }_\mathrm{p}r_\mathrm{p})=0.05`$. This profile gives a constant value of $`H/r`$ throughout the unperturbed disc. It is numerically convenient when using a spherical grid, since the vertical resolution of the disc is constant. We set the kinematic disc viscosity to $`\nu =10^5`$ in our dimensionless units, which corresponds to $`\alpha =4\times 10^3`$ at $`r=r_\mathrm{p}=1`$ (equation 6).

The underlying initial disc density profile was chosen to be axisymmetric and follow $`\rho (r,\theta ,\varphi )r^{3/2}\mathrm{exp}[(\theta \pi /2)^2r^2/(2H^2)]`$. Thus, the disc surface density $`\mathrm{\Sigma }(r)r^{1/2}`$. For planets with masses $`M_p0.1`$MJ, we imposed an initial gap near the planet. The initial gap size and structure were estimated by an approximate torque balance condition between viscous and tidal torques near the planet. For the Jupiter-mass planet, the density at the middle of the gap was 1% of the unperturbed density. For the 0.3 MJ and 0.1 MJ planets the initial gap densities were 4.2% and 35% of the unperturbed density, respectively.

In order to confirm that the presence of an initial gap did not affect the final results, we performed two calculations of the 0.1 MJ planet, one with an initial gap and one without. The accretion and migration rates of the calculation with the initial gap reached steady values after about 30 orbits. The calculation without an initial gap initially had much higher accretion and migration rates, but these converged to the values given by the other calculation after approximately 100 orbits. Thus, the presence of an initial gap for the high-mass planets does not affect the final results, but it significantly reduced the computational time.

## 3 Results

### 3.1 Calculations

We performed calculations of 6 planets with masses of 1, 0.3, 0.1, 0.03, 0.01, and 0.003 Jupiter masses, MJ (i.e. 330, 100, 33, 10, 3.3, and 1 Earth masses, M). The results are scaled so that the planet is at a distance of 5.2 AU from a 1 M star. The disc mass, $`M_\mathrm{d}`$, between 1.56 and 20.8 AU (the boundaries of the grid) is taken to be 7.5 Jupiter masses (i.e. 0.0075 M). This gives an unperturbed disc surface density of 75 g cm-2 at the radius of the planet.

The calculations were run until the accretion rate on to the planet reached a quasi-steady state (Figure 1). The number of orbits required depended on the mass of the planet. Lower-mass planets reached a steady state after only about 5 orbits, while the 1 MJ planet required approximately 100 orbits. The calculations never quite reach a true steady state because the protoplanetary disc is evolving due to the Navier-Stokes viscosity. Thus, the accretion rates of some of the high-mass planets can be seen to decrease slightly with time after the quasi-steady-state accretion is reached.

The calculations were performed on the United Kingdom Astrophysical Fluids Facility (UKAFF), a 128-processor SGI Origin 3800 computer, and on GRAND, a 24-processor SGI Origin 2000 computer. The computational time required for each calculation was relatively independent of the planet’s mass and required approximately 160 CPU hours per orbit on UKAFF. The total CPU time required for all the calculations was approximately 50000 CPU hours.

### 3.2 Density structure and gas flow

The interaction of the planet with the disc alters the density and flow of the gas in the vicinity of the planet. This interaction was analysed by means of two-dimensional simulations for high-mass planets by Lubow et al. . Recently, D’Angelo et al.  studied the gas flow near low-mass planets by using two-dimensional simulations. Here, we study the problem in three dimensions for both high and low-mass planets.

#### 3.2.1 Global features

Figure 2 plots the azimuthally averaged surface density as a function of radius for the planets with masses ranging from 1 to 0.01 MJ. Also plotted (with the thick solid line) is the two-dimensional result for a 1 MJ planet after 150 orbits from Figure 1 of Lubow et al. (1999). The surface density in the vicinity a 1 MJ planet is well modelled by two-dimensional calculations (see also Kley et al. 2001). Only planets with masses $`M_\mathrm{p}\text{ }>0.1`$ MJ produce significant perturbations in the disc’s surface density; $`0.03`$ MJ (10 M) is insufficient. This result is consistent with the expectation that a gap will start to open when the Roche radius of the planet is comparable to the disc scaleheight $`H`$. The Roche radii of planets with masses 0.3, 0.1, and 0.03 MJ are $`r_\mathrm{R}/r_\mathrm{p}=0.046,0.032,`$ and 0.022, respectively. Comparing these to $`H/r=0.05`$, we find that $`r_\mathrm{R}>H/2`$ is required for a significant surface density perturbation. Even for $`H=r_\mathrm{R}`$ (the 0.3 MJ case), the gap only drops to 15% of the unperturbed density.

In Figure 3, we plot the global surface density of the disc at the end of the simulations for each of the 6 planets. The (partial) clearing of gaps by the planets is clearly visible in the global surface density plots for planets with masses $`M_\mathrm{p}0.1`$ MJ. Spiral density waves launched by the planet propagate inwards and outwards from the planet’s radius into the disc. For masses $`M_\mathrm{p}0.1`$ MJ, these density waves are strong enough to appear even in the azimuthally averaged surface density profiles (Figure 2). The strength of the surface density perturbations decreases with lower mass planets, but their form is essentially independent of the planet’s mass. For the lowest-mass planet, the spiral density perturbations are too small to be visible over the underlying density gradient in the disc. These waves are essentially two-dimensional in nature and have been studied in the past using two and three-dimensional nonlinear calculations (e.g. Artymowicz 1992; Kley 1999; Bryden et al. 1999; Lubow et al. 1999; Nelson et al. 2000; Kley et al. 2001) as well as analytically and through linearized numerical calculations (Ogilvie & Lubow 2002).

#### 3.2.2 The flow in the vicinity of the planet

In Figure 4, we plot the surface density in the vicinity of the planet. Overlaid on the plots are the Roche lobes of the planets and their inner (L1) and outer (L2) Lagrangian points (crosses). In Figure 5, we plot the density of the disc (greyscale) and streamlines on the disc midplane in the vicinity of the planet. From these figures, we see that the waves are present in the circumstellar discs, as was found in the previous two-dimensional simulations. They are initiated near the planet and propagate as shocks. Even with the 1 M planet, the density jump across the shock is of order 10%, and the more massive planets have greater density jumps. The discontinuities in the streamlines across the shocks are obvious for planets with masses $`M_\mathrm{p}10`$ M.

The shape of the shocks is independent of the planet’s mass so long as a deep gap is unable to form (i.e. for planets with masses $`M_\mathrm{p}\text{ }<0.1`$ MJ). For the higher-mass planets, the shocks become more curved, the outer shock moves forward, and the inner shock moves back. Unlike the 1 MJ case, the shocks generated by the low-mass planets do not extend radially all the way to the planet’s Roche lobe. Instead, their radial extent is determined by the disc thickness. The shocks reach no closer to the planet than a distance of approximately $`H`$ (i.e. they begin at approximately at $`r=r_\mathrm{p}\pm H`$). This fact will be important when we come to consider the torque exerted on the planet by the disc. Although the disturbances are non-linear, their form is similar to the surface density perturbations expected from linear resonance theory (see Figure 6 of Tanaka et al. 2002).

Figure 6 provides the density of the gas in the vicinity of the planet, but at height $`H`$ above the midplane. Figure 7 gives the density and velocity vectors in an $`rz`$ slice through the disc at the location of the planet (i.e. $`\varphi =\pi `$). The results are plotted by using a Cartesian coordinate system $`(x,y,z)`$, such that the the planet is located at $`(1,0,0)`$ and disc midplane is defined by $`z=0`$.

Notice that we have not plotted streamlines in these figures. It is difficult to plot a set of three-dimensional streamlines because they cross each other in projection on to a graph. Streamlines can be plotted at the disc midplane where the vertical velocity is zero (by symmetry). The velocity vectors in Figure 7 cannot be connected to form streamlines because information about the velocity in the deprojected direction ($`y`$) is required. Consequently, one cannot conclude that material rapidly drops to the midplane near $`x=1`$. Instead, material off the midplane flows past a low-mass planet in the $`y`$direction.

These figures allow us to determine the structure of the shocks outside the Roche lobe of the planet. Comparing the locations of the shocks in similar panels between Figures 5 and 6, we find that the locations are similar for the low-mass planets but, for 1 and 0.3 MJ, the shock fronts on the midplane lead the shock fronts off the midplane. Thus, a section perpendicular to the shock front for the high-mass planets would have a ‘bowshock’ shape, whereas for the low-mass planets the shocks are essentially vertical planes. These shock shapes can be seen clearly in Figure 7 with the transition from nearly vertical (0.1 MJ) to bowshock-shaped (1 MJ). The discontinuity in the velocity vectors across these shocks is also clear in Figure 7.

Returning to the density structure and streamlines at the midplane (Figure 5), we see that gas at the same radius as the planet but outside the planet’s Roche lobe moves on horseshoe orbits, as was reported in earlier studies (e.g., Bryden et al. 1999; Kley 1999; Lubow et al. 1999). For the high-mass planets, these horseshoe orbits occupy part of the gap in the disc. The radial extent of the horseshoe orbits decreases as the mass of the planet is decreased.

Between the outer disc and the portion of the horseshoe orbits ahead of the planet, and between the inner disc and the portion of the horseshoe orbits behind the planet, there are two streams of material that enter the Roche lobe of the planet. This is the material (along with material from above and below the disc midplane that is harder to visualise) that is accreted by the planet (see also Lubow et al. 1999). For planets with masses $`M_\mathrm{p}0.1`$ MJ, we notice that the breadth of these streams generally increases as the planet’s mass increases, with the breadth of the stream being of order the Roche radius of the planet. In Section 3.3, we use this observation to develop a model for the accretion rate of the low-mass planets. For the 1 and 0.3 MJ planets, the streams are significantly narrower than the planet’s Roche lobe. In fact, for the 0.3 MJ planet, no streamlines on the midplane enter the planet’s Roche lobe. Since this planet has the second highest accretion rate among the cases we consider, material is still being accreted, but the flow must be intrinsically three-dimensional with the accretion coming from above and below the midplane. The 0.3 MJ planet may be a special case, since the its Roche radius $`r_\mathrm{R}=0.046`$ is almost identical to the unperturbed scaleheight of the disc. However, even with the 1 MJ planet, the streams entering the Roche lobe on the midplane are significantly narrower than in the two-dimensional calculations of Lubow et al. , again implying that the flow into the planet’s Roche lobe is three-dimensional. To illustrate further this difference between the two-dimensional and three-dimensional accretion streams on to a Jupiter-mass planet, we have performed a two-dimensional 1 MJ calculation with identical resolution in $`r`$ and $`\varphi `$ to our three-dimensional calculation. The density and streamlines in the vicinity of the planet are given in Figure 8 and should be compared to the top-left panel of Figure 5. The accretion streams are much broader in the two-dimensional case.

Finally, we note that flow near the 1 M planet is only marginally resolved, since the evacuated zones at the location of the planet can be seen to affect the gas flow at the boundary of the Roche lobe. Whereas the critical streamlines that mark the boundaries of the accretion streams graze the Roche lobe of the 3.3 and 10 M planets, gas is sucked into the Roche lobe of the 1 M planet by the artificial pressure gradients. This effect leads to a slight over-estimate of the accretion rate for the 1 M planet (Section 3.3).

#### 3.2.3 The circumplanetary disc

Once inside the planet’s Roche lobe, the gas settles into a circumplanetary disc. In our calculations, the discs are resolved only for planets with masses $`M_\mathrm{p}0.1`$ MJ; the resolution within the Roche lobe of lower-mass planets is insufficient to follow material around complete orbits.

The scaleheight of the circumplanetary disc is much smaller than that of the circumstellar disc, $`H`$, due to the planet’s vertical gravity. This can be deduced from the fact that there is no sign of the discs in the plots of density at height $`H`$ above the midplane (Figure 6). Furthermore, Figure 7 clearly shows the cross sections of the discs for the 1 and 0.3 MJ cases. As in two-dimensional calculations, the discs rotate in a prograde manner (Lubow et al. 1999; D’Angelo et al. 2002). These circumplanetary discs may lead to satellite formation around the planets.

Lubow et al.  studied the flow inside the Roche lobe of an accreting planet. In two-dimensional calculations of a Jupiter-mass planet, they found that the circumplanetary disc contained strong shocks that led to rapid accretion of the gas through the planet’s disc. They concluded that the shocks formed from the collision of the two streams passing into the Roche lobe, one from the inner disc and one from the outer disc. D’Angelo et al.  studied lower mass planets with a two dimensional nested grid code and found these shocks persisted down to planet masses of approximately $`5`$ M.

We find that such strong shocks in the circumplanetary disc do not occur in three-dimensional calculations. They are an artifact of the two-dimensional calculations in which the streams from the inner and outer discs are forced to collide without vertical motions. The 0.3 MJ case is a particularly good example of the difference between two and three-dimensional structures, since in three dimensions there are no streamlines at the midplane that even enter the planet’s Roche lobe.

The disappearance of the spiral shocks in the circumplanetary disc when going from two to three-dimensional simulations is further illustrated for a 1 MJ planet by comparing the density and streamlines in a two-dimesional calculation (Figure 8) with those in the three-dimensional calculation (Figure 5, top-left panel). Stream collisions do occur at the midplane in both cases. But the strong spiral shocks present in the two-dimensional calculation are greatly diminished in three dimensions. This weakening in three dimensions is evidenced by the fact that streamlines within the Roche lobe are more tightly wrapped in three-dimensional case. In two dimensions, the spiral shocks drive the accretion through the circumplanetary disc with the gas losing angular momentum on each passage through a shock.

We make two suggestions for the effects of the absence of strong spiral shocks in the three-dimensional circumplanetary discs. First, accretion through the circumplanetary disc may be driven in a similar manner to accretion through the circumstellar disc, rather than by spiral shocks. Second, the more quiescent three-dimensional flow might be more conducive to satellite formation in circumplanetary discs.

### 3.3 Accretion rates

In Figure 9, we plot the final accretion rates in Jupiter masses per year. Lubow et al. (1999) performed two-dimensional calculations to determine the relative accretion rates of planets with masses $`M_\mathrm{p}1`$ MJ. The open circles in Figure 9 give their accretion rates, scaled to match our Jupiter-mass case. Two-dimensional calculations are adequate for modelling accretion on to high-mass planets which open well-defined gaps in the disc (Kley et al. 2001). Figure 9 then gives accretion rates over 3.3 orders of magnitude in planet mass, from 1 Earth mass to 6 Jupiter masses. As discussed in Section 3.2.2, the accretion rate on to the 1 M is slightly over-estimated because of the lack of resolution in the vicinity of the planet’s Roche lobe. Thus, the derived accretion rate is plotted as an upper limit in Figure 9.

For low-mass planets up to $`M_\mathrm{p}10`$ M (0.03 MJ), the accretion rate is proportional to the planet’s mass. The accretion rates peak at $`M_\mathrm{p}0.1\text{M}\text{J}`$, just as the planet starts to open a gap in the disc (Figure 2). For higher masses, the accretion rate drops rapidly with increasing mass as the gap gets wider (Figure 2). As discussed by Lubow et al. , the rapid decline of accretion rate at high masses provides a natural limit of about $`10`$ MJ for the mass of a planet on a circular orbit (see Section 4).

The gas mass accretion process is always dominated by three-body effects and the relevant capture radius is the Roche lobe radius $`r_\mathrm{R}`$. This statement can be justified by considering the ratio of the Roche lobe radius to the Bondi-Hoyle accretion radius which is approximately $`(H/r_\mathrm{R})^2`$. This ratio is greater than unity for planets whose mass is less than 1 MJ for typical disc parameters ($`H/r0.05`$). The mass capture rate for a low-mass planet can then be estimated by a simple argument that $`\dot{M}_\mathrm{p}\pi r_\mathrm{R}^2\rho v`$, where $`\rho `$ is the gas density along the planet’s orbit and we have used the fact that the breadth of the accretion streams into the planet’s Roche lobe scales with the Roche radius, $`r_\mathrm{R}`$ (Section 3.2.2). But since the velocity of the gas relative to the planet $`v\mathrm{\Omega }_pr_\mathrm{R}`$, using equation 7, we can express the accretion rate as

$$\dot{M}_\mathrm{p}=b\frac{M_\mathrm{p}}{M_{}}\rho \mathrm{\Omega }_\mathrm{p}r_\mathrm{p}^3,$$
(8)
where $`b`$ is a constant of order unity. Consequently, we recover the numerical result that $`\dot{M}_\mathrm{p}M_\mathrm{p}`$, with $`b=2.30`$ in the above equation (dotted lines, Figure 9).

The mass accretion efficiency was defined in Lubow et al. (1999) as

$$=\frac{\dot{M}_\mathrm{p}}{3\pi \nu \mathrm{\Sigma }},$$
(9)
where $`\mathrm{\Sigma }`$ is the disc surface density just outside the disc gap. The efficiency measures the ratio of the accretion rate on to the planet to the accretion rate that would occur in the disc if the planet was absent. As in Lubow (1999), we obtain accretion rates of order unity. For a 1 MJ planet, $``$ is about a factor of 2 larger than the value obtain in two dimensions. We have also calculated the mass flow rate through $`r=r_\mathrm{p}`$ (i.e. past the planet). In all cases, this rate is much less than the accretion rate onto the planet. For the 1 MJ planet, which opens a well-defined gap in the disc, the flow through $`r=r_\mathrm{p}`$ is less than 2% of the planet’s accretion rate. For all the lower-mass planets, the flow through $`r=r_\mathrm{p}`$ is less than 7% of the planets’ accretion rates. Together, these results imply that mass freely flows into the gap, but that most is captured by the planet. Naturally, the flow past the planet would be expected to increase if the disc’s viscosity were increased.

### 3.4 Migration rates

A planet experiences torques due to its interaction with the disc. Resonant torques likely play an important role (Goldreich & Tremaine 1980). In addition, torques may arise from the gas that flows in the gap, including material within the planet’s Roche lobe. For low-mass planets that undergo Type I migration, resonant torques peak in value at a radial distance of order $`H`$ from $`r_\mathrm{p}`$. Consequently, it is necessary to resolve the gas density structure at distances of order $`H`$ from the planet. Similarly, for planets that open gaps in the disc and undergo Type II migration (i.e. they migrate due to the disc’s viscous evolution), it is necessary to resolve the flow at distances of order $`r_\mathrm{R}H`$ from the planet.

The migration rates due to resonant torques have been the subject of many linear analyses (e.g. Goldreich & Tremaine 1980; Hourigan & Ward 1984; Ward 1986; Korycansky & Pollack 1993; Ward 1997; Tanaka et al. 2002). Ward (1997) considered the motion of the planet relative to the disc material, treated the disc as being two dimensional, and did not consider corotation resonances. More recently, Tanaka et al. calculated the Type I migration rate expected from linear theory taking into account the three-dimensional nature of the disc and corotation resonances.

Tanaka et al.  give the Type I radial migration velocity to be

$$v_\mathrm{I}=f\frac{M_\mathrm{p}}{M_{}}\frac{r_\mathrm{p}^2\mathrm{\Sigma }}{M_{}}\left(\frac{H}{r_\mathrm{p}}\right)^2r_\mathrm{p}\mathrm{\Omega }_\mathrm{p},$$
(10)
where $`\mathrm{\Sigma }`$ and $`H`$ are evaluated at $`r_\mathrm{p}`$, and $`f`$ is a value of order unity that depends on the radial variation of both the disc surface density profile and the scale height. For our disc parameters, we have $`f=3.00`$. The Type II radial migration velocity is simply the viscous radial velocity of the disc

$$v_{\mathrm{II}}=\frac{3\nu }{2r_\mathrm{p}}=\frac{3}{2}\alpha \left(\frac{H}{r_\mathrm{p}}\right)^2r_\mathrm{p}\mathrm{\Omega }_\mathrm{p}.$$
(11)
The migration timescales are then given by $`\tau =r_\mathrm{p}/|v_\mathrm{I}|`$ and $`r_\mathrm{p}/|v_{\mathrm{II}}|`$.

#### 3.4.1 Non-linear results

Both Ward  and Tanaka et al.  assume that the material near the planet (i.e. inside the planet’s Roche lobe) does not exert a net torque on the planet. In this section, we aim to compare their results with those of our simulations. Thus, when evaluating the torque on the planet from our numerical calculations, we only include torques from material at distances greater than some cutoff radius from the planet, $`r_\mathrm{c}`$. Our default is to set the radius of this sphere to be equal to the Roche radius of the planet, $`r_\mathrm{c}=r_\mathrm{R}`$. We discuss the torques inside the Roche lobe in Section 3.4.2.

For each of the 6 planet masses, we calculate the average torque exerted on the planet over the last simulated orbital period (i.e. when the calculations have reached steady accretion). The resulting migration timescales $`\tau `$, are plotted with filled circles in Figure 10. Notice that because the gas flow near the Roche lobe of the 1M planet is affected by the evacuated zones (Section 3.2.2), we have increased $`r_\mathrm{c}`$ to $`2r_\mathrm{R}`$ for this case.

As seen in Figure 10, our numerical results are in excellent agreement with the linear theory of Tanaka et al.  for $`M_\mathrm{p}0.1\text{M}\text{J}`$. Although the linear theory of Tanaka et al. is only formally valid for $`r_\mathrm{R}H`$, the agreement is excellent up to $`r_\mathrm{R}(2/3)H`$. Higher-mass planets migrate more slowly than predicted by equation 10, with a timescale than converges toward the Type II prediction at $`M_\mathrm{p}\text{ }>1\text{M}\text{J}`$. Our migration timescales are about a factor of 3 longer than those reported from the two-dimensional numerical calculations of D’Angelo et al. .

We point out, however, that the corotation resonance is subject to saturation. For our disc parameters saturation may occur for planets of mass greater than about 8 M (Ward 1992). Saturation would result in a reduction of the Type I migration timescale by less than 30%. This level of change does not significantly impact the level of agreement between the simulations and the theory on the scale of Figure 10.

The errorbars in Figure 10 do not represent the uncertainties in our measurements. Rather, they give the migration timescales that are obtained if we change the value of the cut-off radius to $`r_\mathrm{c}=0.5r_\mathrm{R}`$ or $`1.5r_\mathrm{R}`$ ($`r_\mathrm{c}=1r_\mathrm{R}`$ or $`3r_\mathrm{R}`$ for the 1M case). Thus, they indicate the sensitivity of the migration timescale to the material just inside or outside of the Roche lobe.

Our numerical results show a transition from Type I to Type II migration that is qualitatively, but not quantitatively, similar to that predicted by Ward . In Figure 11, we provide an empirical fit to our results as the migration rate $`\tau =r_\mathrm{p}/|v|`$, where

$$v=\frac{v_\mathrm{I}}{1+(M_\mathrm{p}/M_\mathrm{t})^3}+\frac{v_{\mathrm{II}}}{1+(M_\mathrm{t}/M_\mathrm{p})^3}.$$
(12)
Parameter $`M_\mathrm{t}`$ is the transition mass between Type I and Type II migration. We find that a good fit is obtained with $`3/5`$ of the mass for which $`r_\mathrm{R}=H`$, i.e., $`M_\mathrm{t}=1.8M_{}(H/r_\mathrm{p})^3=0.23\text{M}\text{J}`$ (Figure 11). As can be seen from Figures 10 and 11, the transition from Type I to Type II migration involves a much smaller shift in timescales (about a factor of 2), than was suggested by Ward (1997) (about a factor of 10). Note that this shift in timescales is quite sensitive to $`H/r_\mathrm{p}`$ (equations 10 and 11) and increases for thinner discs.

#### 3.4.2 The torque distribution

Linear theory predicts that the strongest resonances occur at a radial distance of order $`H`$ from a planet. But since the resonances have a non-zero width, the strongest torques are exerted over a region within radial distance of order $`H`$ from the planet. In Figure 12, we plot the radial distribution of the torque and the cumulative torque as a function of radius (solid lines). For the low-mass planets that do not begin to open a gap in the disc ($`M_\mathrm{p}10\text{M}\text{}`$), almost the entire torque comes from a region $`r=r_\mathrm{p}\pm 2H`$. For higher-mass planets, the radial extent increases. For 1MJ, the torque comes from $`r=(1\pm 0.25)r_\mathrm{p}`$, which is consistent with a region $`r=r_\mathrm{p}\pm 2(r_\mathrm{R}+H)`$ for the high-mass planets. The resulting cumulative torque distributions are sharply peaked for low-mass planets and broad for the high-mass planets. Also visible in most of the cumulative plots are low-order resonances (e.g. the 2:1 at $`r1.59r_\mathrm{p}`$), although they do not contribute significantly to the total torque. In all cases, the cumulative torque at large radius is negative, indicating inward migration.

As discussed above, in order to compare our results with linear theory, we have specifically excluded torques from material inside the Roche lobe of the planet ($`r_\mathrm{c}=r_\mathrm{R}`$). In any case, as demonstrated by the gas flows in Figure 5, we only resolve the flow inside the Roche lobe for $`M_\mathrm{p}0.1\text{M}\text{J}`$. The question arises as to whether torques from material inside the Roche lobe can affect the migration rate of the planet.

In Figure 12, the dashed lines give the torque exerted on the planet including all zones on the numerical grid. For $`M_\mathrm{p}0.1\text{M}\text{J}`$ (for which the flow inside the Roche lobe is resolved), the net torque inside the Roche lobe is positive. In fact, for the 1MJ case, the total torque is slightly positive (i.e., outward migration). The problem with including torques from deep inside the Roche lobe is that near the planet the torque per unit mass becomes very large. Therefore, small departures from axisymmetry (either real or due to finite numerical resolution) can produce a large net torque that is comparable to the net torque from outside the Roche lobe. As shown by the errorbars in Figures 10 and 11, including torques into $`r_\mathrm{c}=0.5r_\mathrm{R}`$ only results in changes at the 20% level (for $`M_\mathrm{p}10\text{M}\text{}`$). Here the flow is still reasonably well resolved. Going deeper than this, however, gives very unreliable results.

The same problem has been encountered in other studies. Lubow et al.  found that the contributions to the torque from radii near the planet were in the opposite sense to the torque from further out in the disc (i.e. radii outside but near the planet gave a positive net torque and radii inside but near the planet contributed a negative torque). D’Angelo et al.  found that taking into account torques near the planet usually resulted in the migration switching from inward to outward. They neglected the contribution from zones within approximately $`0.2r_\mathrm{R}`$ when calculating the net torque. The problem is worse in these two dimensional calculations than in three dimensions because of the strong spiral shocks within the circumplanetary disc – any asymmetry due to numerical resolution results in a large net contribution.

A related issue is that the material that is accreted by the planet and removed from the calculation carries some angular momentum. If this angular momentum were put back into the orbital angular momentum of the planet it would reduce the migration rate of the planet. The magnitude of this effect depends on the radius of the region from which the gas is removed. In most of our calculations, the amount of angular momentum that is removed is low. If it were all added to the orbital angular momentum of the planet it would reduce the migration rates by approximately 2% in the cases of the 1 and 0.3 MJ planets, by approximately 10% percent for the 0.1, 0.03 and 0.01 MJ planets, and by approximately 50% percent in the case of the 1 Earth-mass planet. The large effect for the Earth-mass planet is partially due to its accretion rate being overestimated (Section 3.2.2).

To accurately evaluate the torques from inside the Roche lobe will require higher resolution, three-dimensional calculations. It may also require a realistic equation of state, since low-mass planets may be unable to accrete gas at the rate the disc can provide it because the thermal energy of the gas cannot be radiated away quickly enough (see Section 4). For the present, we have presented an accurate determination of the torque from outside the Roche lobe and find excellent agreement with the linear theory of Tanaka et al. .

### 3.5 Further tests of Type I migration

The discs in the above calculations had unperturbed scaleheights $`H/r=0.05`$ and surface density profiles $`\mathrm{\Sigma }r^{1/2}`$. Tanaka et al. (2002) predicted how the migration timescales should vary with the disc’s thickness and surface density profile. To further test their predictions, we performed calculations of 10 M planets in discs with different thicknesses and surface density profiles.

We performed a calculation with $`H/r=0.10`$. According to equation 10, the planet should migrate 4 times slower than with $`H/r=0.05`$. We find that the planet migration timescale increases by a factor of about $`4.9\pm 1.4`$, where the estimated error comes from using values of $`r_\mathrm{c}`$ ranging between $`r_\mathrm{c}=0.5r_\mathrm{R}`$ and $`1.5r_\mathrm{R}`$. We also note that the torque density and cumulative torque distributions are similar to those in the $`M_\mathrm{p}=0.03\text{M}\text{J}`$ panels of Figure 12, except they are spread over twice the radial range around the planet. Thus, the torque is generated almost entirely within the radial range $`r=r_\mathrm{p}\pm 2H`$, as discussed in the previous section.

We performed a calculation with $`\mathrm{\Sigma }r^{3/2}`$. According to Tanaka et al. , the planet should migrate 36% faster than with $`\mathrm{\Sigma }r^{1/2}`$. We find that the planet migration timescale decreases by a factor of about $`1.8\pm 0.2`$. The distribution of the torque is similar to the $`M_\mathrm{p}=0.03\text{M}\text{J}`$ panels of Figure 12.

Given the difficulties in measuring the torques, these results are in reasonable agreement with the predictions of Tanaka et al. .

## 4 Discussion

The accretion and migration rates from the previous sections can be used to investigate the timescale for the gas accretion phase of giant planet formation. In the core-accretion model of giant planet formation, giant planets are thought to form through gas accretion on to a solid core of about 10 M (e.g., Mizuno 1980; Hayashi, Nakazawa & Nakagawa 1985). Runaway gas accretion is thought to occur at higher masses, greater than about 50 M (Pollack et al. 1996). The overall formation timescale may be dominated by the phase between the onset of gas accretion and the beginning of runaway gas accretion. During this stage, the accretion on to the planet is limited by the energy losses in the contracting envelope. This model is subject to several major assumptions, and there are issues with the apparently small observationally inferred core mass of Jupiter (Pollack et al. 1996; Wuchterl et al. 2000). In any case, the accretion rates during this intermediate stage are typically much smaller than the rates we plotted in Figure 3.3. Consequently, the accretion rates we obtain are valid, within the framework of the current calculations of the core-accretion model, only during the runaway gas accretion phase.

As a protoplanet accretes, it will migrate through the disc. Its rate of migration is given by Figure 11. Thus, we can investigate the timescale and result of giant planet formation by starting with a core and integrating its mass and orbital radius forward in time. The results of such integrations are given in Figure 13. The evolution of the protoplanet could be stopped at any point along the curves by dispersal of the disc. We plot results for cores beginning at three different orbital radii and for three different disc surface densities (i.e. three different disc masses). We have assumed that the disc surface density varies as $`\mathrm{\Sigma }r^{1/2}`$, but the results are very insensitive to the index; $`\mathrm{\Sigma }r^{3/2}`$ gives almost identical results.

For convenience, we assume the initial core mass is 1 M and ignore the possibility that the planet may not be able to accept mass at the rates we obtain. Notice, however, the results we obtain are insensitive to the value of the threshold mass for rapid gas accretion, as long as it is less than approximately 50 M. The time required for the gas accretion to produce a giant planet is almost independent of the initial core mass, because the mass doubling timescale for a low-mass core is much shorter than that for a high-mass planet (Figure 9, right panel). Thus, only a small fraction of the total time (Figure 13) is spent while it is a low-mass core. Similarly, the orbital radius of the planet is independent of the initial core mass because the migration timescale of a low mass object is much longer than its mass doubling timescale. Thus, the timescale of runaway gas accretion phase and the final mass and orbital radii of a giant planet are essentially independent of whether one assumes that runaway gas accretion begins at 1 M or 50 M.

We find the time required for a core to accrete to 1 MJ is only $`t_{\mathrm{acc}}=2.0\times 10^32.0\times 10^4`$ years. None of the objects migrates significantly during this period, and the final orbital radius of the 1 MJ planet simply depends on the initial orbital radius of the core.

Much more of a problem is the formation of very massive planets. To form a 4 MJ planet requires an order of magnitude longer ($`t_{\mathrm{acc}}=2.0\times 10^42.0\times 10^5`$ years). Cores that begin rapid gas accretion at 2 AU migrate into the star before they reach 4 MJ, unless the disc is massive ($`\mathrm{\Sigma }\text{ }>200`$ g cm-2 at 5.2 AU). Cores that begin at larger radii and survive all migrate significantly during this period. The formation of giant planet with $`M_\mathrm{p}>7\text{M}\text{J}`$ is only possible for cores that begin at large distances in massive discs. To form a giant planet of 10 MJ requires either a very massive circumstellar disc, a large initial orbital radius, or both.

The above rates represent the fastest possible accretion and the least possible migration. Even then, the formation of planets with masses $`M_\mathrm{p}10`$ MJ is difficult or impossible. Thus, there is a natural limit to the mass of a giant planet of $`10`$ MJ, as suggested in Lubow et al. . This is in good agreement with the mass distribution of extra-solar planetary systems (e.g. Basri & Marcy 1997; Mayor, Queloz, & Udry 1998; Mazeh, Goldberg, & Latham 1998; Halbwachs et al. 2000; Zucker & Mazeh 2001; Jorissen, Mayor, & Udry 2001; Tabachnik & Tremaine 2002). However, it is worth pointing out that our results apply only to planets with circular orbits. Planets with highly eccentric orbits may accrete significantly even for masses $`M_\mathrm{p}>10`$ MJ. Thus, the eccentric orbits of many observed extrasolar planets may assist them in growing to large masses.

## 5 Conclusions

We simulated the three-dimensional interaction of a young planet, ranging in mass from 1 M to 1 MJ, with a gaseous disc by means of the ZEUS hydrodynamics code. The disc was vertically isothermal with an unperturbed disc thickness ratio was $`H/r=0.05`$ and disc turbulent viscosity parameter was $`\alpha =4\times 10^3`$. We have analyzed the flow patterns, the accretion rates, and migration rates. Each planet was assumed to remain in a circular orbit and accrete gas without expansion on the scale of its Roche lobe. To incorporate the latter assumption into the simulations, we removed mass from the grid zones that immediately surround the planet at each timestep.

Only planets with masses $`M_\mathrm{p}\text{ }>0.1`$ MJ produced significant perturbations in the disc’s surface density; $`0.03`$ MJ (10 M) is insufficient (see Figure 2). The flow near the planet is fully three-dimensional (see Figure 7). The flow at the disc midplane generally involves material that passes by the planet, two gas streams that penetrate the Roche lobe and supply material to a circumplanetary disc, and material on horseshoe orbits (Figure 4). These features are similar to those found previous two-dimensional studies (e.g., Bryden et al. 1999; Kley 1999; Lubow et al. 1999), but with some important differences.

The shocks in the circumplanetary flow are much weaker in the three-dimensional case than in the two-dimensional case (see Figures 5 and 8). The circumplanetary disc in three dimensions is likely to behave more like a standard accretion disc, rather than being subject to strong shock-driven accretion, as is the case in two-dimensions.

The gas streams at midplane are narrower in the three-dimensional case. However, the overall efficiency of accretion in three dimensions is still high. The accretion rate peaks at approximately 0.1 MJ, but a 1 MJ planet still accretes mass at a rate greater than the usual local viscous rate. This result suggests that the lower accretion occurring at the midplane in the three dimensional calculations is compensated by accretion occurring from above the midplane. For small mass planets, the accretion occurs with a cross-section whose length scale is of order the size of the Roche lobe (see Figures 4 and 9). This leads to the accretion rate increasing in proportion to the planet’s mass (equation 8).

We investigated the disc torques on planets and the resulting migration timescales. The migration timescales obtained from torques exerted outside the planet’s Roche lobe are in excellent agreement with recent linear theory (Tanaka et al. 2002), as seen in Figure 10. The transition from Type I (non-gap) migration to Type II (gap) migration occurs at $`M_\mathrm{p}\text{M}\text{J}/4`$. The transition is smooth with only about a factor of 2 difference in rates. The torque outside the Roche lobe is concentrated in a region of distance of order $`H`$ from the planet, as expected from linear theory.

It is not yet possible to determine the torques from within the Roche lobe of a planet. The problem is that the torque per unit mass increases as the distance from the planet decreases. Consequently, small numerical fluctuations in the density close to the planet give rise to numerical noise in the torque (see Figure 12). On the other hand, it is not clear that strong torques can arise close to the planet, since the outer portions of the circumplanetary disc in three dimensions appear to be quite smooth.

We considered the orbital and mass evolution of a planet, based on the orbital migration and mass accretion rates we obtained (Figure 13). For low-mass protoplanets, less than about 50 M, the accretion rate may be limited by the response of the planet (see discussion in Section 4). In the absence of this limit, a planet could gain considerable mass with little migration. Starting with a higher mass core of say 50 M, higher mass planets, up to about 2 MJ, can also grow with little migration. Masses as high as 10 MJ do not appear possible (as was also suggested in Lubow et al. 1999), because the planet will migrate inward over a large distance before accreting enough material.

## Acknowledgments

The computations reported here were performed using the U.K. Astrophysical Fluids Facility (UKAFF) and the GRAND computer. MRB is grateful for support from the Space Telescope Science Institute’s Visitor’s Programme and the Institute of Astronomy’s Visitor’s Programme. SHL acknowledges support from NASA grants NAG5-4310 and NAG5-10732.


text:
# Jet quenching and event-wise mean-𝑝_𝑡 fluctuations in Au-Au collisions at √𝑠_{𝑁𝑁}= 200 GeV in Hijing-1.37

## I Introduction

QCD lattice gauge calculations  suggest that a phase transition from hadronic nuclear matter to a color-deconfined medium or quark-gluon plasma (QGP) may occur in relativistic heavy ion collisions at RHIC. Measurements of multiplicity and transverse-momentum fluctuations have been proposed to establish the existence and study the properties of a QGP. It was suggested that in case of a second-order phase transition fluctuations may be reduced dramatically , whereas for a first-order phase transition large fluctuations may occur due to localized QGP droplet formation , or more gradual density or temperature fluctuations may occur. The past two years have seen significant progress in the number and quality of preliminary fluctuation and correlation measurements at both SPS and RHIC. Further efforts are now required to draw quantitative conclusions about the phase transition and the physical characteristics of the QGP.

In this Letter we consider the possibility that a major source of event-wise mean-$`p_t`$ or $`p_t`$ fluctuations in relativistic heavy ion collisions is jet/minijet production, and that these fluctuations can therefore be used to study the early collision process and the properties of the QGP, in a manner complementary to studies of inclusive event properties such as transverse momentum distributions .

According to  minijets (semi-hard parton scattering with momentum transfers of a few GeV/c) should be copiously produced in the initial stage of ultra-relativistic heavy ion collisions at RHIC energies. A QGP formed in the early stage of such collisions could act as a dissipative medium, causing jets/minijets to lose energy through induced gluon radiation , a process conventionally called jet quenching in the case of higher-$`p_t`$ partons. The properties of the dissipative medium should determine the specific energy loss of jets and minijets . We explore the role of fluctuation measurements in understanding the dissipation properties of a color-deconfined medium, using Hijing-1.37  to study jet-quenching and jet/minijet contributions to $`p_t`$ fluctuations.

## II Hijing Event Types

Three sets of Hijing-1.37 events were generated for this study: Hijing default with jet quenching on, with quenching switched off and with jet/minijet production switched off. Because jet quenching effects should be most significant in central collisions this study was applied to Hijing-1.37 events with impact parameter less than 3 fm, equivalent to the top 5% most central collisions.

Fluctuation analysis was applied to all generated charged particles subject to the conditions that the pseudorapidity and transverse momentum satisfied $`|\eta |<1`$ and 0.15 GeV/c $`<p_t<`$2 GeV/c. The total number of events, the mean $`\widehat{p}_t`$ and variance $`\sigma _{\widehat{p}_t}`$ of the charged-particle inclusive $`p_t`$ distributions, and the mean $`\overline{n}`$ and variance $`\sigma _n`$ of the charged-particle multiplicity distributions for the three types of central Hijing events are listed in Table I. We analyze $`p_t`$ fluctuations for these three event classes using graphical and numerical methods.

## III Graphical Analysis of Fluctuations

Graphical analysis of $`p_t`$ fluctuations directly compares the data distribution with a reference. The gamma distribution models statistically-independent particle emission . Any difference between data and reference would imply that emitted particles are correlated and/or event-wise fluctuations in collision dynamics (such as minijets) are present. Effects of fluctuating event multiplicity in this method require folding the gamma distribution with a negative binomial distribution . The reference distribution is formulated as:

$$f(p_t)=\underset{n_{min}}{\overset{n_{max}}{}}f_{NBD}(n,1/k,\overline{n})f_\mathrm{\Gamma }(p_t,np,nb)$$
(1)
where $`f_{NBD}(n,1/k,\overline{n})`$ and $`f_\mathrm{\Gamma }(p_t,np,nb)`$ are the negative binomial distribution and n-fold gamma distribution respectively

$`f_{NBD}(n,1/k,\overline{n})`$ $`=`$ $`{\displaystyle \frac{(n+k1)!}{n!(k1)!}}{\displaystyle \frac{(\frac{\overline{n}}{k})^n}{(1+\frac{\overline{n}}{k})^{n+k}}}`$ (2)
$`f_\mathrm{\Gamma }(p_t,np,nb)`$ $`=`$ $`{\displaystyle \frac{nb}{\mathrm{\Gamma }(np)}}(nbp_t)^{np1}\mathrm{exp}^{nbp_t}`$ (3)
The parameters in Eq. (2) are defined as $`p=\widehat{p}_t^2/\sigma _{\widehat{p}_t}^2`$, $`b=\widehat{p}_t/\sigma _{\widehat{p}_t}^2`$ and $`\overline{n}/k=\sigma _n^2/\overline{n}1`$, where $`\widehat{p}_t`$ and $`\sigma _{\widehat{p}_t}^2`$ are the mean and variance respectively of the inclusive $`p_t`$ distribution. This model is related to the central limit theorem result that the cumulants of a distribution of independent $`n`$-sample means from a fixed parent distribution are related to the cumulants of the parent by inverse powers of the sample number. The gamma distribution cumulants satisfy this relation, and the inclusive $`p_t`$ distribution is approximately described by a gamma distribution with $`n=2`$. The relations above then follow. The graphical comparison is important as a qualitative indicator of excess variance and of significant deviations from central limit behavior in cumulants higher than the second.

### A Basic histograms

Distributions of $`p_t`$ for each of three Hijing event types are shown in Fig. 1 as histograms together with reference distributions based on Eq. (1) plotted as curves. The reference distributions are normalized to have the same peak values as the corresponding data distributions to emphasize the width comparison. It can be seen from Fig. 1 that the Hijing data histograms deviate substantially from the reference distributions, indicating that nonstatistical fluctuations are significant. Fig. 1 indicates that excess fluctuations in events with jets switched off are substantially less than those for the other two event classes.

Comparison of the top and bottom panels of Fig. 1 indicates qualitatively that jet production plays a significant role in $`p_t`$ fluctuations in the Hijing model. To investigate the effect of jet quenching we can compare nonstatistical fluctuations in top and middle panels of Fig. 1. However, with these histograms it is not possible to detect visually a significant difference in the magnitude of nonstatistical $`p_t`$ fluctuations. We require a more differential approach.

### B Difference histograms

The graphical analysis can be made more sensitive by plotting the difference between data histogram and reference distribution, facilitating quantitative comparisons between Monte Carlo configurations and references. The normalized event-number difference is defined as

$$\delta N/\sqrt{N}\frac{N_{data}N_{reference}}{\sqrt{N_{data}}}$$
(4)
The independent random variable $`\sqrt{\overline{n}}(p_t\widehat{p}_t)/\sigma _{\widehat{p}_t}`$ normalizes the distribution width relative to the inclusive width according to a central limit expectation. The result is a universal plotting format; the horizontal and vertical axis units are standard deviations.

Figure 2 shows results for the three Hijing-1.37 configurations. Significant deviations from zero reveal nonstatistical fluctuations in all three event types. The dramatic increase from jet-off events (lowest curve) to events with jets on is consistent with the qualitative difference in Fig. 1. A significant difference with and without jet quenching (upper two curves) is also now apparent. This difference indicates that according to the Hijing model transverse momentum fluctuations are sensitive to jet quenching, which may in turn reflect the properties of a dense colored medium or QGP formed in relativistic heavy ion collisions. We conclude that while the original graphical method is an indicator for fluctuations associated with jet production, the difference histograms are more sensitive to the effects of jet quenching.

## IV Numerical Analysis of Fluctuations

Numerical analysis of nonstatistical fluctuations requires comparison of a variance from data with the variance of a reference representing only statistical fluctuations. The linear variance comparison measure $`\mathrm{\Delta }\sigma _{p_t}`$  was developed to isolate pure $`p_t`$ fluctuations from effects of multiplicity fluctuations which may appear as a systematic bias in some statistical measures.

Several quantities have been proposed to measure event-wise $`p_t`$ fluctuations in relativistic heavy ion collisions: $`f_{p_t}`$, $`F_{p_t}`$, $`\mathrm{\Phi }_{p_t}`$, and $`\sigma _{p_t,dynamical}^2`$. A detailed study of $`\mathrm{\Phi }_{p_t}`$ is given in . In the present analysis the closely-related $`\mathrm{\Delta }\sigma _{p_t}`$ is applied to event-wise $`p_t`$ fluctuations in Hijing-1.37 events. This difference factor for $`p_t`$ fluctuations is defined as

$$\mathrm{\Delta }\sigma _{p_t}\left\{\overline{n(p_t\widehat{p}_t)^2}\sigma _{\widehat{p}_t}^2\right\}/2\sigma _{\widehat{p}_t}$$
(5)
where the overline means average over all events in a data set. Values of $`\mathrm{\Delta }\sigma _{p_t}`$ for the three Hijing event types are included in Table I. Consistent with Fig. 1 the results show very significant contributions from jets, and they also show that fluctuations are significantly reduced when jet quenching is applied, consistent with Fig. 2.

As expected, numerical analysis of $`p_t`$ fluctuations reveals more precisely than graphical analysis that jet/minijet production in the Hijing model contributes significantly to nonstatistical event-wise $`p_t`$ fluctuations. The increase in fluctuations for jets on compared to jets off is a factor 3-5. The effect of jet quenching is to reduced $`\mathrm{\Delta }\sigma _{p_t}`$ by about 30% relative to no quenching, 20 times the statistical error. These results suggest that properties of the QGP or dissipative medium formed in the initial stage of heavy ion collisions may be studied precisely with numerical analysis of event-wise $`p_t`$ fluctuations.

## V Centrality Dependence

However, we observe an excess of binary-collision scaling in Hijing compared to RHIC data which affects the interpretation of this fluctuation analysis. Hijing particle production increases with number of binary collisions

(centrality) at twice the rate observed in RHIC Au-Au collisions, as indicated by the multiplicity entries in Table I. The jets-off multiplicity reflects soft processes, whereas for the same number of participant pairs the jets-on multiplicities, which include hard-component scaling with binary collisions, represent a 2.5-3-fold increase over participant scaling for central events. For RHIC data we observe a corresponding 1.5-fold increase . The fraction of particle production from hard processes in Hijing events is thus about twice that observed for RHIC events. We find that $`p_t`$ fluctuations in Hijing (jets on) are dominated by the hard component, by a factor 5. It may be advisable therefore to reduce by a factor of two the $`p_t`$ fluctuations attributed by this study to jets/minijets in Hijing in comparisons with RHIC collisions.

## VI Hadronic Rescattering

Hadron rescattering between chemical and kinetic de-coupling might attenuate nonstatistical fluctuations generated in the early stage of the collision, which could then be confused with effects of dissipation in the prehadronic medium. Recent Monte Carlo calculations suggest that hadronic rescattering may be negligible : the time for rescattering may be so brief  that fluctuation attenuation is not significant. Nonstatistical $`p_t`$ fluctuations due to jets/minijets would thus survive hadronic rescattering unscathed. This is an important point for interpreting attenuation of fluctuations observed in the final state. Reduced $`p_t`$ fluctuations as reported in this Letter would then be an unambiguous indication of prehadronic dissipation (minijet/jet quenching).

## VII Conclusions

In summary, event-wise $`p_t`$ fluctuations in central Au-Au collisions at $`\sqrt{s_{NN}}=\mathrm{\hspace{0.25em}200}`$ GeV have been studied with Hijing-1.37 using graphical and numerical methods. The study shows that jet/minijet production in the initial stage of heavy ion collisions may contribute substantially to nonstatistical $`p_t`$ fluctuations. We have also shown that in the Hijing model $`p_t`$ fluctuations are sensitive to jet quenching, thus suggesting that fluctuation analysis could serve as a probe of dissipation properties of a dense colored medium formed in the initial stages of ultrarelativistic heavy ion collisions. Studies of fluctuation centrality dependence in models and data may help clarify questions on hadronic dissipation and Hijing binary-collisions scaling raised by this study.

We appreciate helpful discussions with J.G. Reid (University of Washington). This work was supported in part by USDOE contract DE-FG03-97ER41020.


text:
# VLT observations of the solitary millisecond pulsar PSR J2124-3358

## Introduction

The about 1500 radio pulsars known to date are interpreted as rapidly spinning and strongly magnetized neutron stars which are radiating at the expense of their rotational energy. Millisecond pulsars form a separate group among the rotation-powered pulsars. They are distinguished by their small spin periods ($`20`$ ms) and their high spin stability ($`\dot{P}10^{18}10^{21}`$ s s-1), with corresponding spin-down ages $`P/2\dot{P}`$ of typically $`10^910^{10}`$ years. Only $`20\%`$ of the galactic millisecond pulsars are solitary, including PSR B1257+12 which is in a planetary system. The rest are in binaries, usually with a low-mass white dwarf companion. Optical observations of binary millisecond pulsars have so far only allowed to detect the companion star (see, e.g., Danziger et al. 1995; Lundgreen et al. 1996) and to constrain the global evolution of the binary systems. In no case it was possible to detect the fainter optical emission from the neutron star which is buried by the brighter companion. Optical observations of solitary millisecond pulsars thus are the only chance to explore the optical emission properties of these peculiar objects.
The optical emission of solitary millisecond pulsars can be ascribed either to magnetospheric or thermal radiation from the neutron star surface. However, unless some reheating mechanism is at work, optical thermal radiation is expected to be virtually undetectable for old ($`10^9`$ yrs), cooled off, neutron stars and magnetospheric radiation appears the most likely emission process. From the model developed by Pacini (1971), the magnetospheric emission is expected to depend on the pulsar’s spin parameters according to the relation $`L_{opt}B^4P^{10}`$. Although valid for young pulsars such a dependance, as well as any breakdown with the pulsar’s age, is still unclear for older neutron stars. Detecting optical emission from millisecond pulsars, i.e. short period and low $`\dot{P}`$ objects, would thus be of paramount importance to define a general framework.
As for classical, i.e., isolated non-millisecond pulsars, a detection at soft X-ray energies represents an important indicator that the object is likely detectable at optical wavelength. Currently, soft X-ray emission have been discovered from a number of galactic solitary millisecond pulsars but only for a few of them X-ray pulses have been detected (see Becker & Trümper 1999; Becker & Aschenbach 2002).

## PSR J2124$``$3358

One of the most interesting objects is PSR J2124$``$3358. The pulsar was identified in radio during the Parkes 436 MHz southern survey (Bailes et al. 1997). Its measured $`\dot{P}`$ of $`1.09610^{20}`$ s s-1 gives an age of 7.2 billion years, a magnetic field of $`2.310^8`$ G and a rotational energy loss of $`3.510^{33}`$ ergs s-1. PSR J2124$``$3358 was discovered in soft X-ray energies during a recent study of galactic isolated millisecond pulsars using the ROSAT HRI (Becker & Trümper 1999). Its X-ray emission was found to be pulsed at the radio period, with a sharp double-peaked pulse profile. The shape of the X-ray pulse profile was considered a strong evidence that the X-ray emission from this millisecond pulsar is dominated by non-thermal emission (cf. Becker & Trümper 1999). Indeed, the X-ray lightcurve exhibits a double peak structure remarkably similar to the pulse profile observed in the radio domain at 436 MHz (Bailes et al. 1997; Becker & Trümper 1999). Figure 1 shows a comparison between the PSR J2124$``$3358 pulse profiles as oberved in the X-ray and radio bands. The interpretation in terms of magnetospheric emission has been later supported by ASCA observations which showed that the pulsar spectrum in the 0.5-10 keV interval is compatible with a power law (Sakurai et al. 2001).
PSR J2124$``$3358 is therefore one of the most promising millisecond pulsars for optical studies as it turns out to be one of the brightest solitary millisecond pulsars detected at X-ray energies. Its close distance of about 270 pc (Cordes & Lazio 2002) together with its small Hydrogen column absorption of $`2510^{20}\text{cm}^2`$ (Becker & Trümper 1999) also make this object a natural target for subsequent optical observations.

## Observations

The first optical observations of the PSR J2124$``$3358 field were performed in July 1998, shortly after the detection of the pulsar in X-rays. Observations were performed with the ESO NTT using the ESO Multi Mode Instrument (EMMI) in its Blue Arm configuration with a pixel size of 0.37” (field of view $`6.2^{}\times 6.2^{}`$). A total of six exposures (1200s each) were acquired through a $`B`$ filter. Single images were reduced using standard MIDAS recipes and combined using a median filter algorithm. Photometric calibration was computed through short exposures of Landolt’s fields. Unfortunately, since all the exposures of the PSR J2124$``$3358 field were affected by bad atmospheric conditions and by a seeing constantly above 1.5” it was not possible to reach a detection limit fainter than $`B25`$ (Becker & Mignani, unpublished).
The field of PSR J2124$``$3358 was observed again between August and September 2001 using the ESO VLT/Antu telescope operated in Service Mode from the Paranal Observatory (Becker et al. 2002). Imaging was performed using the FOcal Reducer and Spectrograph 1 (FORS1) camera operated in its standard resolution mode, with a pixel size of 0.2” and a corresponding field of view of $`6.8^{}\times 6.8^{}`$. Images were taken through the $`U`$ ($`\lambda =3660`$Å, $`\mathrm{\Delta }\lambda =360`$Å), $`B`$ ($`\lambda =4290`$Å, $`\mathrm{\Delta }\lambda =880`$Å) and $`V`$ ($`\lambda =5540`$Å, $`\mathrm{\Delta }\lambda =1110.5`$Å) passbands to allow for a significant spectral coverage. For each filter, sequences of 12 repeated exposures were obtained for a total integration time of 9360 s in $`U`$ and 6000 s in $`B`$ and $`V`$. The detailed summary of the observations is listed in Table 1. To allow for relative adjustments in the flux calibration, one exposure of each sequence was obtained under photometric conditions. The pulsar field was always observed close to the minimum airmass with average values of 1.10 ($`U`$), 1.09 ($`B`$) and 1.03 ($`V`$). Seeing conditions across all the exposures varied between $`0.6\mathrm{"}`$ and $`1.0\mathrm{"}`$.
Standard reduction steps were applied through the FORS1 image reduction pipeline. For each night, master bias and sky flats were used for debiassing and flatfielding. Flux calibration was computed using the extinction and color-corrected photometric zero-points routinely computed for each night. For each passband, sequences of repeated exposures were finally co-added by applying a median filter algorithm to reject cosmic ray hits.

## Results

Since the first hint to claim the optical identification of the pulsar comes from a positional coincidence with a potential counterpart, an accurate astrometry of the field is a key issue. In order to register the very accurate (down to few milliarcseconds) pulsar position (Gaensler, Jones and Stappers 2002) on the FORS1 images, we have recomputed the FORS1 astrometry using as a reference the positions of stars selected from the Guide Star Catalogue II, which have an intrinsic absolute astrometric accuracy of $`0.35^{\prime \prime }`$ per coordinate. After identifying a number of well suited GSC-II objects in the co-added $`V`$-band image, the pixel-to-sky coordinate transformation has been computed using the program ASTROM111http://star-www.rl.ac.uk/Software/software.htm. The final precision on the astrometric fit was $`0.09^{\prime \prime }`$ in both RA and Dec. To this value we then added an error of $`0.17^{\prime \prime }`$ due to the propagation of the intrinsic absolute errors on the GSC-II coordinates.

Since the error on the pulsar position is negligible, the overall uncertainty in the registration of the radio coordinates on the optical frames is dominated by the accuracy of the astrometric catalogue used as a reference and by the rms of the astrometric fit. Taking both factors into account, we finally ended up with an uncertainty of $`0.2^{\prime \prime }`$ in both RA and Dec. Any additional uncertainty on the pulsar position due to its known proper motion of $`\mu =52.6`$ mas yr-1 (Cordes, Jones and Stappers 2002) is certainly within our global error budget. We note that the GSC-II astrometry was tied to the extragalactic radio source frame (ICRF), thus our position should not be significantly affected by any systematic offset between the radio and optical reference systems. Our astrometry is shown on the co-added $`UBV`$-band images (Figure 2). No object is detected at the pulsar position down to $`3\sigma `$ limiting magnitudes of $`U26`$, $`B27.7`$ and $`V27.8`$. Thus, we assume these values as upper limits on the optical flux of the pulsar. By scaling for the assumed distance of 270 pc (Cordes & Lazio 2002), our upper limits imply an optical luminosity $`L_{opt}2.510^{26}`$ ergs s-1, which makes PSR J2124$``$3358 intrinsically fainter than the old ($`10^610^7`$ yrs) ordinary pulsars PSR B1929+10 and PSR B0950+08 (Pavlov et al. 1996).

## Conclusions

We have reported on the first, deep, optical observations of the field of the millisecond pulsars PSR J2124$``$3358 performed with the ESO VLT. The pulsar is undetected down to a flux limit of $`V27.8`$, the deepest obtained so far for an object of this class. Unfortunately, this result adds to the list of non-detections of millisecond pulsars deriving from (shallower) observations recently performed with the VLT, namely PSR J1024$``$0719, PSR J1744$``$1134 (Sutaria et al. 2002) and PSR J0030+0451 (Koptsevich et al. 2002). Future observations of these objects in the near-UV with the HST will probably offer higher chanches of detection.