text:
# Local Inhomogeneity in Asymmetric Simple Exclusion Processes with Extended Objects

## 1 Introduction

In recent years, asymmetric simple exclusion processes (ASEP) have become a subject of increasing scientific interest because of their crucial role in the investigations of numerous dynamic phenomena in chemistry, physics and biology . They are important for understanding mechanisms of biopolymerization , reptation polymer dynamics , diffusion through membrane channels , traffic problems , and protein synthesis . ASEP are non-equilibrium one-dimensional lattice models where particles interact only through hard-core exclusion and move preferentially in one direction.

Although dynamic rules of asymmetric simple exclusion processes are very simple, they show a very rich dynamic phase behavior. One of the most striking features of ASEP with open boundaries is the occurrence of non-equilibrium phase transitions between stationary states that have no analogs in equilibrium systems . The physics of these phase transitions can be explained by utilizing a phenomenological domain-wall theory . According to this theory, the entrance, the bulk, and the exit of the system define their own stationary domains with specific uniform densities and currents. Domain walls exist in the border region between different domains, and the dynamics of these domain walls, which depends on the system parameters, in the stationary-state limit will determine the dominant phase. In the maximal-current phase, the steady-state density profile and current are due to the bulk dynamics, while the low-density (high-density) phase is enforced by the entrance (exit) rate.

Recently, a new class of ASEP with particles occupying more than one lattice site has been investigated using various mean-field and continuum approaches . These models provide a more realistic description of many biological processes such as RNA translation, vesicle locomotion along filaments, and proteins sliding along DNA. For example, in RNA translation a ribosome typically covers 10-12 codon sites but moves only one codon at a time . However, one more important dynamic feature, the non-uniformity of hopping rates, should also be considered for ASEP with extended objects in order to give a more realistic description of one-dimensional biological transport phenomena. This issue has been only briefly considered by using crude approximation methods and Monte Carlo simulations . More thorough theoretical investigations of effects of inhomogeneity are needed. The main goal of this paper is to develop a theoretical and computational description of inhomogeneous ASEP with particles of arbitrary size. We consider here the totally asymmetric simple exclusion processes (TASEP), where particles can only move in one direction, although our theoretical arguments and our method can be extended to more general exclusion models.

We present a theoretical study of the simplest inhomogeneous asymmetric exclusion model with particles of arbitrary size. In this model, the only nonuniform hopping rate is at the site which is in the middle of the lattice. We investigate this system by using a simple approximate theory, based on the domain-wall approach, and extensive computer simulations. The paper is organized as follows. In section 2 we outline our model and known results for homogeneous ASEP with extended objects, and we present calculations from a simple approximate theory and domain- wall arguments. Monte Carlo simulations are discussed in section 3. We summarize and conclude in section 4.

## 2 TASEP with extended objects

### 2.1 Model

We consider identical particles moving on a one-dimensional lattice consisting of $`N`$ sites, as shown in Fig. 1. Each particle may cover $`l1`$ sites. For convenience, we associate the position of each particle with the position of its left edge. In the bulk of the system, a particle located at site $`i`$ will move to the next site $`i+1`$ with the rate 1, given that site $`i+l`$ is empty. Furthermore, when the particle is at the special site $`k`$, far away from the boundaries, the particle jumps forward to site $`k+1`$ with the rate $`q`$ (provided that site $`k+l`$ is empty). Thus particles travel from the left to the right with uniform rates except at the special site $`k`$, which is the place of local inhomogeneity. Particles also can enter the system from the left with rate $`\alpha `$ if the first $`l`$ sites are empty. When a particle reaches site $`N`$, it can exit with rate $`\beta `$. We consider $`N`$ to be very large, and in this thermodynamic limit the exact details of entrance and exit dynamics are not very important . As in , we define the density of the particles in the system as the coverage density, namely, for $`M`$ particles on the lattice the density is given by $`\rho =Ml/N`$. Since we discuss only the totally asymmetric simple exclusion model with extended objects, particles can move in only one direction, from left to right as shown in Fig. 1.

### 2.2 Results for homogeneous TASEP with extended objects

Analytical mean-field calculations and domain-wall arguments, supported by extensive Monte Carlo simulations, indicate that homogeneous TASEP with extended objects has the same three-phase diagram as homogeneous standard TASEP with $`l=1`$ .

For $`\alpha \frac{1}{\sqrt{l}+1}`$ and $`\beta \frac{1}{\sqrt{l}+1}`$, the dynamics of the system is determined in the bulk, and we have a maximal-current phase with the stationary current and bulk density given by

$$J=\frac{1}{(\sqrt{l}+1)^2},\rho =\frac{\sqrt{l}}{\sqrt{l}+1}.$$
(1)
When the particle entry is the rate-limiting step, which is realized for $`\alpha <\frac{1}{\sqrt{l}+1}`$ and $`\alpha <\beta `$, the system is in a low-density stationary phase with the following current and bulk densities:

$$J=\frac{\alpha (1\alpha )}{1+(l1)\alpha },\rho =\frac{\alpha l}{1+(l1)\alpha }.$$
(2)
The conditions $`\beta <\frac{1}{\sqrt{l}+1}`$ and $`\beta <\alpha `$ induce a high-density phase, which is governed by the exit dynamics. In this case the current and bulk densities are given by

$$J=\frac{\beta (1\beta )}{1+(l1)\beta },\rho =1\beta .$$
(3)
There are two types of phase transitions in this system. A first-order phase transition, with discontinuous change in density, is observed when $`\alpha =\beta `$. However, the transitions between low-density (high-density) and maximal-current phases are continuous.

### 2.3 Approximate solutions of inhomogeneous TASEP with extended objects

For simplicity, let us assume that the size of the system $`N`$ is a large even number and that the special site with jumping rate $`q`$ is positioned at $`k=N/2`$. The exact position of local inhomogeneity is not important as long as it is far away from both boundaries .

The special site $`k`$ breaks the translational symmetry of original homogeneous TASEP. However, it also divides the system of size $`N`$ into two homogeneous translationally invariant sublattices of size $`N/2`$. This observation allows us to consider our model with local inhomogeneity as two coupled homogeneous TASEP with extended objects. The sublattices are coupled by a condition that the stationary currents in the left and right subsystems should be the same. As a result, to calculate the properties of TASEP with extended objects and with local inhomogeneity, we can use the known results for homogeneous TASEP. This approach has been used successfully before in different inhomogeneous asymmetric exclusion models . Note that the results of Ref.  correspond to our special case $`l=1`$, i.e., when each particle occupies only one site.

The main assumption of our theoretical method is that TASEP with local inhomogeneity can be viewed as two independent homogeneous TASEP coupled only by the requirement to have the same steady-state currents. Then each of the sublattices may exist in one of three different stationary states, and there are nine possible phases in the overall system. However, since the particle current in the left subsystem should be equal to the current in the right subsystem, it is impossible to have the maximal-current phase in one of the sublattices while the other sublattice is in the low-density or high-density phase. This observation eliminates four possible stationary phases and leaves only low-density/low-density (ld/ld), high-density/high-density (hd/hd), low-density/high-density (ld/hd), high-density/low-density (hd/ld) and maximal-current/maximal-current (mc/mc) stationary phases. In this notation, the first term corresponds to the state of the left sublattice, while the second term describes the right sublattice.

Consider in a more detail the possibility of existence of an ld/hd phase in our system. This phase may exist only when $`\alpha \beta `$ and

$$\alpha <\frac{1}{\sqrt{l}+1},\beta <\frac{1}{\sqrt{l}+1}.$$
(4)
Then, because the currents in the sublattices are the same,

$$J=\frac{\alpha (1\alpha )}{1+(l1)\alpha }=\frac{\beta (1\beta )}{1+(l1)\beta },$$
(5)
which yields

$$\beta =\frac{1\alpha }{1+\alpha (l1)}.$$
(6)
However, using the condition (4) for $`\alpha `$ we may conclude that $`\beta >\frac{1}{\sqrt{l}+1}`$, which violates the condition for existence of the high-density phase on the right sublattice. Thus the ld/hd phase also cannot be realized in our system.

The local inhomogeneity in our model is realized for different values of $`q`$. When $`q>1`$ the fast jumping rate is introduced. In this case, the dynamics near the local inhomogeneity will not affect the overall dynamics in the system since crossing the special site $`k`$ will not be a rate-limiting step. As a result, we have exactly the same phase diagram as for homogeneous TASEP with extended objects with ld/ld, hd/hd and mc/mc phases. The presence of local inhomogeneity will only modify the density profiles near the site $`k`$.

The more interesting case is $`q<1`$, when the dynamics near local inhomogeneity may determine the overall behavior of the system. Now consider in more detail the process of a particle crossing from the left sublattice to the right one. The particle first approaches the right sublattice when its left edge is on the site $`kl+1`$. Then it makes $`l1`$ jumps with rate 1 and one jump with the rate $`q`$, if these motions are allowed, and the particle transfers completely into the right subsystem. Thus, the effective rate of moving from the left sublattice to the right sublattice is given by

$$q_{eff}=\frac{l}{1/q+(l1)}=\frac{lq}{1+q(l1)}.$$
(7)
If the entrance to the left subsystem is the rate-limiting step, then the system will be in the ld/ld phase with the particle current and bulk density given in (2). This situation is realized for $`\alpha `$ less than some $`\alpha ^{}`$. Similarly, when the exit from the right subsystem determines the overall dynamics, the hd/hd phase exists with the particle current and bulk density given in (3), again for $`\beta <\beta ^{}`$. Parameters $`\alpha ^{}`$ and $`\beta ^{}`$ will be calculated explicitly below.

When the transition from the left sublattice to the right sublattice becomes the rate-limiting step, the hd/ld phase is realized in our system. The parameters for existence of this phase can be found from the condition that the current through left sublattice is equal to the current through the right sublattice and equal to the current passing through the local inhomogeneity. Following the expressions (2) and (3), the stationary current and bulk density in the left sublattice are given by

$$J_{left}=\frac{\beta _{eff}(1\beta _{eff})}{1+\beta _{eff}(l1)},\rho _{left}=1\beta _{eff},$$
(8)
while in the right sublattice,

$$J_{right}=\frac{\alpha _{eff}(1\alpha _{eff})}{1+\alpha _{eff}(l1)},\rho _{right}=\frac{l\alpha _{eff}}{1+\alpha _{eff}(l1)},$$
(9)
where $`\beta _{eff}`$ is the effective rate to exit from the left sublattice and $`\alpha _{eff}`$ is the effective rate to enter the right sublattice. The current passing through the local inhomogeneity can be written as

$$J_0=q_{eff}\frac{\rho _{left}}{l}\frac{(1\rho _{right})}{(1\rho _{right}+\rho _{right}/l)}.$$
(10)
This expression can be understood as follows. The parameter $`q_{eff}`$, which is given explicitly in (7), is an effective rate of crossing the local inhomogeneity. The factor $`\frac{\rho _{left}}{l}`$ is the probability to find a particle at site $`kl+1`$, i.e., the particle is crossing from the left sublattice to the right sublattice. Finally, $`\frac{(1\rho _{right})}{(1\rho _{right}+\rho _{right}/l)}`$ is the conditional probability for the corresponding sites in the right sublattice to be empty (as in ). The condition of stationary state implies that $`J_{left}=J_0=J_{right}`$, which yields

$`\rho _{right}`$ $`=`$ $`{\displaystyle \frac{\rho _{left}(1q_{eff})[l\rho _{left}(l1)]}{1(1q_{eff})[l\rho _{left}(l1)]}},`$
$`\rho _{right}`$ $`=`$ $`q_{eff}\rho _{left}.`$ (11)
Solving these equations leads to expressions for the densities in the right and left sublattices:

$`\rho _{left}`$ $`=`$ $`{\displaystyle \frac{(1+q_{eff})l}{2q_{eff}(l1)}}\left[1\sqrt{1{\displaystyle \frac{4q_{eff}(l1)}{l(1+q_{eff})^2}}}\right],`$
$`\rho _{right}`$ $`=`$ $`q_{eff}\rho _{left}.`$ (12)
The explicit formula for the current can be found from (10). This phase can only exist when $`\alpha >\beta _{eff}`$ and $`\beta >\alpha _{eff}`$, and using (8), (9) and (2.3) we obtain that hd/ld phase is stable for $`\alpha >\alpha ^{}`$ and $`\beta >\beta ^{}`$, where

$$\alpha ^{}=\beta ^{}=1\frac{(1+q_{eff})l}{2q_{eff}(l1)}\left[1\sqrt{1\frac{4q_{eff}(l1)}{l(1+q_{eff})^2}}\right].$$
(13)
Our theoretical results can be easily checked in some special limiting cases. For $`q=q_{eff}=1`$, and $`\alpha >\alpha ^{},\beta >\beta ^{}`$, we obtain that the bulk densities at the left and right sublattices are the same and equal to $`\rho _{left}=\rho _{right}=\sqrt{l}/(\sqrt{l}+1)`$, i.e., the system is in the maximal-current phase, in agreement with known results . In the limit of $`q0`$, the effective rate of crossing the local inhomogeneity is also very small, $`q_{eff}lq`$. In this limit the bulk densities reduce to $`\rho _{left}1q`$ and $`\rho _{right}lq`$, in agreement with intuitive expectations. These results can be understood as follows. Particles will move into the right sublattice with a rate of $`q`$, so that the spacing between the particles will be about $`1/q`$ (distance traveled before the next particle comes into the sublattice), giving a coverage density of $`lq`$. Now consider the left sublattice. It will be nearly full, so consider the backwards current of holes moving through the slow site. Holes will enter the left sublattice at a rate $`q`$, so the hole density will be $`q`$. Hence the coverage density in the left sublattice will be $`1q`$. Another interesting limit is $`l1`$, which corresponds to the case of very large particles. In this limit we obtain $`J1/l`$, as was found before for the homogeneous model .

It is interesting to note that our theory predicts a “mixed” nature of phase transitions from hd/ld to ld/ld or hd/hd phases. At one sublattice the transition will be continuous, while at the other sublattice there is a jump in density, which corresponds to a first-order phase transition.

## 3 Monte Carlo simulations and discussions

Our theoretical approach gives correct results in limiting cases. However, to check the overall validity of our approximate theory we performed Monte Carlo computer simulations.

Monte Carlo simulations were implemented similarly to those in . Random sequential updating was used. Suppose that at some time the lattice contains $`M`$ particles. Then in the next Monte Carlo step (MCS), $`M+1`$ particles are chosen at random, in sequence, to attempt moves. They are selected from a pool containing the $`M`$ particles on the lattice plus a free particle that can enter the system with probability $`\alpha `$ if the first $`l`$ sites are empty. Particles on lattice site $`k`$ advance with probability $`q`$ if site $`k+l`$ is empty. Particles on site $`N`$ leave the system with probability $`\beta `$. All other particles advance if they have room to move. Simulations were begun with an empty lattice and run until steady state was reached. For particles with $`l=12`$, the size studied in greatest detail, the system size was $`N=2500`$, and $`3\times 10^5`$ MCS were allowed to reach steady state. For other particle sizes, the system was run for at least $`100N`$ MCS to reach steady state. After steady state was attained, systems were simulated for an additional $`1.5\times 10^6`$ Monte Carlo steps, during which the current leaving the system was recorded constantly, and the particle density at each lattice site was sampled every 100 MCS. Continuous time Monte Carlo  was used to reduce the computation time required.

The phase diagram was determined in detail from simulations for $`l=12`$, $`q=0.2`$, and $`N=2500`$ (see Fig. 2). From the theory, we expected the ld/ld phase to occur for $`\alpha <\alpha ^{}=0.1338`$, $`\alpha <\beta `$, the hd/hd phase to occur for $`\beta <\beta ^{}=0.1338`$, $`\beta <\alpha `$, and the hd/ld phase to occur for $`\alpha ,\beta >0.1338`$. Simulations were consistent with these expectations. Sample density profiles from each phase are shown in Fig. 3. The theoretical predictions for currents and bulk densities in each half of the system also matched the simulations to within $`1\%`$ (see Fig 4).

The phase transitions were identified as follows. Near the transition between the hd/hd phase and the hd/ld phase, in the right half of the system, $`\alpha _{eff}\beta `$. This situation is expected to produce shock waves in the right half of the system (see ). We detected these shock waves in simulations by observing the approximately linear density profile that they produce in the right half of the system. An example of a density profile for the hd/hd to hd/ld transition is shown in Fig. 5. Similarly, the left half of the system exhibits an approximately linear density profile near the ld/ld to hd/ld transition, and both halves of the system have approximately linear density profiles near the ld/ld to hd/hd transition. Visual inspection of the density profiles allows the transition lines to be localized to within about $`0.02`$ units. Other techniques would be required to find the transition lines more precisely. For example, we could determine the boundaries of the hd/ld phase by measuring the current as a function of $`\alpha `$ or $`\beta `$ and determining at which $`\alpha `$ or $`\beta `$ the current saturates.

Simulations were conducted in each phase for other values of $`l`$ ($`l=2,5,20`$) to confirm the theoretical predictions. System sizes $`N`$ were always at least $`200l`$, in keeping with the assumption of large $`N`$. The theoretical predictions for the current and bulk densities were generally more accurate for larger particles. The greatest discrepancies between theory and simulations were seen for $`l=2`$ in the hd/ld phase, but observed discrepancies were always less than $`10\%`$. These results show that our approximate theory becomes more exact for larger-size particles. This is because the particle-particle correlations decrease with increasing $`l`$. Consider a lattice consisting of $`N=nl`$ sites with $`M`$ particles (each of size $`l`$). If we assume that particles are uniformly distributed, we can easily estimate the average number of empty sites $`x_l`$ between two neighboring particles: $`x_l=\frac{nlMl}{M}=lx_1`$. Thus the number of empty sites between the particles grows linearly with the particle size $`l`$. As result, particles correlate to lesser degree with each other, and our theoretical mean-field approach becomes more accurate. Similar conclusions have been reached in Ref. .

The phase diagram and phase transitions in TASEP with extended objects and with a local inhomogeneity can also be easily explained with the help of the domain-wall theory . As we discussed above, the introduction of local inhomogeneity divides the system into the two homogeneous sublattices, and the phase behavior in each sublattice is determined by the dynamics of domain walls. Thus there are always two domain walls present in the system. In the hd/ld phase at stationary state the domain walls are localized near the entrance and the exit of the system, while in ld/ld stationary phase the domain walls can be found diffusing near the local inhomogeneity and the exit. Similarly, in hd/hd stationary phase the domain walls are near the entrance and the local inhomogeneity. The phase transitions can be associated with the change of position of domain walls. For example, the transition from ld/ld phase to hd/ld phase (when $`q<1`$) is the consequence of the motion of the domain wall from the position of local inhomogeneity to the entrance, while the second domain wall stays at the exit. This observation also explains the “mixed” nature of this phase transition.

## 4 Summary and conclusions

We investigated the effect of inhomogeneity in totally asymmetric simple exclusion processes with particles of arbitrary size. Specifically, the simple model with one inhomogeneous jumping rate in the middle of the system was considered. The model was solved using a simple approximate method. Our analytical approach was based on two assumptions. The first one is the fact that the local inhomogeneity divides the system into the two coupled homogeneous systems, for which very precise mean-field solutions already exist. In the second approximation we assumed that the dynamics of particles in each sublattice is independent, i.e., we neglected the correlations between the particles in the two coupled subsystems. Using these assumptions, we calculated explicitly the phase diagrams, bulk densities, and stationary currents.

Our theoretical results indicate that for a fast jumping rate ($`q1`$), the phase diagrams, currents and bulk densities do not change as compared with homogeneous TASEP with extended objects. However, for the slow jumping rate, the situation is different. As in the homogeneous TASEP, three stationary phases are found, although phase boundaries, stationary currents and bulk densities change significantly. Our theoretical predictions are well supported by extensive Monte Carlo simulations. Some small deviations between the analytical theory and computer simulations are due to the neglect of correlations near local inhomogeneity in our theoretical approach. However, the precision of our theoretical predictions increases with the size of the particles. The microscopic nature of phase diagram was analyzed using phenomenological domain-wall theory.

Although we considered only totally asymmetric exclusion processes, our theoretical arguments can be easily extended to more general partially asymmetric exclusion processes where particles can hop in both directions. However, the important parameter of effective crossing rate from the left sublattice to the right should be calculated using more general expressions, $`q_{eff}=l/\tau `$, where $`\tau `$ is mean first-passage time for the particle to move to the right sublattice. It can be calculated using standard procedures .

ABK acknowledges support from the Camille and Henry Dreyfus New Faculty Awards Program (under Grant No. NF-00-056), from the Welch Foundation (under Grant No. C-1559), and from the US National Science Foundation through the grant CHE-0237105. KHL acknowledges support from the NSF through grant BES-0120315. LBS was supported by a Corning Foundation Fellowship and IBM Fellowship. This research was conducted using the resources of the Cornell Theory Center, which receives funding from Cornell University, New York State, federal agencies, foundations, and corporate partners.

## References


format: text

text:
# Is there a dichotomy in the Dark Matter as well as in the Baryonic Matter properties of ellipticals?

## 1. Background and new evidence

There are several lines of evidence for a dichotomy in the properties of early-type galaxies: fainter systems have pointed (disky) isophotes and central power-law surface brightness profiles, while bright galaxies are boxy and show central cores (Nieto & Bender 1989, Faber et al. 1997). This dichotomy has been interpreted in an evolutionary framework: disky/faint systems have not experienced merger events in the recent past (Nieto & Bender 1989), or alternatively are remnants of gas-rich merging events (Faber et al. 1997), while bright/boxy systems are probable merger remnants (Nieto & Bender 1989, Faber et al. 1997). This scheme is supported by X-ray properties of early-types (Pellegrini (1999) showed that faint/disky/power-law early-type galaxies are also fainter in X-ray luminosity, while bright/boxy/core galaxies are X-ray bright) and GCs number densities (Kissler-Patig 1997). What then is the actual mechanism which has triggered the evolution of both the stellar and hot gas components in galaxies?

In Fig. 1 we plot the global M/L radial gradients, $`\mathrm{\Delta }\mathrm{\Gamma }/\mathrm{\Delta }\mathrm{R}`$ ($`\mathrm{\Gamma }=M/L_B`$), based on planetary nebulae kinematics and long-slit spectroscopy archive data, as a function of the intrinsic absolute magnitude, the isophotal shape parameter $`a_4`$, and the $`\gamma `$ parameter, i.e. the slope of the surface brightness profile in the galaxy core ($`R^\gamma `$). Fig. 1 suggests a general regularity of the $`M/L`$ gradients with respect to the structural parameters for the majority of the galaxies in the sample, except for a few cases (open symbols): these are noted in literature as interacting candidates since they show dynamical peculiarities suggesting they are not in equilibrium. If we exclude this subsample, with very steep “apparent” $`M/L`$ gradients, we see that smaller gradients ($`\mathrm{\Delta }\mathrm{\Gamma }_B/\mathrm{\Delta }R0.8`$) are found for systems with faint total magnitudes ($`M_B>20`$), mostly disky ($`100\times a_4/a>0.2`$) and power-law ($`\gamma >0.15`$), while bright/boxy/core galaxies show larger gradients ($`0.8<\mathrm{\Delta }\mathrm{\Gamma }_B/\mathrm{\Delta }R2.7`$). We have found these trends significant at better than 95% c.l. via the Spearman Rank test.

## 2. Discussion and Conclusions

Following an earlier suggestion (Capaccioli et al. 2002), we have found that the dichotomy in the early-type galaxies with respect to their structural parameters and X-ray properties, see Pellegrini 1999) seems to correspond to a trend in the $`M/L`$ global radial gradients. This possibly indicates that the dark matter has triggered the evolution of both the stellar and hot gas components in galaxies. For instance, Eskridge et al. (1995) and Matsushita (2001) have suggested that the correlation of the hot gas assembly in early-type galaxies and the depth of the potential wells could explain the correlation of the X-ray luminosity with the shape and dynamical parameters ($`a_4`$, axial ratio and central velocity dispersion). Here we can confirm that faint/disky galaxies do, indeed, have shallower potential wells when compared to bright/boxy galaxies.

## References

Capaccioli, M., Napolitano, N. R., Arnaboldi, M. 2002, Sakharov Conference of Physics, Moskow, June 2002, in press, \[preprint:astro-ph/0211323\]

Eskridge, P.B., Fabbiano, G., & Dond- Woo, K. 1995, ApJ, 442, 523

Faber, S.M., et al. 1997, AJ, 114, 1771

Kissler-Patig, M., et al., 1997, A&A, 319, 83

Matsushita, K. 2001, ApJ, 547, 693

Nieto, J.L., & Bender, R. 1989, A&A, 215, 266

Pellegrini, S. 1999, A&A, 351, 487


text:
# The distance to the giant elliptical galaxy NGC 5128 Based on observations collected at the European Southern Observatory, Paranal, Chile, within the Observing Programmes 63.N-0229, 65.N-0164, 67.N-503, 68.B-0129 and 69.B-0292 and at La Silla Observatory, Chile, within the Observing Programme 64.N-0176(B).

## 1 Introduction

Thanks to its proximity NGC 5128 (Centaurus A) has attracted lots of attention. It is the closest representative of the active radio galaxy and giant elliptical galaxy class of objects and is among the nearest AGNs. The latest review (Israel israel98 (1998)) summarises its characteristics and it is a good starting point for the vast literature about this galaxy. It also gives a summary of distance determinations, which started converging in the last decade from the early highly uncertain values ranging between 2.1 and 8.5 Mpc (Sersic sersic58 (1958), Sandage & Tammann sandage&tammann74 (1974)) towards a less uncertain range between 3.2 Mpc (Hui et al. hui+93 (1993)) and 4.2 Mpc (Tonry et al. tonry+01 (2001)).

The recent determinations of distance to NGC 5128 applied a range of methods. Hui et al. (hui+93 (1993)) used the planetary nebula luminosity function to derive a distance modulus of $`(\mathrm{m}\mathrm{M})_0=27.73\pm 0.14`$ mag. The globular cluster luminosity function was analysed by Harris et al. (harris+88 (1988)) yielding $`(\mathrm{m}\mathrm{M})_0=27.53\pm 0.5`$ mag. More recent globular cluster searches in this galaxy (e.g. Rejkuba rejkuba01 (2001), Peng pengPhD (2003)) will allow a more precise determination of distance with this method through a much better sampled globular cluster luminosity function.

The luminosity of the red giant branch (RGB) tip stars in the $`I`$-band is a recognised distance indicator (e.g. Lee et al. lee+93 (1993)). In NGC 5128 it was first used by Soria et al. (soria+96 (1996)) who resolved the stellar halo using HST+WFPC2. They derived a distance modulus of $`(\mathrm{m}\mathrm{M})_0=27.86\pm 0.16`$ mag for WF chips and $`(\mathrm{m}\mathrm{M})_0=27.76\pm 0.16`$ mag for the PC chip, and adopted a mean distance modulus of $`(\mathrm{m}\mathrm{M})_0=27.8\pm 0.2`$ mag. More recently, Harris et al. (harris+99 (1999)) used deeper HST+WFPC2 photometry in a less crowded field to derive the distance. Their RGB tip luminosity analysis resulted in $`(\mathrm{m}\mathrm{M})_0=27.98\pm 0.15`$ mag. Moreover, the same authors adjusted the Hui et al. (hui+93 (1993)) distance modulus to $`(\mathrm{m}\mathrm{M})_0=27.97\pm 0.14`$, increasing it by 0.2 mag in order to correct for the contemporary Local Group distance scale and the M31 distance of $`(\mathrm{m}\mathrm{M})_{M31}=24.5`$ (van den Bergh vandenbergh95 (1995), Harris harris99 (1999)).

More than the planetary nebula luminosity function, the surface brightness fluctuations (SBF) method received several revisions of its calibration zero point. Tonry & Schechter (tonry&schechter90 (1990)) first derived $`(\mathrm{m}\mathrm{M})_0=27.48\pm 0.06`$ using I-band SBF. This value was subsequently revised to $`(\mathrm{m}\mathrm{M})_0=27.71\pm 0.10`$ by Israel (israel98 (1998)) who used the results from Tonry (tonry91 (1991)) and then to $`(\mathrm{m}\mathrm{M})_0=28.18\pm 0.07`$ by Marleau et al. (marleau+00 (2000)) after Tonry et al. (tonry+97 (1997)). Most recently, Tonry et al. (tonry+01 (2001)) report $`I`$-band SBF yielding a distance modulus of $`(\mathrm{m}\mathrm{M})_0=28.12\pm 0.15`$. In their list of nearby galaxies with SBF distance measurements, NGC 5128 occupies $`9^{th}`$ place and is the nearest giant elliptical galaxy. The distance modulus determinations from the literature are summarized in Table 1.

In this paper I use Mira variables from the long period variable star catalogue in NGC 5128 (Rejkuba et al. 2003a ) and the $`H`$ and $`K`$-band luminosity functions to derive independent measurements of distance to NGC 5128. Data are briefly described in Sect. 2, and in Sect. 3 $`J`$, $`H`$, and $`K`$-band luminosity functions are analysed. The distance to NGC 5128 is derived from the RGB tip magnitude in Sect. 4.1. The Mira period-luminosity (PL) diagram is constructed in Sect. 4.2 and used to determine the distance to NGC 5128. Finally the results are summarized in Sect. 5.

## 2 The data and photometry

Data used here were described in detail in Rejkuba et al. (rejkuba+01 (2001) and 2003a ). They consist of a set of 20 $`K_s`$-band epochs, 1 $`J_s`$, and 1 $`H`$-band image of a field located $`17\mathrm{}`$ north-east from the center of NGC 5128 (Field 1). The second field (Field 2) was observed once in the $`J_s`$ and $`H`$-bands and 24 times in the $`K_s`$-band. It is located $`9\mathrm{}`$ south from the galactic center. All the observations, except one $`K_s`$-band epoch of Field 2, were taken with the ISAAC near-IR array at UT1 (Antu) telescope at ESO Paranal Observatory in service mode. The additional $`K_s`$-band observation of Field 2 was secured during an observing run with the SOFI near-IR array on the NTT at ESO La Silla Observatory under exceptional seeing conditions.

For the details about the photometry and completeness simulations the reader is referred to the above mentioned papers. The photometry is complete more than 50% in the $`K_s`$ and $`H`$-bands for stars brighter than 22.5 for Field 1 and 21.6 for Field 2. The 50% completeness limit for the $`J_s`$-band is 23.25 and 22.5 mag for Field 1 and 2, respectively.

Multi-epoch $`K_s`$-band observations span a 1197 day interval and were used to search for variable stars. A total of 1504 red variables were detected in the two halo fields. For 1146 variables with at least 10 good measurements, periods, amplitudes and mean $`K_s`$ magnitudes were determined using Fourier analysis and non-linear sine-curve fitting algorithms. Almost all these variables belong to the class of long period variables (LPVs) with Mira and semiregular variable stars. The complete catalogue of all LPVs is presented in Rejkuba et al. (2003a ) to which reader is referred for discussion of completeness and accuracy of period and amplitude determination.

The ISAAC filters used are $`J_s`$, $`H`$ and $`K_s`$. The $`J_s`$ filter was preferred over $`J`$, due to the red leakage of the latter. It is centered at 1.24 $`\mu `$m and has a width of 0.16 $`\mu `$m. The transformation to the $`J`$-band of the LCO photometric system (Persson et al. persson+98 (1998)) was obtained through the following transformation (Chris Lidman private communication):

$$J_{LCO}=Js_{ISAAC}+0.033(JsKs)_{ISAAC}0.022$$
(1)
The ISAAC $`H`$ and $`K_s`$ filters are on the LCO photometric system.

It should be noted that to all the $`K_s`$-band magnitudes from the catalogue (Rejkuba et al. 2003a ) a constant of 0.1 mag has been subtracted. A comparison of the color-color diagrams with the intrinsic colors of Galactic dwarfs and giants (Bessell & Brett bessell+brett88 (1988)) has revealed an error in the aperture correction of the $`K_s`$-band photometry. For reference the $`JH`$ vs. $`HK`$ color-color diagram is shown in Fig. 1 with the Bessell & Brett (bessell+brett88 (1988)) fiducials overplotted. The solid line indicates the intrinsic colors of early to late type giants and the long-dashed line is for dwarfs. Foreground Galactic stars are found along these lines. The large majority of the stars in NGC 5128 are found close to the locus of late type giants and redwards from there, where long period variables and carbon stars are located. NGC 5128 photometry has been de-reddened adopting $`\mathrm{E}(BV)=0.11`$ (Schlegel et al. schlegel+98 (1998)) and the Cardelli et al. (cardelli+89 (1989)) reddening law and it has been transformed to the Bessell & Brett (bessell+brett88 (1988)) $`JHK`$ photometric system using the transformation equations available from the 2MASS web page111http://www.astro.caltech.edu/$``$jmc/2mass/v3/transformations/.

## 3 Luminosity functions

The $`K_s`$-band luminosity functions for all the stars in Field 1 (left) and 2 (right) are shown in Fig. 2 as solid line histograms. These are not corrected for incompleteness nor for Galactic foreground contamination. The number of foreground Milky Way stars can be obtained from the color-color diagram (Fig. 1) or from optical-near-IR color-magnitude diagrams published by Rejkuba et al. (rejkuba+01 (2001)). There are 194 stars with near-IR colors consistent with dwarf or early type giants in the Field 1 color-color diagram. Similarly, $`192`$ sources have blue colors in the $`VK_s`$ CMD (see Rejkuba et al. rejkuba+01 (2001)). They are plotted as a dotted histogram in Fig. 2. These stars are mostly foreground Galactic stars as well as young blue and red supergiants in Field 1. The variable star luminosity function is plotted as a shaded histogram. It is clear that majority of the variables have luminosities brighter than the RGB tip (see below) and are in the thermally pulsing asymptotic giant branch (TP-AGB) phase. In the lower panels the total luminosity functions are compared with the luminosity functions of the sums of the contributions of foreground stars and LPVs.

There are $`2176`$ stars in Field 1 above the tip of the RGB ($`K_s<21.24`$, see next section) detected in at least 3 $`K_s`$-band images. $`198`$ of these sources have not been detected in $`J_s`$ and/or $`H`$-band frames. According to the $`VK_s`$ CMD, 192 stars have blue colors and they belong either to the young population in Field 1 or to the Galactic foreground. Among the remaining 1984 stars 426 are variable stars (accounting for 71% of variables detected in this field).

In Field 2 there are 6072 sources brighter than $`K_s=21.24`$ detected in at least $`3`$ $`K`$-band images. $`351`$ of these sources have not been detected in $`J`$ and/or $`H`$-band frames. Of these $`217`$ objects have blue colors and belong to the Galactic foreground contamination. Of the remaining $`5855`$ stars brighter than the RGB tip $`829`$ are asymptotic giant branch (AGB) variables (accounting for 92% of detected variables in this field).

Some of the bright non-variable AGB giants could be blends of 2 RGB tip stars. Renzini (renzini98 (1998)) provided the equations that enable the calculation of the expected number of stars in a given evolutionary phase ($`N_j`$) in an image in which the total sampled luminosity is $`L_T`$:

$$N_j=B(t)L_Tt_j$$
(2)
$`B(t)`$ is the specific evolutionary flux of the population or the number of stars entering or leaving any post main sequence evolutionary stage per year and per solar luminosity of the population. It is a very weak function of age and ranges from $`0.5\times 10^{11}`$ to $`2.2\times 10^{11}`$ stars L$`{}_{}{}^{1}{}_{}{}^{}`$ yr-1, as age increases from 10 Myr to 15 Gyr. $`t_j`$ is the duration of the evolutionary phase. For a solar-metallicity, 15 Gyr old population $`L_T0.36L_K`$. We can estimate the total sampled luminosity in the two fields by summing the observed counts on the dark subtracted, flat-fielded, sky-subtracted and calibrated frames. Using the combined frames for the reference epoch observed on JD=2451734, the total sampled luminosity is $`5\times 10^7`$ L☉ and $`7\times 10^7`$ L☉ for Field 1 and 2, respectively. From the surface brightness profile of Dufour et al. (dufour+79 (1979)) the total luminosities are $`5.9`$ and $`28.7\times 10^7`$ L☉, while from that of Mathieu et al. (mathieu+96 (1996)) they are $`9.5`$ and $`31\times 10^7`$ L☉. Our measured value for Field 1 is close to that of Dufour et al., but the measured luminosity of Field 2 is much lower.

From Table 1 of Renzini (renzini98 (1998)) the expected number of LPVs in an old near-solar metallicity population is $`5/10^6`$ L☉. The total sampled luminosity of Field 1 is 50-100 times higher, yielding 250-500 LPVs in agreement with the 437 LPVs detected in that field (Rejkuba et al. 2003a ). In Field 2 our measured luminosity would suggest $`350`$ LPVs, but the surface brightness measurements from the literature predict $`1500`$ LPVs. For comparison, we have detected 903 variable stars and measured periods for 709 LPVs in Field 2.

The probability that a pixel contains two stars in phase $`j`$ is $`N_j^2`$ and the number of 2-star blends in a frame with a total number of pixels $`N^{pix}`$ is:

$$N_{2j}=N_j^2N^{pix}=[B(t)L_Tt_j]^2/N^{pix}$$
(3)
The resolution of ground based images is worse than 1 pixel and the sampled luminosity has to be divided by the available number of resolution elements. Each resolution element is $`0\stackrel{}{.}31`$ for Field 1 and $`0\stackrel{}{.}36`$ for Field 2, corresponding to the best seeing epochs used to detect sources. Hence the expected number of blends of two RGB tip giants in Field 1 is between 100 and 400 for the total sampled luminosity between 5 and $`9.5\times 10^7`$ L☉. In Field 2 it is $`250`$ if the measured luminosity from the image is assumed, but it rises to $`4900`$ for $`31\times 10^7`$ L☉ from Mathieu et al. The expected number of blended stars in the vicinity of the RGB tip is low enough to measure the RGB tip discontinuity in the luminosity function of Field 1, but may smear out that feature in Field 2. Taking the highest probable number of blends, there are at least 1150 non-variable bright AGB giants in Field 1 and some 150 in Field 2. Rejkuba et al. (2003a ) made extensive simulations and estimated the completeness of the detection of LPVs as a function of their magnitude, period and amplitude. These simulations show a similar completeness level of the LPV catalogue in the two fields implying that the number of bright extended giant branch non-variable stars is larger in Field 1. However, these stars could be semi-regular small-amplitude ($`\mathrm{\Delta }K0.3`$) variables for which the LPV catalogue is less complete. Their near-IR properties will be discussed together with near-IR properties of the LPVs in a forthcoming paper (Rejkuba et al. in preparation).

LPVs account for at least 26% and 70% of the extended giant branch population in Fields 1 and 2, respectively. Their high luminosity cut-off extends beyond $`\mathrm{M}_\mathrm{K}=8.7`$. The maximum brightness achieved by AGB stars is evidence for the presence of an intermediate-age population.

## 4 NGC 5128 distance

In this section the distance to NGC 5128 is measured using two different methods: (i) the brightness of the RGB tip in the $`K`$ and $`H`$-bands and (ii) the Mira period-luminosity relation.

### 4.1 RGB tip

The $`J_s`$, $`H`$ and $`K_s`$-band luminosity functions of Field 1 stars are presented in Fig. 3. These are the deepest luminosity functions for a stellar population belonging to a galaxy beyond the Local Group. The solid lines are used for completeness-corrected luminosity functions and dotted lines for the observed ones. Vertical dotted line indicates 50% completeness limit of the photometry. Short thick marks indicate the measured RGB tip magnitudes in Field 1.

In the bottom panels of Fig. 3 I show the luminosity functions convolved with the Sobel edge detection filter with a $`(2,0,+2)`$ kernel (Lee et al. lee+93 (1993)). They are used as localised slope estimators with 2-point smoothing.

There is a very sharp peak in the Field 1 response function curve at $`J=22.54\pm 0.06`$, $`H=21.54\pm 0.06`$, and $`K_s=21.24\pm 0.05`$. For easier comparison with the models and literature, the $`J`$-band measurement is transformed to the LCO photometric system according to Eqn. 1. These maxima appear at the position where there is a strong change in the slope of the RGB luminosity function and it corresponds to the magnitude of the tip of the RGB. The error in the RGB tip magnitude is obtained by averaging over 100 different measurements of Sobel kernel peak magnitude, where luminosity function sampling is varied by 0.001 mag over 0.08 to 0.2 mag.

The RGB tip feature is smoothed over a couple of tenths of magnitude in Field 2 (Fig. 4). From the 2-point Sobel kernel, the RGB tip in Field 2 is centered around $`K_s=21.22\pm 0.15`$ mag. The broadening and smoothing of the peak is due to much shallower completeness limits (note the short thick marks indicating the RGB tip magnitude measured in Field 1), resulting in larger photometric errors, as well as due to contribution of blends, the presence of AGB stars and a larger metallicity spread in this field.

The tip of the RGB in the $`K`$-band is a function of metallicity. Its dependence has been measured empirically by Ferraro et al. (ferraro+00 (2000)) using a set of Galactic globular clusters with a range of metallicities. They obtained the following relation:

$$\mathrm{M}_\mathrm{K}^{\mathrm{TRGB}}=(0.64\pm 0.12)[\mathrm{M}/\mathrm{H}](6.93\pm 0.14)$$
(4)
The mean metallicity for the stars in the NGC 5128 halo has been inferred to be around $`[\mathrm{M}/\mathrm{H}]=0.4`$ dex (Marleau et al. marleau+00 (2000), Harris & Harris harris+harris00 (2000), Rejkuba rejkubaPhD (2002)). With this estimate, and assuming a reddening of $`\mathrm{E}(BV)=0.11`$ (Schlegel et al. schlegel+98 (1998)), corresponding to $`\mathrm{A}_\mathrm{K}=0.039`$ (Cardelli et al. cardelli+89 (1989)), the empirical calibration of the RGB tip $`K`$-band magnitude can be used to derive the distance to NGC 5128. The difference between the photometric bands is of the order of 0.01 mag (Persson et al. persson+98 (1998)). This yields a distance modulus of $`(\mathrm{m}\mathrm{M})_0=27.87\pm 0.16`$. The quoted error does not include the uncertainty in reddening and in mean metallicity of the stars. A higher mean metallicity by 0.1 dex implies a larger distance modulus by 0.06 mag. The $`K`$-band is rather insensitive to extinction, and reddening as high as $`\mathrm{E}(BV)=0.35`$ would lower the distance modulus by 0.088 mag.

The above measurements can be compared to that of the RGB tip brightness in the LMC and in the Galactic Bulge. For example, Cioni et al. (cioni+00 (2000)) have measured $`K_s^{TRGB}(LMC)=11.94\pm 0.04`$ and $`J^{TRGB}(LMC)=13.06\pm 0.02`$; these values are not corrected for metallicity. The difference between the RGB tip magnitudes in the LMC and NGC 5128 imply distance moduli differences of $`\mathrm{\Delta }(\mathrm{m}\mathrm{M})_{K_s}=9.26\pm 0.06`$ mag and $`\mathrm{\Delta }(\mathrm{m}\mathrm{M})_J=9.39\pm 0.06`$ mag. Adopting a distance modulus of 18.5 for the LMC (Alves et al. alves+02 (2002)) yields distance moduli for NGC 5128 of $`27.76`$ and $`27.89`$ using $`K_s`$ and $`J`$ bands, respectively. While this is within errors consistent with the above quoted $`27.87\pm 0.16`$, it should be remembered that the higher metallicity of the NGC 5128 stars with respect to the LMC, makes the RGB tip magnitude brighter in $`K`$-band (e.g. Zoccali et al. zoccali+03 (2003)), hence the relative difference of distance modulus is a lower limit in the $`K_s`$-band.

The $`J`$-band RGB tip magnitude is expected to decrease slowly as the metallicity increases, while the tip luminosity in the $`H`$-band appears to be fairly constant in the metallicity range $`1.0[Fe/H]0.0`$ (Zoccali et al. zoccali+03 (2003)). Hence the latter could be used for distance determination. Schulte-Ladbeck et al. (schulte-ladbeck+99 (1999)) have investigated the behaviour of the $`H`$-band tip luminosity with age and metallicity. They recommend $`M_H^{TRGB}=5.5\pm 0.1`$ for $`2.3<[\mathrm{Fe}/\mathrm{H}]<1.5`$. Bertelli et al. (bertelli+94 (1994)) models suggest $`M_H^{TRGB}6`$ mag for higher metallicity, while $`M_H^{TRGB}=6.4\pm 0.2`$ in the Galactic Bulge (Zoccali et al. zoccali+03 (2003)). The latter yields a distance modulus to NGC 5128 of $`27.9\pm 0.2`$. Due to high extinction, small number of stars and possible differential reddening in the Galactic Bulge field observed, the $`M_H^{TRGB}`$ value is rather uncertain. However, it yields a distance modulus to NGC 5128 of $`27.9\pm 0.2`$ in excellent agreement with that derived from the $`K`$-band RGB tip magnitude.

### 4.2 Mira PL relation

In Fig. 5 I plot the mean magnitudes $`vs.`$ logarithm of periods for all the variables with 10 or more data points for which derived periods had a significance parameter from Fourier analysis $`<0.7`$ (Rejkuba et al. 2003a ). Field 1 variables are shown with filled and Field 2 with open symbols. The mean magnitudes are corrected for extinction and transformed to the SAAO photometric system (Carter carter90 (1990)) in which most Mira studies are published.

All the variables with $`K_0`$ magnitudes brighter than $`19.0`$ have been carefully checked on the reference (best seeing) image. The two Field 1 LPVs with periods around 430 days, both of which are brighter than $`K_s<18.0`$ are located in a highly crowded areas. The fainter of the two is actually blended with a background galaxy and the brighter might be a blend of more than two stars. Also the third brightest Field 1 LPV could be too bright due to contributing light from the nearby stars. LPVs which are blended with nearby stars should have smaller amplitudes due to larger relative contribution of the neighbour at minimum phases. Only one of the eight Field 2 bright LPVs ($`K_0<19.0`$) is a possible blend of two stars.

There are two sequences in the PL diagram in Fig. 5. Multiple parallel sequences in the PL diagram of LPVs were reported by Wood & Sebo (wood&sebo96 (1996)), Wood (wood00 (2000)) and Cioni et al. (cioni+01 (2001)) in the Magellanic Cloud. The more populated, longer period sequence corresponds to sequence C of Wood (wood00 (2000)) where Mira variables are found. The shorter period sequence B discussed by Wood (wood00 (2000)) is approximately 1.3 mag brighter, while the shorter sequence in Fig. 5 lies only 0.5-1 mag brighter. Moreover, it is is somewhat steeper, indicating a possible problem of aliasing periods.

Actually, a close inspection, shows that most of the stars on the shorter period sequence have periods close to $`1/2`$ year and can be fitted similarly well with periods that are twice as long, and which place them on the Mira sequence. It is interesting to note that this increases also the amplitude of the fit. An example is shown in Fig. 6. In the upper panel the light curve is folded with the original, shorter, period and with the longer period on the bottom. Carefully phased additional observations of these variables would be necessary to determine their periods unambiguously. The period distribution of the LPVs and the comparison between the two fields is discussed by Rejkuba et al. (2003b ).

In Table 2 all the variables that could be fitted with 2 aliasing periods are listed. Columns 1 to 6 list: field identification and ID number, original period, semi-amplitude, reduced $`\chi ^2`$ of the sine-curve fit, significance of the original period and the mean $`K_s`$ magnitude obtained from the sine-curve fit. The new period, semi-amplitude, reduced $`\chi ^2`$ of the sine-curve fit with this new period and the associated average $`K_s`$ magnitude, are listed in columns 7, 8, 9 and 10, respectively. Additionally I have checked for possible aliasing periods all the LPVs that are located below the Mira PL relation. In most cases no shorter period could provide a satisfactory fit to the data, except for 2 LPVs in Field 1 and 5 in Field 2. Fainter average $`K`$-band magnitudes for the LPVs that lie below the Mira PL relation are due to extinction within their circumstellar shells (see also Kiss & Bedding kiss+bedding03 (2003)). In another 6 cases $`\chi ^2`$ of the sine-curve fit was improved and a new amplitude and mean magnitude were calculated, although the period remained very close to the original value. Apparently running the sine-curve fitting algorithm in an automatic way failed to find the best fitting amplitudes in a few cases. These are also listed in Table 2. In the following distance determination these LPVs are not used. Fig. 7 displays the PL diagram in which aliasing periods ($`P_{new}`$) from Table 2 are used.

A well-defined PL relation has been found for Miras in the Large Magellanic Cloud (Glass & Lloyd Evans glass&lloydevans81 (1981), Wood wood+99 (1999)), the Small Magellanic Cloud (Cioni et al. cioni+03 (2003)), the Galactic Bulge (Glass et al. glass+95 (1995)), the solar neighbourhood (van Leeuwen et al. vanleeuwen+97 (1997)) and in Galactic globular clusters (Feast et al. feast+02 (2002)). The relation holds for both $`\mathrm{M}_{\mathrm{bol}}`$ and $`\mathrm{M}_K`$. Since Miras are very luminous, their tight PL relation makes them interesting for distance determination to other galaxies.

Calibration of the $`K`$-band PL relation relies on the LMC PL relation for Miras. The Feast et al. (feast+89 (1989)) fit to the LMC Mira PL relation is:

$$M_K=3.47\mathrm{log}P+\beta $$
(5)
The zero point of $`\beta =0.91`$ has been derived by Feast (feast96 (1996)) using LMC Miras. More recently Whitelock et al. (whitelock+00 (2000)) used Hipparcos parallaxes for Solar neighbourhood Miras to derive $`\beta =0.84\pm 0.14`$ and Feast et al. (feast+02 (2002)) obtained $`\beta =0.93\pm 0.14`$ for globular cluster Miras using the cluster distances determined from subdwarf fitting with Hipparcos parallaxes of subdwarfs (Caretta et al. caretta+00 (2000)). The straight mean of the two Hipparcos based zero points yields $`\beta =0.88\pm 0.10`$ and a Large Magellanic Cloud distance modulus of $`18.60\pm 0.10`$. Preferring instead a somewhat shorter LMC distance, $`(mM)_0^{LMC}=18.50\pm 0.04`$ (Alves et al. alves+02 (2002)), I use a zero point of $`\beta =0.98\pm 0.11`$ (where uncertainty in the LMC distance modulus has been added to the zero point error in quadrature).

Glass & Lloyd Evans (glass+lloydevans03 (2003)) have re-reanalysed the Mira PL relation using the MACHO data for the same stars as Feast et al. (feast+89 (1989)) in the LMC. They conclude that the periods of these variables did not change significantly and remained essentially constant over 2–3 decades. Their best fit to the $`K`$-band Mira PL relation is:

$$K=3.52(\pm 0.21)\mathrm{log}P+19.64(\pm 0.49)\sigma =0.13$$
(6)
In order to determine the distance to NGC 5128 only Miras with colors $`J_sK_s<1.4`$ have been selected. This ensures that reddening from circumstellar shells that may be present in redder LPVs would not dim the $`K`$-band magnitudes and also matches the colors of Mira variables used to derive the Mira PL relation in the LMC ($`(JK)<1.5`$, Feast et al. feast+89 (1989)). Moreover, only those LPVs that had the most regular periods, and could be fitted with a sine-curve with $`\chi ^2<5`$ were chosen. They are plotted in Fig. 8 together with the best linear least-square fit to the data (thick solid line):

$$K_0=3.37(\pm 0.11)\mathrm{log}P+28.67(\pm 0.29)(N=240)\sigma =0.20$$
(7)
The fits to Field 1 and Field 2 LPVs are identical within the errors. Largely, the error in the constant term is due to the extrapolation of the relation to the period of 1 day. The mean $`\mathrm{log}P`$ is $`2.5`$. If instead the fit is made at the mean period with the following functional form:

$$K_0=A(\mathrm{log}PX)+B$$
(8)
the error of the constant term drops to only $`0.02`$.

The slope of the PL relation in NGC 5128 agrees well with that in the LMC. It is within $`1\sigma `$ of the slope of the Feast et al. (feast+89 (1989)) relation and $`1.5\sigma `$ of the Glass & Lloyd Evans (glass+lloydevans03 (2003)) PL relation. This is yet another piece of evidence that the slope of the PL relation in the $`K`$-band is universal, lending confidence to the distance determination through the PL relation.

In order to determine the distance modulus of NGC 5128 from the Mira PL relation, one has to assume that the slope is universal and use the same value as in the calibrating PL relation (i.e. that in the LMC). It is straightforward then to calculate the distance modulus from the zero point difference of the two PL relations. Restricting the slope to $`3.47`$, and determining the zero-point of the relation at $`(\mathrm{log}P2.3)`$, approximately the mean period of the LMC stars used to determine the relation (Feast et al. feast+89 (1989), Glass & Lloyd Evans glass+lloydevans03 (2003)) the best fitting zero point is $`28.94\pm 0.03`$ (dotted green line in Fig. 8) with the same RMS of the fit. Using the above discussed zero point of $`\beta =0.98`$, I derive a distance modulus to NGC 5128 of $`(mM)_0=27.96\pm 0.11`$. If instead the slope of the PL relation is fixed to $`3.52`$ (Glass & Lloyd Evans glass+lloydevans03 (2003)), the resulting distance modulus to NGC 5128 is $`(mM)_0=27.93\pm 0.11`$.

Feast et al. (feast+89 (1989)) as well as Glass & Lloyd Evans (glass+lloydevans03 (2003)) noted the group of stars with periods in excess of $`420`$ days which lie some $`0.7`$ mag above the PL relation. Due to the fact that a few of them have been found to be Li-rich, it was suggested that they are more massive, hence younger and are in the Hot-Bottom Burning phase (Smith et al. smith+95 (1995)). So, restricting the LPV data in NGC 5128 to LPVs with periods shorter than 400 days, where the above mentioned LMC PL relations have been calibrated, distance moduli of $`27.95\pm 0.12`$ and $`27.92\pm 0.12`$ are derived adopting slopes of $`3.47`$ and $`3.52`$, respectively. These distance moduli are virtually indistinguishable from that determined from the full data set and they point out that most of the long period stars do follow the same Mira PL relation as shorter period Miras. They are not overluminous and in the hot-bottom burning phase.

There are a few dozen stars in Fig. 5 that are $`2\sigma `$ above the PL relation. Some of them might be brighter due to contributing light from their neighbours, but some show indications of humps in their light-curves. Similar humps have been found in MACHO light curves of the LMC Li-rich stars by Glass & Lloyd Evans (glass+lloydevans03 (2003)). However, before being able to place firm conclusions, accurate light curves over several periods should be determined for a larger statistical sample of Li-rich hot-bottom burning AGB stars.

## 5 Conclusions

The tip of the RGB was detected in Field 1 at $`K_s=21.24\pm 0.05`$, yielding a distance modulus of NGC 5128 of $`(\mathrm{m}\mathrm{M})_0=27.87\pm 0.16`$. The comparison of the $`H`$-band RGB tip luminosity in the Galactic Bulge and NGC 5128 implies a similar distance modulus (Table 3). The RGB tip in Field 2 is not a sharp feature due to a brighter completeness limit, larger photometric errors and the presence of blends and AGB stars.

A large population of stars above the tip of the RGB contains 2176 stars in the outer halo field (Field 1) and 6072 stars in the inner halo field (Field 2). Subtracting these foreground sources, detected LPVs, as well as maximum probable number of blends of two RGB tip stars, there are some 1150 and 150 non-variable stars brighter than the first ascent giant branch tip in the two fields. LPVs account for 26% and 70% of the AGB population in Fields 1 and 2, respectively. The high luminosity ($`M_K8.7`$) achieved by AGB stars is a sign of an intermediate-age population.

The first Mira period-luminosity relation outside the Local Group is presented. Miras with the best determined and most regular periods that do not have red colors ($`J_sK_s<1.4`$) were used to determine the distance of NGC 5128 from a period-luminosity relation. I derive the distance modulus of $`27.96\pm 0.11`$, adopting the LMC distance modulus of 18.50 (Alves et al. alves+02 (2002)).

The mean distance of the Miras, obtained from fields on both sides of the center of NGC 5128, is slightly larger than the mean distance of the red giants in Field 1. This indicates that the orientation of NGC 5128 is such that the north-eastern part of the galaxy is closer to us. However, the relative distance of the two fields is smaller than $`0.03`$ Mpc and negligible compared to the distance from the Milky Way. The mean of the two methods yields a distance to NGC 5128 of $`(mM)_0=27.92\pm 0.19`$ ($`D=3.84\pm 0.35`$ Mpc), very close to the $`(mM)_0=27.98\pm 0.14`$ measurement of Harris et al. (harris+99 (1999)) and within $`1\sigma `$ of the distance modulus $`27.8\pm 0.2`$ result of Soria et al. (soria+96 (1996)).

###### Acknowledgements.

I would like to thank M.-R. Cioni for useful discussions, M. Feast for suggestions regarding the Mira PL relation and to Chris Lidman for communicating the transformation equations between different near-IR photometric systems. Thanks also go to Tim Bedding for carefully reading the manuscript and for making a number of useful comments. Comments from H. Ford and the second anonymous referee improved the clarity of the presentation. I am indebted to many ESO staff astronomers who took the data presented in this paper in service mode operations at Paranal Observatory.


text:
# The Early Optical Afterglow of GRB 030418 and Progenitor Mass Loss

## 1. Introduction

Around half of all well-localized gamma-ray bursts (GRBs) have resulted in the detection of optical counterparts. This low success rate is partly due to the difficulty in obtaining prompt coordinates and so it has been argued that many GRB afterglows fade too rapidly for discovery by late time follow-up observations. It is also possible that extinction from dense circumburst environments may cut optical emission below detectable levels (Klose et al., 2003). GRB 990123 remains unique as the only burst from which prompt optical emission was detected during gamma-ray emission (Akerlof et al., 1999), despite much effort from small rapidly responding telescopes such as ROTSE-I and LOTIS (Akerlof et al., 2000; Kehoe et al., 2001; Park et al., 1999). Larger fast-slewing telescopes such as ROTSE-IIIa have since come online in an effort to achieve deeper imaging at early times. To date, only two other afterglows have been detected within 10 minutes of the burst—GRB 021004 (Fox et al., 2003c) and GRB 021211 (Fox et al., 2003b; Li et al., 2003)—and both of these were detected only after the afterglow began to decay.

In this paper, we report on early-time optical observations of GRB 030418 with the ROTSE-IIIa (Robotic Optical Transient Search Experiment) telescope and the SSO 40-inch telescope, both located at Siding Spring Observatory, Australia. We also report on early-time optical observations of GRB 030723 with the ROTSE-IIIb telescope at McDonald Observatory, Texas. Despite rapid responses to each of these bursts (211 s and 47 s respectively), we have no evidence of prompt optical counterparts. We present here a physical model that ascribes the afterglow rise to extinction in the local circumburst environment.

The ROTSE-III array is a worldwide network of 0.45 m robotic, automated telescopes, built for fast ($`6`$ s) responses to GRB triggers from satellites such as HETE-2. They have a wide ($`1\stackrel{}{\mathrm{.}}85\times 1\stackrel{}{\mathrm{.}}85`$) field of view imaged onto a Marconi $`2048\times 2048`$ back-illuminated thinned CCD, and operate without filters. The ROTSE-III systems are described in detail in Akerlof et al. (2003).

The SSO 40-inch telescope has an f/8 direct imager at a Cassegrain focus. The field of view has a $`20.8\mathrm{}`$ diameter on a Tek $`2048\times 2048`$ CCD with 24 micron pixels. The telescope can be operated unfiltered or with a range of filters. For these observations the CCD was used with $`2\times 2`$ binning giving pixels of $`1\stackrel{}{\mathrm{.}}2`$.

## 2. Observations and Analysis

On 2003 April 18, HETE-2 detected GRB 030418 (HETE-2 trigger 2686) at 9:59:18.85 UT. The first determination of its position was distributed as a Gamma-ray Burst Coordinates Network (GCN) notice at 10:02:54 UT, with a $`28\mathrm{}`$ radius error box, 205 s after the start of the burst (HETE-2 2686;2). Ground analysis improved the error radius to $`18\mathrm{}`$, and a second GCN notice was distributed at 11:43:01 UT (HETE-2 2686;4). The burst was determined to have lasted 135 s, with a fluence of $`1.2\times 10^6\mathrm{erg}\mathrm{cm}^2`$ (2–25 keV) and $`2.5\times 10^6\mathrm{erg}\mathrm{cm}^2`$ (30–400 keV), classifying it as a long, X-ray rich GRB (Shirasaki et al., 2003). The burst occurred while the Moon was bright (96% illumination) and very few telescopes reported follow-up observations.

ROTSE-IIIa responded automatically to the first GCN notice in under 6 s with the first exposure starting at 10:03:00 UT, 211 s after the burst and only 76 s after the cessation of gamma-ray activity. The automated scheduling software began a program of ten 5-s exposures followed by 90 20-s exposures. Longer exposures were not taken because the bright Moon would have saturated the images. The second GCN notice triggered ROTSE-IIIa to repeat the same sequence of observations. Analysis of the individual frames in near real-time did not reveal any new source brighter than the unfiltered limiting magnitude of $`16`$.

After the receipt of the second GCN notice, Price et al. (2003) initiated a burst response on the SS0 40-inch telescope, beginning 7139 s after the burst, with an unfiltered 100-s image, followed by a sequence of 10 300-s images taken with a Johnson $`R`$-band filter. Price et al. (2003) identified a new object at $`\alpha =10^h54^m33\stackrel{\mathrm{s}}{\mathrm{.}}674`$, $`\delta =7\mathrm{°}01\mathrm{}40\stackrel{}{\mathrm{.}}75`$ (J2000.0), at $`18.8`$ magnitude that was not on the Digitized Sky Survey red plates (Price et al., 2003).

Co-adding sets of 10 ROTSE-IIIa images revealed the optical counterpart reported by Price et al. (2003). As can be seen in Figure 1, the afterglow is barely detected in the first ROTSE-IIIa image, while a nearby $`18^{\mathrm{th}}`$ magnitude star at $`\alpha =10^h54^m32\stackrel{\mathrm{s}}{\mathrm{.}}6`$, $`\delta =07\mathrm{°}03\mathrm{}38\stackrel{}{\mathrm{.}}0`$ (J2000.0) is clearly visible. The afterglow then increases in brightness, exceeding the comparison star, before fading below our detection threshold. Using an identical analysis as for the afterglow itself (to be described below), we measured the fluctuations of the comparison star to be $`<0.1`$ magnitudes. The comparison star is detected with a S/N of $`7.7`$ in the first images, and a S/N of $`>10`$ in subsequent images.

The ROTSE-IIIa data and the first SSO 40-inch observation were taken without filters. Further SSO 40-inch observations were taken with an $`R`$-band filter. Therefore, it is important to bring both sets of measurements to the same standard photometric system. We compared each image to the standard $`R`$-band photometric calibration from the USNO-1m telescope (Henden, 2003). Unfortunately, we do not have any color information for the afterglow at the early times. Table 1 shows the results of our photometry for the afterglow, with the light curve plotted in Figure 2.

The ROTSE-IIIa images were bias-subtracted and flat-fielded in the standard way. The flat-field image was generated from 30 twilight images. We used SExtractor (Bertin & Arnouts, 1996) to perform the initial object detection, and to determine the centroid positions of the stars. We used a robust PSF fitting code to measure the photometry. We used six bright, well-measured stars from the Henden (2003) list to fit our PSF, as well as to determine an $`R`$-band magnitude zero-point. The intrinsic variations in object color generate an RMS dispersion of $`0.23`$ between our unfiltered ROTSE magnitudes and the Henden (2003) $`R`$-band magnitudes.

The SSO 40-inch images were bias-subtracted and flat-fielded with a twilight flat in the standard way. We used weighted-aperture photometry and set the magnitude zero point for each frame using between 30 and 90 Henden (2003) stars. The systematic error estimated from comparison to the Henden (2003) stars is about 0.25 magnitudes.

On 2003 July 23, at 06:28:17.45 UT, HETE-2 detected another X-ray bright GRB (HETE trigger 2777), 23 s in duration (Prigozhin et al., 2003). The ROTSE-IIIb instrument responded automatically and began taking images within 5 s of the GCN notice distribution. The first ROTSE exposure began 47 s after the burst trigger time. The system took ten 5-s images, ten 20-s images, and 40 60-s images of the burst field. Much like GRB 030418, the burst counterpart was not found in these early images to limiting magnitudes of 17–18 mag. Later observations by larger telescopes revealed a faint, fading source at $`\alpha =21^h49^m24\stackrel{\mathrm{s}}{\mathrm{.}}40`$, $`\delta =27\mathrm{°}42\mathrm{}47\stackrel{}{\mathrm{.}}4`$ (J2000.0) at $`20^{\mathrm{th}}`$ magnitude (Fox et al., 2003a). We co-added sets of ten images and applied the same PSF-fitting technique as described above. The object was not detected in our four earliest co-added images, but the last two images yield marginal detections. We derive $`19.5\pm 0.4`$ for the fifth image (S/N=2.7) and $`19.3\pm 0.4`$ for the sixth image (S/N=3.1). The resulting light curve is shown in Figure 3.

## 3. Results

Figure 2 shows the early-time light curve of GRB 030418 from the ROTSE-IIIa and SSO 40-inch observations. After our first detection, the afterglow brightness is rising. After the rise, the afterglow slowly varies around $`17.3`$ mag for about an hour, before fading following a power law.

Most GRB afterglows have been seen to decline with a sequence of one or more power laws. Two hours after the burst, this is also the case for GRB 030418. Extrapolating this late time power law decay back to the early time provokes the question of what happened to the missing optical flux. There are two main alternatives: intrinsically, there were fewer optical photons emitted by the source or, extrinsically, the optical photons were absorbed after they were created.

The first possibility is not easy to describe without making extensive assumptions about the physics of the shock front or the density of the ambient medium. We have investigated the possibility that the light curve results from the spectral break frequency coming through the optical band (Sari et al., 1998). We modeled the flux spectrum as a two component power-law as in Sari et al. (1998) for the slow cooling regime, with $`F_\nu \nu ^{1/3}`$ for $`\nu <\nu _m`$, and $`F_\nu \nu ^{(p1)/2}`$ for $`\nu >\nu _m`$, where $`\nu _m`$ is the characteristic synchrotron frequency, and $`p`$ is the spectral index of the electrons, fit to the late time power-law decline. The synchrotron frequency decreases as $`\nu _mt^{3/2}`$, with a best fit value of $`\nu _m=9.7\times 10^{14}\mathrm{Hz}`$ at 1500 s after the burst. After integrating the flux in our optical passband, we found the predicted optical peak is much too sharp and does not reproduce the smooth rollover we see in the light curve of GRB 030418, as can be seen as the dotted line in Figure 4. This fast transition is not caused by the sharpness of the frequency break, but rather the rapidity with which the break moves through the optical band.

We next investigated the possibility that we are seeing the transition from fast cooling to slow cooling, as in Sari et al. (1998). In addition to our model described above, we added the cooling frequency, $`\nu _c`$, which decays as $`\nu _ct^{1/2}`$. By fitting the initial values of $`\nu _c`$ and $`\nu _m`$, the modeled light curve is plotted as the dashed line in Figure 4. This model succeeds in describing the flat peak of the light curve, but does not describe the rise at early times. In addition, the best fit frequency values at 1500 s for $`\nu _c=4.6\times 10^{14}\mathrm{Hz}`$ and $`\nu _m=1.8\times 10^{15}\mathrm{Hz}`$ are not physically reasonable, according to the prescription in Granot & Sari (2002).

Recently, evidence has emerged linking GRBs to core-collapse supernovae of massive stars, including the detection of a SN spectrum in the afterglow of GRB 030329 (Stanek et al., 2003; Hjorth et al., 2003). A key consequence of a massive star progenitor is that the GRB occurs inside a massive stellar wind (Meszaros et al., 1998; Chevalier & Li, 1999). To date, most of the literature has focused on the generation of the shock front in a wind medium without consideration of other physical effects. Here, we investigate the possibility that the absorption of the optical photons by this circumburst environment can explain the early behavior of the GRB 030418 light curve. Although a prompt optical/UV flash might sublimate the dust (Waxman & Draine, 2000), if this were the entire story then we would likely have seen the afterglow decay at the early time.

The stellar wind density profile is $`\rho =\rho _0r_0^2/r^2`$, assuming a constant mass loss rate from the massive star. The optical depth scale then becomes $`\tau =\beta _r/r`$, and the optical flux absorbed by the circumburst medium is attenuated as $`e^{\beta _r/r}`$. If we assume that at the early time the emitting shell is moving at a roughly constant velocity with a bulk Lorentz factor of $`\mathrm{\Gamma }`$, then the distance traveled is related to the time in the Earth’s frame of reference as $`r\frac{1}{2}\mathrm{\Gamma }^2ct`$. Therefore, the attenuation as a function of time goes as $`e^{\beta _t/t}`$. As the standard model has the afterglow fading as a power law after the initial energy injection into the circumburst medium, we assume that the power law decline applies at all relevant times. We can then fit the rise and fall of the afterglow with an attenuated power law function,

$$F_\nu =F_0t^\alpha e^{\beta _t/t}.$$
(1)
We derive the three free parameters of this fitting function, $`F_0`$, $`\alpha `$, and $`\beta _t`$ empirically via a linear regression fit to the observed afterglow light curve. The best fit function is plotted in Figure 4, with $`\alpha =1.36\pm 0.02`$, and $`\beta _t=1.81(\pm 0.05)\times 10^3\mathrm{s}`$, and $`F_0`$ is an arbitrary normalization factor. The formal errors quoted are from the linear regression fit. The fit has a $`\chi ^2/\mathrm{DOF}=8.7`$. This large value is mostly due to local effects in the data at 8000 s.

From Figure 4 it can be seen that this physical model is able to represent the gross features of the light curve. We will now show that the observed value for the attenuation time-scale, $`\beta _t`$, is consistent with reasonable assumptions about the circumburst medium. To estimate the mass loss implied by $`\beta _t`$, we assume that the primary absorbing medium is dust grains similar in composition to the ISM. The opacity of the ISM can be approximately modeled as $`\kappa _\nu =(5\times 10^{13}\mathrm{cm}^2\mathrm{g}^1\mathrm{s})\nu `$ (Draine & Lee, 1984; Adams & Shu, 1985). Although there is no redshift measurement for this burst, we assume a typical redshift of $`z1`$, which puts the peak response of our CCD at $`\lambda 325\mathrm{nm}`$ in the GRB local frame of reference.

We estimate the wind velocity as $`v_{wind}100\mathrm{km}\mathrm{s}^1`$, the escape velocity from a $`50M_{\mathrm{}}`$ star with a radius $`2000R_{\mathrm{}}`$, a typical value for a red supergiant. Finally, we assume that the bulk Lorentz factor of the GRB ejecta is 100. The derived mass loss rate for the progenitor star is then:

$$\frac{dM}{dt}=1.2\times 10^3\frac{M_{\mathrm{}}}{\mathrm{yr}}\left(\frac{\mathrm{\Gamma }}{100}\right)^2\frac{\beta _t}{1800}\frac{450\mathrm{cm}^2\mathrm{g}^1}{\kappa }\frac{v}{100\mathrm{km}\mathrm{s}^1},$$
(2)
where $`\mathrm{\Gamma }`$ is the bulk Lorentz factor, $`\kappa `$ is the opacity and $`v`$ is the velocity of the stellar wind. With our assumptions, a mass loss rate of $`1.2\times 10^3M_{\mathrm{}}\mathrm{yr}^1`$ is high compared to that of late stage, high mass stars, which is typically $`10^5`$$`10^4M_{\mathrm{}}\mathrm{yr}^1`$ (Garcia-Segura et al., 1996a, b). However, such a mass loss rate may not be out of line for the extreme stellar masses required by the prevailing collapsar/hypernova models of GRB progenitors. This estimate is also dependent on the accuracy of our assumed value for the bulk Lorentz factor – a factor of three lower would reduce the mass loss by a factor of ten. Additionally, our estimate is only reasonable with a mass loss rate typical of a red supergiant, not a Wolf-Rayet star as discussed as a possible progenitor for GRB 021004 (Schaefer et al., 2003).

## 4. Discussion

GRB 030418 is one of the earliest afterglows yet imaged, with the initial detection only 76 s after the cessation of the gamma-ray activity. Unlike GRB 990123, this burst does not appear to have a prompt optical counterpart that can be attributed to the reverse shock. However, if our model of local extinction is correct, we would not expect to see any prompt emission; this model implies optical extinction of roughly 20 magnitudes at $`100\mathrm{s}`$, near the time gamma-ray emission ceased.

A backward extrapolation of the late power law decline overestimates the optical emission from GRB 030418. Our model of local dust absorption in a stellar wind medium is a useful way of characterizing the data. Unlike frequency break models, our model is able to describe the steep rise and slow rollover of the light curve. This early time behavior is far from universal, as several bursts, including GRB 990123, and GRB 021211, had more emission than predicted from the late power law decline. However, we already have evidence that the light curve behavior of GRB 030418 is not unique. Our early-time observations of GRB 030723 show a similarity to the light curve of GRB 020418. Given the observational biases against detecting such dim fading objects, it is not surprising that this class of GRB afterglows is just now being discovered as a consequence of more accurate coordinate determinations in space and more sensitive optical detectors on the ground.

One of the main consequences of our absorption model in a stellar wind medium is that some afterglows will rise very steeply in the early time. It is at this very early time that the degeneracy between our model and the frequency break models is broken. Another consquence of a dusty local environment is that the extinction in the optical bands will be much greater than in the near infrared. Prompt multi-color observations will therefore be invaluable to firmly establish if this type of initial behavior is due to optical absorption as described above.

This work has been supported by NASA grants NAG5-5281 and F006794, NSF grants AST-0119685 and 0105221, the Australian Research Council, the University of New South Wales, and the University of Michigan. Work performed at LANL is supported by NASA SR&T through Department of Energy (DOE) contract W-7405-ENG-36 and through internal LDRD funding. Special thanks to the staff at Siding Spring and McDonald Observatories, to Scott Barthlemy for maintaining the GCN and to Fred Adams for his illuminating discussion of dust absorption.