format: text

text:
# Nature of singularities in anisotropic string cosmology

## I Introduction

Superstring theory continues to be of interest as a possible candidate to unify all fundamental interactions in nature . It is known that there exist five supersymmetric perturbative string theories which are classified as the type I, type IIA, type IIB, SO(32) heterotic and $`E_8\times E_8`$ heterotic string. Recently it was found that these theories are connected by dual symmetries, which leads to the conjecture that each theory appears as one of five branches of a unified theory, called M-theory . In particular Hořava and Witten  showed that the 10-dimensional $`E_8\times E_8`$ heterotic string theory is equivalent to an 11-dimensional M-theory compactified to $`𝐌^{10}\times 𝐒^1/Z_2`$. Then the 10-dimensional spacetime is expected to be compactified into $`𝐌^4\times \mathrm{𝐂𝐘}^6`$, in which case the standard model particles are confined on the 3-dimensional brane. This gives rise to the well-known brane world scenario  where the extra dimension is noncompact and gravity is effectively 3-dimensional.

For cosmologists it is very important to test the viability of string theories by extracting various cosmological implications from them . One of such attempts is the pre-big-bang (PBB) scenario  based on the low energy effective action of string theory. In this scenario there exist two branches of solutions by assuming a $`T`$-duality. One of which ($`t<0`$) corresponds to the stage of pole-like inflation driven by the kinetic term of the dilaton field, and another ($`t>0`$) is the stage where the curvature continues to decrease. However it is difficult to smoothly connect these two branches without singularity in tree-level string action .

One is required to take into account quantum loop or derivative corrections in order to overcome such singularity problems. In fact, Antoniadis, Rizos, and Tamvakis  included a Gauss-Bonnet term to the tree-level string effective action with dilaton and modulus fields, and showed the existence of nonsingular cosmological solutions. In this case nonsingular behavior of solutions is mainly determined by the evolution of the modulus field. Therefore allowed ranges of parameters were analyzed in the absence of the dilaton in the flat Friedmann-Robertson-Walker (FRW) background  (see also ref. ). Since it is important to confirm the generality of singularity avoidance even starting from an anisotropic spacetime, several authors analyzed nonsingular cosmological solutions in the Bianchi type-I spacetime without dilaton  and with dilaton . The presence of the modulus coupled to the Gauss-Bonnet term allows the existence of nonsingular solutions unless the dilaton controls the dynamics of the system.

In order to understand how nonsingular or singular solutions appear, it is necessary to classify nature of singularities in an anisotropic background. In particular, the main determinant $`D`$ of the system is an important quantity to describe the singularities. When only the dilaton field $`\varphi `$ is coupled to Gauss-Bonnet term in the Bianchi I background, it was conjectured in ref.  that nonsingular cosmological solutions in future crosses the determinant singularity ($`D=0`$) in past when $`\dot{\varphi }`$ is positive. While this singularity was found more than ten years ago , only now we begin to understand its importance. The similar kind of singularity also appears in the context of black hole inner solutions in the presence of dilaton coupled to gravity via the Gauss-Bonnet term . In this paper we shall make detailed analysis about nature of singularities both in dilaton- and modulus-driven cosmologies in the Bianchi I background. We do not include an axion field in our analysis, but it is important to emphasize that its effect is generally vital as studied in Refs. . We will classify other kinds of singularities where at least one expansion rate diverges. These investigations are important to understand how nonsingular solutions emerge in the modulus-driven case. In addition our analysis will be useful to construct more complicated nonsingular string-inspired models in the presence of other fields such as axion.

This paper is organized as follows. In Sec. II we show background equations in anisotropic string-inspired models with dilaton or modulus fields. In Sec. III we study nature of singularities in dilaton-driven cosmology both for positive and negative $`\dot{\varphi }`$ cases. Sec. IV is devoted to the modulus-driven cosmology where both of nonsingular and singular solutions exist. We present summary and discussions in the final section.

## II The model and background equations

We begin with the action ,

$`S={\displaystyle d^4x\sqrt{g}\left[\frac{1}{2}R\frac{1}{2}(\varphi )^2+f(\varphi )R_{\mathrm{GB}}^2\right]},`$ (1)
written in the Einstein frame. Here $`R`$ is the scalar curvature and $`\varphi `$ denotes a scalar field which is either dilaton or modulus. $`f(\varphi )`$ depends on string theories, whose explicit forms are given later. We do not consider the multi-field system of dilaton and modulus fields  induced from the one-loop effective action of heterotic string theory. In addition we neglect the anti-symmetric tensor $`H_{\mu \nu \lambda }`$ and the curvature terms higher than the second order. The Gauss-Bonnet term, $`R_{\mathrm{GB}}`$, is defined as

$`R_{\mathrm{GB}}^2=R^24R^{\mu \nu }R_{\mu \nu }+R^{\mu \nu \alpha \beta }R_{\mu \nu \alpha \beta }.`$ (2)
We normalize time and spatial coordinates by the string length scale $`\sqrt{\lambda _s}`$ as $`\overline{x^\mu }=x^\mu /\sqrt{\lambda _s}`$, and the scalar fields as $`\overline{\varphi }=\varphi \sqrt{\lambda _s}`$. Hereafter we drop bars for simplicity.

Let us consider the Bianchi type-I spacetime whose metric is given by

$`ds^2=dt^2+a^2(t)dx^2+b^2(t)dy^2+c^2(t)dz^2,`$ (3)
where $`a(t),b(t),c(t)`$ are the scale factors in an anisotropic background. We define the anisotropic expansion rates, $`p(t),q(t),r(t)`$, as

$`p(t)={\displaystyle \frac{\dot{a}}{a}},q(t)={\displaystyle \frac{\dot{b}}{b}},r(t)={\displaystyle \frac{\dot{c}}{c}},`$ (4)
where a dot denotes the derivative with respect to $`t`$. It is also convenient to introduce new variables $`h,\alpha ,\beta `$ :

$`p(t)=h+\alpha +\sqrt{3}\beta ,q(t)=h+\alpha \sqrt{3}\beta ,r(t)=h2\alpha .`$ (5)
Here $`h=(p+q+r)/3`$ is an average expansion rate, which is the generalization of the Hubble parameter in an isotropic case, and $`\alpha `$ and $`\beta `$ correspond to parameters of anisotropy. The triangle in $`(\alpha ,\beta )`$-plane,

$`\alpha +\sqrt{3}\beta >h,\alpha \sqrt{3}\beta >h,\alpha <h/2,`$ (6)
which extends around the isotropic point $`\alpha =\beta =0`$, represents the regions where the Universe expands in all directions. In the outside region the Universe is contracting at least in one direction.

The dynamical equations for the background are written as

$`(1+8r\dot{f})(\dot{q}+q^2)+(1+8q\dot{f})(\dot{r}+r^2)+(1+8\ddot{f})qr+{\displaystyle \frac{1}{2}}\dot{\varphi }^2=0,`$ (7)
$`(1+8r\dot{f})(\dot{p}+p^2)+(1+8p\dot{f})(\dot{r}+r^2)+(1+8\ddot{f})rp+{\displaystyle \frac{1}{2}}\dot{\varphi }^2=0,`$ (8)
$`(1+8q\dot{f})(\dot{p}+p^2)+(1+8p\dot{f})(\dot{q}+q^2)+(1+8\ddot{f})pq+{\displaystyle \frac{1}{2}}\dot{\varphi }^2=0,`$ (9)
$`\ddot{\varphi }+(p+q+r)\dot{\varphi }8f^{}\left\{\dot{p}qr+p\dot{q}r+pq\dot{r}+pqr(p+q+r)\right\}=0,`$ (10)
together with the constraint equation,

$`pq+qr+rp+24pqr\dot{f}{\displaystyle \frac{1}{2}}\dot{\varphi }^2=0.`$ (11)
Defining a 4-dimensional vector, $`𝐱=(\dot{p},\dot{q},\dot{r},\ddot{\varphi })`$, eqs. (7)-(10) can be written in the matrix form,

$`Z𝐱=𝐲,`$ (12)
where $`𝐲=𝐲(p,q,r,\varphi ,\dot{\varphi })`$ and

$`Z=\left[\begin{array}{cccc}0& 1+8f^{}\dot{\varphi }r& 1+8f^{}\dot{\varphi }q& 8f^{}qr\\ 1+8f^{}\dot{\varphi }r& 0& 1+8f^{}\dot{\varphi }p& 8f^{}rp\\ 1+8f^{}\dot{\varphi }q& 1+8f^{}\dot{\varphi }p& 0& 8f^{}pq\\ 8f^{}qr& 8f^{}rp& 8f^{}pq& 1\end{array}\right].`$ (17)
The determinant of $`Z`$ yields

$`D=2`$ $`+`$ $`16f^{}\dot{\varphi }(p+q+r)64f^2(p^2q^2+q^2r^2+r^2p^2)+128f^2pqr(p+q+r)+128f^2\dot{\varphi }^2(pq+qr+rp)`$ (18)
$`+`$ $`1024f^3\dot{\varphi }pqr(pq+qr+rp+\dot{\varphi }^2)+12288f^4\dot{\varphi }^2p^2q^2r^2.`$ (19)
In the case of $`D0`$, the solutions of eqs. (7)-(10) are given by $`𝐱=Z^1𝐲`$. When $`D`$ vanishes, however, we cannot proceed numerical calculations further. This “determinant singularity” plays an important role in the anisotropic background .

From Eq. (11) we find the constraint

$`\alpha ^2+\beta ^2h^2+96(pqrf^{})^2,`$ (20)
in which case $`\dot{\varphi }`$ is solved as

$`\dot{\varphi }=24pqrf^{}\pm \sqrt{(24pqrf^{})^2+2(pq+qr+rp)}.`$ (21)
When $`f=0`$ anisotropy parameters are restricted in the circle, $`\alpha ^2+\beta ^2h^2`$. If the Gauss-Bonnet term is taken into account, we have the wider allowed range of anisotropy parameters given by eq. (20).

## III Dilaton-driven case

Firstly we consider the dilaton-driven case with

$`f(\varphi )={\displaystyle \frac{\lambda }{16}}e^{2\varphi },`$ (22)
where the string coupling, $`\lambda `$, takes a positive value. We set $`\lambda =1`$ in our numerical simulations. In the scenario (22) nature of singularities was analyzed in ref.  in the Bianchi I background in the case of the plus sign in the rhs of Eq. (21). Hereafter we shall make detailed analysis about the property of singularities in both signs of eq. (21).

The asymptotic behavior of solutions in past and future can be analyzed by assuming the following power-law forms for the expansion rates:

$`p=c_1|t|^s,q=c_2|t|^s,r=c_3|t|^s.`$ (23)
In order for the $`\dot{f}`$ term in eqs. (7)-(10) to have a power-law dependence, the dilaton is required to take the form

$`\varphi =\varphi _0+c_4\mathrm{ln}|t|.`$ (24)
When the contribution from the Gauss-Bonnet term is negligible, one has $`\dot{\varphi }^2=2(c_1c_2+c_2c_3+c_3c_1)|t|^{2s}`$ from eq. (11). Comparing this with eqs. (7)-(10) and (24), we find

$`s=1,c_1+c_2+c_3=\mathrm{sign}(t),c_1^2+c_2^2+c_3^2+c_4^2=1,c_4^2=2(c_1c_2+c_2c_3+c_3c_1).`$ (25)
In the absence of the dilaton ($`c_4=0`$), the solution (25) for $`t>0`$ represents the vacuum Kasner solution where the Universe is expanding in two directions and contracting in one direction. The interpretation of this solution is that large anisotropies are required at least in one dimension as $`t0`$, in order to make the spacetime curved by anisotropies.

The situation is different when the dilaton is taken into account. For example, when $`\sqrt{1/2}c_4\sqrt{2/3}`$ the Universe is expanding in all directions for $`t\mathrm{}`$, while there exist both Friedmann- and Kasner-type solutions for $`0c_4<\sqrt{1/2}`$. Note that in the limit where the Gauss-Bonnet term is negligible ($`f^{}0`$) the determinant approaches a constant value $`D=2`$ from eq. (19).

When the Gauss-Bonnet term is dominant in eq. (11) ($`|pq+qr+rp||24pqr\dot{f}|`$), one has $`\dot{\varphi }=6\lambda c_1c_2c_3|t|^{3s}e^{2\varphi }`$ by using (23). Integrating this equation with respect to $`t`$, we easily find that $`\dot{\varphi }=\mathrm{sign}(t)(3s+1)/(2|t|)`$. Combining this with eqs. (7)-(10) gives

$`s=2,c_4={\displaystyle \frac{5}{2}},c_1c_2c_3=\mathrm{sign}(t){\displaystyle \frac{5e^{2\varphi _0}}{12\lambda }}.`$ (26)
Since $`c_1c_2c_3<0`$ for $`t<0`$, the Universe is either contracting in all directions or expanding in two directions and contracting in one direction. When $`t>0`$, the Universe is either expanding in all directions or expanding in one direction and contracting in two directions. For the asymptotic solution (26), we have that $`\dot{\varphi }|t|^1`$ and $`f^{}|t|^5`$, in which case the determinant is given by $`D|t|^6\mathrm{}`$ for $`|t|\mathrm{}`$. In spite of this divergent behavior of the determinant, the solutions are nonsingular with $`p,q,r|t|^20`$ for $`|t|\mathrm{}`$.

We shall classify the cases where the solutions of eqs. (7)-(10) exhibit singular behavior. When the system passes through the determinant singularity ($`D=0`$), eq. (12) indicates that $`𝐱=(\dot{p},\dot{q},\dot{r},\ddot{\varphi })`$ diverge. This singularity appears in an anisotropic background where three expansion rates are multiple-valued functions of time . It is also a physical singurality where the curvature invariant, $`R^{\mu \nu \alpha \beta }R_{\mu \nu \alpha \beta }`$, diverges due to the divergence of the time derivative of the expansion rates. Near the determinant singularity, the expansion rates and the scalar field can be expanded as

$`h_i`$ $`=`$ $`h_{is}+h_{i1}\sqrt{|tt_s|}+h_{i2}(\sqrt{|tt_s|})^2+\mathrm{},`$ (27)
$`\varphi `$ $`=`$ $`\varphi _s+\varphi _1(\sqrt{|tt_s|})^2+\varphi _2(\sqrt{|tt_s|})^3+\mathrm{},`$ (28)
where $`h_i=p,q,r(i=1,2,3)`$, and $`t_s`$ is the time at singularity. This means that $`\dot{h}_i`$ and $`\ddot{\varphi }`$ diverge as $`tt_s`$, while $`p,q,r,\dot{\varphi }`$ are finite. This property is different from the ordinary kind of singularity where $`p,q,r`$ do not stay finite. The determinant singularity plays a crucial role in an anisotropic background.

The ordinary kind of singularities can be classified as

$`(\mathrm{i})`$ $`pp_0/(tt_s),qq_0,rr_0,`$ (29)
$`(\mathrm{ii})`$ $`pp_0/(tt_s),qp_0/(tt_s),rr_0,`$ (30)
$`(\mathrm{iii})`$ $`pp_0/(tt_s),qp_0/(tt_s),rp_0/(tt_s),`$ (31)
where $`p_0,q_0`$, and $`r_0`$ are constants with $`p_0>0`$. If the time-direction is futurewards ($`tt_s0`$) one has $`p\mathrm{}`$, while $`p\mathrm{}`$ for $`tt_s+0`$. The asymptotic forms of the determinant $`(\text{19})`$ depend upon the cases presented above.

In the case (i) with a plus sign in rhs of Eq. (21), the signs of $`q_0`$ and $`r_0`$ are the same and the asymptotic form of $`\dot{\varphi }`$ is given by $`\dot{\varphi }(q_0+r_0)/(3\lambda e^{2\varphi _s}q_0r_0)`$ for $`tt_s+0`$, with $`\varphi _s`$ being a constant. Then the determinant $`D`$ yields from eq. (19):

$`D`$ $``$ $`{\displaystyle \frac{4\lambda ^2e^{4\varphi _s}p_0^2}{3(tt_s)^2}}\left[\left(q_0{\displaystyle \frac{r_0}{2}}\right)^2+{\displaystyle \frac{3}{4}}r_0^2\right]\mathrm{}(\mathrm{for}q_0r_0>0).`$ (32)
In the minus sign of eq. (21), one has $`D\mathrm{}`$ for $`q_0r_0<0`$ and $`tt_s+0`$.

In the case (ii) with a plus sign of of eq. (21), one has $`\dot{\varphi }1/(3\lambda r_0e^{2\varphi _s})`$ for $`r_0>0`$ and $`tt_s+0`$. Then the determinant should asymptotically takes the form

$`D`$ $``$ $`{\displaystyle \frac{4\lambda ^2p_0^4e^{4\varphi _s}}{3(tt_s)^4}}\mathrm{}(\mathrm{for}r_0>0).`$ (33)
In numerical analysis we did not find this case for the parameter ranges and initial conditions we adopt. Typically ordinary singularities are dominated by the cases where one or three expansion rates tend to diverge. When one chooses the minus sign of eq. (21), the asymptotic behavior is $`D\mathrm{}`$ for $`r_0<0`$ and $`tt_s+0`$.

In the case (iii) with a plus sign of of eq. (21), we find $`\dot{\varphi }t/(\lambda p_0e^{2\varphi _s})`$ and

$`D0,`$ (34)
for $`tt_s+0`$. Strictly speaking this holds only for the isotropic case ($`\alpha =\beta =0`$) where all expansion rates are the same. In this case the solutions do not cross the determinant singularity, although they approach $`D=0`$ as $`tt_s+0`$. When small anisotropies are included, the trajectories can pass through $`D=0`$. For the minus sign of of eq. (21) the asymptotic behavior of the determinant is not described by eq. (34), as we will see later.

In another limit, $`tt_s0`$, the signs of diverging expansion rates in eqs. (29)-(31) are reversed, in which case the asymptotic forms of $`D`$ are altered. Nevertheless the determinant of the cases (i) and (ii) generally approaches the asymptotic value, $`D=\mathrm{}`$. We shall confirm this by numerical investigations in subsequent sections.

### A Plus sign of eq. (21)

We first analyze the case of the plus sign in eq. (21). When $`\varphi `$ is largely positive, the term $`f^{}=(\lambda /8)e^{2\varphi }`$ is negligible in (21), implying that $`\dot{\varphi }`$ is positive as long as $`pq+qr+rp>0`$. In this case $`\varphi `$ increases toward the future, which results in nonsingular asymptotic solutions with determinant, $`D2`$. When we go back to the past, there are two possibilities for the evolution of the determinant. One is the case where the solutions pass through the determinant singularity ($`D=0`$) and another is the one where $`D`$ goes toward infinity ($`D+\mathrm{}`$). Our numerical investigations suggest that the latter case does not occur for the plus sign in eq. (21), which means that solutions nonsingular in future ($`D2`$) meet the determinant singularity in past, irrespective of the initial values of $`\varphi `$.

In the left panel of Fig. III A (a) we plot the evolution of the expansion rates for the anisotropy parameters $`\alpha =0.05,\beta =0.05`$ with initial conditions $`h=0.16`$ and $`\varphi =0`$, corresponding to $`D1.7`$ at $`t=0`$. The determinant continues to grow until it approaches the finite value $`D=2`$ as $`t\mathrm{}`$, which indicates that this solution belongs to the nonsingular solution in future given by eq. (25). Note that for $`t>0`$ the Universe is expanding in all directions. When we solve the equations of motion pastwards ($`t<0`$), the solution meets the determinant singularity around $`t=0.82`$, thereby leading to the divergence of $`\dot{p},\dot{q},\dot{r}`$ and $`\ddot{\varphi }`$. If we introduce a new time parameter, $`\tau `$, defined by

$`\tau {\displaystyle \frac{dt}{D}},`$ (35)
it becomes possible to enter the region of the negative sign of $`D`$ by overpassing the determinant singularity . This does not mean that we can remove the singularity by coordinate transformations. The determinant singularity is a physical one where the divergence of the curvature invariant is unavoidable even in other coordinates. Then the solution turns back futurewards and the determinant begins to decrease rapidly toward $`D\mathrm{}`$ \[see the right panel of Fig. III A (a)\]. From Fig. III A (a) we find that this belongs to the class of the case (i) with $`pp_0,qq_0,rr_0/(tt_s)`$ and $`p_0q_0>0`$ \[see eq. (32)\]. The Universe is rapidly contracting in one direction ($`r\mathrm{}`$) for $`tt_s0`$. As claimed in ref. , these trajectories can be understood as the pair creation of two branches ($`D>0`$ and $`D<0`$) at the determinant singularity. For the nonsingular solutions in future, it is inevitable to cross the determinant singularity in past in an anisotropic background. Notice that in the isotropic case $`D`$ is always positive and decreases toward zero as $`t\mathrm{}`$. In the presence of small anisotropies, however, the solution reaches the determinant singularity at finite past and fall into the $`D=\mathrm{}`$ singularity as shown in Fig. III A (a).

When $`\varphi `$ is large, the allowed anisotropy parameters in $`(\alpha ,\beta )`$-plane lie inside a circle, given by

$`\alpha ^2+\beta ^2\stackrel{<}{}h^2,`$ (36)
which comes from the constraint (20). In Fig. III A (a) we show nature of singularities in past and future in the $`(\alpha ,\beta )`$-plane for $`h=0.16`$ and $`\varphi =0`$ at $`t=0`$. In this case the allowed region is approximately described by eq. (36). We find that all solutions inside this region exhibit the determinant singularity in past (region II in Fig. III A) while they are nonsingular ($`D2`$) in future (region I) .

With the decrease of $`\varphi `$, the allowed region gets larger as found by eq. (20). Nature of singularities becomes more complicated due to the appearance of ordinary singularities. If the system is close to the isotropic case ($`|\alpha |,|\beta |1`$), one has $`D2`$ for $`t\mathrm{}`$ and $`D0`$ at finite past as in the case of Fig. III A (a). However, larger anisotropies alter this picture. For example, we show the evolution of the expansion rates in the left panel of Fig. III A (b) for $`\alpha =0.15`$ and $`\beta =0.05`$ with $`h=0.16`$ and $`\varphi =1`$ at $`t=0`$. In this case $`p`$ is singular at finite past while $`q,r`$ are positive constants, which means that $`D\mathrm{}`$ from eq. (32). The solution exhibits the same kind of singularity ($`D\mathrm{}`$) in future. In this case we have numerically found the determinant is always negative ($`D=1.02`$ at $`t=0`$) and the solution never crosses $`D=0`$ in both past and future \[see the right panel of Fig. III A (b)\]. In Fig. III A (b) the asymptotic property of singularities is presented for $`h=0.16`$ and $`\varphi =1`$ at $`t=0`$. There exists the region III where the determinant is singular ($`D\mathrm{}`$) in both past and future. We also find that the region IV with $`D+\mathrm{}`$ appears in past when anisotropy parameters are large. This corresponds to asymptotic past nonsingular solutions given by eq. (26) where quadratic curvature corrections are dominant. It is important to separate this case from the ordinary singularity with $`D=\mathrm{}`$, although this classification was not done in ref. . From Fig. III A (b) the solutions with $`D=0`$ in future correspond to, in past, the determinant singularity or the regular solutions with $`D+\mathrm{}`$. Notice that the parameter range of future nonsingular solutions are smaller compared to the case of Fig. III A (a).

In Fig. III A (c) we show the density plot for $`h=0.16`$ and $`\varphi =2`$ at $`t=0`$. This suggests that future nonsingular trajectories are restricted to be very narrow near the isotropic point, $`\alpha =\beta =0`$. We also find that the determinant singularity in future corresponds to the regular solutions with $`D+\mathrm{}`$ in past. As one example we plot in Fig. III A (c) the evolution of $`p,q,r`$ for $`h=0.16`$ and $`\varphi =2`$ at $`t=0`$ with anisotropy parameters, $`\alpha =0.1`$, and $`\beta =0.1`$. In this case the past asymptotic solution is categorized in the nonsingular solution given by eq. (26). Fig. III A (c) shows that two expansion rates $`p,q`$ are positive while $`r`$ is negative for $`t<0`$, implying $`c_1c_2c_3<0`$ in eq. (26). The solution comes regularly from the asymptotic past with $`D`$ being decreased toward the future \[see the right panel of Fig. III A (c)\]. It crosses the determinant singularity around $`t=0.3`$, after which the determinant continues to decrease until the solution falls into the ordinary singularity with $`D\mathrm{}`$. Note that $`p`$ diverges as $`tt_s+0`$, while $`q`$ and $`r`$ approach positive constant values, in which case one has $`D\mathrm{}`$ by eq. (32). This trajectory can be regarded as the pair annihilation of two branches with $`D>0`$ and $`D<0`$. From Fig. III A (b) and (c), the solutions nonsingular in past ($`D+\mathrm{}`$) meet the determinant singularity in future. In this case the solutions do not go futurewards beyond the determinant singularity. In Fig. III A (c) the parameter regions with ordinary singularity in future ($`D=\mathrm{}`$) typically correspond to those with the same singularity in past, in which case the asymptotic behavior of the expansion rates are the similar as in Fig. III A (b). Notice that in Fig. III A (c) there exist some parameter ranges where the solutions meet the determinant singularity in both past and future. Although this case is very rare, it is still possible to cross the determinant singularity twice.

From Fig. III A we find that nonsingular solutions in future are of the determinant-type in past. This can be understood that the determinant evolves from $`D=2`$ to $`D=0`$, if we solve the equations of motions from asymptotic future to past \[see the right panel of Fig. III A (a)\]. We also analyzed other cases varying the values of $`h`$ and $`\varphi `$ for $`10^5<h<1`$ and $`15<\varphi <15`$. For a fixed value of $`h`$ there exists a minimal value of $`\varphi `$ which leads to nonsingular solutions in future \[see Fig. III A\]. The allowed range of the dilaton for future nonsingular solutions gets wider with the decrease of $`h`$ as shown in Fig. III A. In all cases analyzed in our numerical simulations, the past singularity for future nonsingular solutions corresponds to the determinant-type in an anisotropic background. We also found that trajectories with nonsingular past asymptotic meet the determinant singularity in future. In addition we showed that the ordinary singularity with (29) appears, in which case the determinant is divergent as $`D\mathrm{}`$. This is different from the case $`D+\mathrm{}`$ where the solutions are nonsingular in asymptotic past. In the next subsection we will analyze how the behavior of the determinant is altered for the different sign in eq. (21).

### B Minus sign of eq. (21)

In the previous subsection solutions nonsingular in future meet the determinant singularity in past, due to the fact that the determinant evolves from the asymptotic value $`D=2`$ toward the determinant singularity. For the minus sign of eq. (21), the determinant can increase toward the past, implying that the trajectories may not cross the determinant singularity. In fact we will show that the past singularity is not necessarily of the determinant-type.

Let us first analyze the asymptotic behavior of the determinant when all expansion rates diverge as eq. $`(\text{31})`$. In eq. (21) when the $`2(pq+qr+rp)`$ term is positive and much larger than $`(24pqrf^{})^2`$, one has $`\dot{\varphi }\sqrt{2(pq+qr+rp)}\sqrt{6}p_0/(tt_s)`$. Therefore the dilaton grows as $`\varphi \sqrt{6}p_0\mathrm{ln}(tt_s)`$ for $`tt_s+0`$, in which case the asymptotic form of the determinant yields

$`D2`$ $`+`$ $`6\sqrt{6}\lambda p_0^2(tt_s)^{2(\sqrt{6}p_01)}+39\lambda ^2p_0^4(tt_s)^{4(\sqrt{6}p_01)}`$ (37)
$`+`$ $`18\sqrt{6}\lambda ^3p_0^6(tt_s)^{6(\sqrt{6}p_01)}+18\lambda ^4p_0^6(tt_s)^{8(\sqrt{6}p_01)}.`$ (38)
This indicates that the determinant grows infinitely ($`D+\mathrm{}`$) for $`p_0<1/\sqrt{6}`$, while it decreases toward a finite value $`D=2`$ for $`p_0>1/\sqrt{6}`$. The latter case corresponds to the one where the average expansion rates are large initially.

When the condition $`(24pqrf^{})^2|2(pq+qr+rp)|`$ is satisfied, the asymptotic form of $`\dot{\varphi }`$ is given by $`\dot{\varphi }6\lambda p_0^3e^{2\varphi }/(tt_s)^3`$, which yields $`e^{2\varphi }6\lambda p_0^3/(tt_s)^2`$. Then one has the following positive finite value of the determinant in the limit of $`tt_s+0`$:

$`D2+{\displaystyle \frac{13}{12p_0^2}}+{\displaystyle \frac{7}{36p_0^4}}+{\displaystyle \frac{5}{432p_0^6}}.`$ (39)
Whether the asymptotic form is given by (38) or (39) depends on the initial values of the expansion rates and the dilaton. In the latter case $`\varphi `$ is generally small so that the condition $`(24pqrf^{})^22(pq+qr+rp)`$ is fulfilled. In both cases we can expect that the trajectories do not cross the determinant singularity in past.

We show in Fig. III B the density plot of the asymptotic behavior of the determinant for $`h=0.16`$ at $`t=0`$. When $`\varphi =0`$ at $`t=0`$ all trajectories satisfying the constraint (20) are nonsingular in future with $`D=2`$. In this case the past solutions are dominated by the singularity with asymptotic positive determinant as found from the left panel of Fig. III B (a). As one example we plot in Fig. III B (a) the evolution of $`p`$, $`q`$, $`r`$, $`\varphi `$, and $`D`$ for $`\varphi =0`$ at $`t=0`$ with anisotropy parameters $`\alpha =0.05`$ and $`\beta =0.05`$. This belongs to the class (38) with $`\alpha <1/\sqrt{6}`$ where all expansion rates diverge with determinant $`D+\mathrm{}`$. The determinant is always positive and continues to grow pastwards. The left panel of Fig. III B (a) shows that there exist some past solutions which meet the determinant singularity when anisotropies parameters are large. In this case the determinant grows until some moment of time pastwards, after which it begins to decrease toward the determinant singularity. This behavior can be understood that large anisotropies prevent all expansion rates from evolving almost similarly as eq. (31). However this region is typically small for the positive values of $`\varphi `$ at $`t=0`$.

With the decrease of the initial $`\varphi `$, the constraint (20) gives wider allowed parameter ranges in $`(\alpha ,\beta )`$-plane. Let us consider the case $`\varphi =1`$ and $`h=0.16`$ at $`t=0`$. As is found from Fig. III B (b) we have additional regions with the $`D=\mathrm{}`$ singularity in both past and future. Nature of singularities around the isotropic point is similar to the $`\varphi =0`$ case explained above. For the solutions nonsingular in future, the past singularity is either the type of $`D>0`$ or $`D=0`$. Fig. III B (b) is the latter case where the trajectory comes regularly from the asymptotic future ($`D=2`$) and meets the determinant singularity in past. This evolution is similar to the case of large anisotropy parameters in Fig. III B (a). From Fig. III B (b) the solutions singular in future ($`D=\mathrm{}`$) correspond to, in past, the singularity with $`D=\mathrm{}`$ or the determinant singularity. The difference between two cases is whether the determinant is always negative \[like Fig. III A (b)\] or it passes through $`D=0`$.

The situation becomes somewhat different with $`\varphi `$ being decreased further. Fig. III B (c) indicates that nonsingular trajectories with $`D+\mathrm{}`$ appear in future around the isotropic point for $`\varphi \stackrel{<}{}2`$ at $`t=0`$. This is the case all expansion rates are finite as described by eq. (26). However these solutions are singular in past with all expansion rates being infinite. In the region V shown in Fig. III B (c), the condition $`(24pqrf^{})^2>|2(pq+qr+rp)|`$ is typically satisfied at $`t=0`$. Therefore the determinant tends to approach the finite value (39) in past. We show in Fig. III B (c) the evolution of the system for $`\alpha =0.02`$, $`\beta =0.02`$, $`\varphi =2`$, and $`h=0.16`$. While the evolution of the expansion rates in Fig. III B (c) looks similar as in Fig. III B (a), the behavior of the determinant is different in both past and future. In Fig. III B (c) the determinant decreases from infinity to a finite positive value toward the past. We have numerically found that two terms $`(24pqrf^{})^2`$ and $`2(pq+qr+rp)`$ become comparable during the evolution. Therefore the asymptotic value of $`D`$ is somewhat different from eq. (39), but it is still a finite positive value. We find from Fig. III B (c) that future nonsingular solutions with either $`D=2`$ or $`D=+\mathrm{}`$ correspond to the singularity (31) with positive determinant in past. Figs. III B (b) and (c) also suggest that as $`\varphi `$ decreases the ordinary singularity with $`D\mathrm{}`$ appears in both past and future for large anisotropy parameters,

Compared to the case of the previous subsection, the range of nonsingular solutions in future is not so narrow due to the presence of nonsingular trajectories with $`D+\mathrm{}`$. However this solution is not appropriate to lead to our present Universe due to the dominance of the quadratic curvature term. In addition these solutions typically approach the ordinary singularity (31) with all expansion rates infinite in past. This property is different from the plus sign of eq. (21) where the past singularity for solutions nonsingular in future corresponds to the determinant-type. In both cases we have found that nonsingular cosmological solutions in both past and future do not exist for the dilaton-driven case for wide ranges of the parameter space ($`10^5<h<1`$ and $`15<\varphi <15`$). However the situation is changed in the modulus-driven case as we will analyze in the next section.

## IV Modulus-driven case

In the modulus-driven case the function $`f(\varphi )`$ in eq. (1) is expressed as

$`f(\varphi )={\displaystyle \frac{1}{16}}\delta \xi (\varphi ),`$ (40)
where the coefficient $`\delta `$ is determined by the 4-dimensional trace anomaly of the $`N=2`$ sector. Here the function $`\xi (\varphi )`$ is defined by

$`\xi (\varphi )=\mathrm{ln}\left[2e^\varphi \eta ^4(ie^\varphi )\right],`$ (41)
where $`\eta (ie^\varphi )`$ is the Dedekind $`\eta `$ function. Since $`\xi (\varphi )`$ is well approximated as $`\xi (\varphi )(2\pi /3)\mathrm{cosh}\varphi `$ , the function $`f`$ takes the form

$`f(\varphi ){\displaystyle \frac{\pi \delta }{48}}\left(e^\varphi +e^\varphi \right).`$ (42)
When $`\varphi `$ is largely negative ($`|\varphi |1`$), eq. (42) reduces to the form (22) by setting $`\delta =3\lambda /\pi `$. Therefore when $`\delta >0`$ solutions nonsingular in both past and future do not exist in this case. However nonsingular cosmological solutions have been found for negative values of $`\delta `$ . Hereafter we shall focus on the negative $`\delta `$ case (setting $`\delta =1`$ for simplicity) using the function $`f(\varphi )`$ given by (42).

When the Gauss-Bonnet term dominates in eqs. (7)-(10), the asymptotic solution is similarly given as eq. (26), i.e.,

$`p=c_1|t|^2,q=c_2|t|^2,r=c_3|t|^2,\varphi =\varphi _0\pm 5\mathrm{l}\mathrm{n}|t|,c_1c_2c_3=\mathrm{sign}(t){\displaystyle \frac{5e^{2\varphi _0}}{4\pi \delta }},`$ (43)
The past asymptotic solutions correspond to $`c_1c_2c_3>0`$ for negative $`\delta `$. This sign is different from the dilaton-driven case with positive $`\lambda `$ which we already analyzed in the previous section. Note that the determinant is divergent ($`D+\mathrm{}`$) for $`|t|\mathrm{}`$.

For another asymptotic solution where the effect of the Gauss-Bonnet term is negligible, the evolution of the background is given by

$`p=c_1|t|^1,q=c_2|t|^1,r=c_3|t|^1,\varphi =\varphi _0+c_4\mathrm{ln}|t|,`$ (44)
$`c_1+c_2+c_3=\mathrm{sign}(t),c_1^2+c_2^2+c_3^2+c_4^2=1,c_4^2=2(c_1c_2+c_2c_3+c_3c_1).`$ (45)
Therefore one has $`f^{}|t|^{|c_4|}`$, $`\dot{\varphi }|t|^1`$, and $`p,q,r|t|^1`$ with $`|c_4|<1`$ in eq. (19), in which case the solutions are nonsingular with an asymptotic value of the determinant, $`D=2`$.

For the nonsingular cosmological solutions found in ref. , $`\dot{\varphi }`$ does not change its sign . Due to the symmetric structure of the function (42) with respect to $`\varphi =0`$, we will consider the positive $`\dot{\varphi }`$ case where the modulus continues to grow from past to future. For negative $`\dot{\varphi }`$ the analysis is essentially the same by changing $`\varphi `$ to $`\varphi `$. When the solutions are singular, they meet the determinant singularity \[see eqs. (27) and (28)\] or the ordinary singularity \[see eqs. (29)-(31)\]. Let us consider the asymptotic behavior of the determinant for the ordinary singularity. In the cases (29)-(31), possible asymptotic behavior of the determinant can be summarized as

$`\mathrm{Case}(\mathrm{i})`$ $`D{\displaystyle \frac{256f_{s}^{}{}_{}{}^{2}p_0^2}{3(tt_s)^2}}\left[\left(q_0{\displaystyle \frac{r_0}{2}}\right)^2+{\displaystyle \frac{3}{4}}r_0^2\right]\mathrm{},\mathrm{with}\dot{\varphi }{\displaystyle \frac{1}{24f_s^{}}}{\displaystyle \frac{q_0+r_0}{q_0r_0}}\mathrm{const},`$ (46)
$`\mathrm{Case}(\mathrm{ii})`$ $`D{\displaystyle \frac{256f_{s}^{}{}_{}{}^{2}p_0^4}{3(tt_s)^4}}\mathrm{},\mathrm{with}\dot{\varphi }{\displaystyle \frac{1}{24r_0f_s^{}}}\mathrm{const},`$ (47)
$`\mathrm{Case}(\mathrm{iii})`$ $`D0,\dot{\varphi }{\displaystyle \frac{tt_s}{8p_0f_s^{}}},`$ (48)
where $`f_s^{}(\pi \delta /48)(e^{\varphi _s}e^{\varphi _s})`$ with $`\varphi _s`$ being a constant. Although the case (ii) is rare, we have numerically checked that this asymptotic solution certainly exists for the modulus-driven case. The asymptotic (iii) corresponds to the isotropic case, as is similar to the dilaton-driven cosmology. It is easy to verify that this trajectory is possible for positive $`\varphi `$ but impossible for negative $`\varphi `$. In what follows we shall analyze how the evolution of the expansion rates and the determinant are different from that of the dilaton-driven case paying particular attention for nonsingular trajectories.

Nonsingular asymptotic solutions for $`\dot{\varphi }>0`$ can be described as $`\varphi =\varphi _05\mathrm{l}\mathrm{n}|t|`$ in past \[see eq. (43)\] and $`\varphi =\varphi _0+c_4\mathrm{ln}|t|`$ with $`0<c_4<1`$ in future \[see eq. (45)\]. Therefore $`\varphi `$ continues to grow from asymptotic past starting from large negative values of $`\varphi `$ toward the future. Such examples are plotted in Fig. IV (a) and (b). In these cases the determinant evolves from $`D=+\mathrm{}`$ (past) to $`D=2`$ (future) without passing through the determinant singularity. In the dilaton-driven case when the past trajectories are nonsingular with $`D=+\mathrm{}`$ they inevitably meet the determinant singularity in future \[see Fig. III A (b) and (c)\]. This is mainly due to the fact that $`f^{}`$ is negative in the dilaton case while its sign is different in the modulus case for $`\delta <0`$ and $`\varphi <0`$. Therefore the determinant (19) is dominated by positive terms in the modulus case, which provides a way not to pass through the determinant singularity. Namely negative $`\delta `$ is crucial for the existence of nonsingular solutions.

In Fig. IV we show past and future asymptotic properties in three different cases. Note that we have defined $`\stackrel{~}{\alpha }\alpha /h`$ and $`\stackrel{~}{\beta }\beta /h`$ in order to compare the cases where the average expansion rate $`h`$ is changed. When $`h=0.05`$ and $`\varphi =5`$, nonsingular solutions in past ($`D=+\mathrm{}`$) are not singular in future with determinant, $`D2`$ \[see Fig. IV (a)\]. The Universe exhibits superinflation with growing expansion rates until the graceful exit around $`t=0`$ \[see Fig. IV (a)\]. Notice that we have $`pqr>0`$ in asymptotic past, as predicted by eq. (43) for negative $`\delta `$. The expansion rates begin to decrease after the graceful exit, whose asymptotic solutions in future are given by eq. (45). In Fig. IV (a) we find that the future solution corresponds to the Kasner-type where the Universe is contracting in one direction.

Fig. IV (a) indicates that some trajectories which are nonsingular in future cross the determinant singularity in past. The evolution of the background is similar to Fig. III A (a) which we already analyzed in the dilaton case. We also find from Fig. IV (a) that when anisotropy parameters are large the solutions meet the ordinary singularity with $`D=\mathrm{}`$ in both past and future. This is the case of (29) or (30) where at least one expansion rate diverges as plotted in Fig. III A (b). Although in Fig. IV past nonsingular solutions ($`D+\mathrm{}`$) do not meet the determinant singularity in future, we have checked that this singular behavior occurs for smaller values of $`\varphi `$ as shown in Fig. III A (c). These results imply that the property of singularities is similar to the dilaton-driven case described in Sec. III A.

For larger initial values of $`\varphi `$, the allowed region can be approximately described as $`\stackrel{~}{\alpha }^2+\stackrel{~}{\beta }^2\stackrel{<}{}1`$. When $`h=0.05`$ and $`\varphi =2`$ shown in Fig. IV (b), all future asymptotic solutions are nonsingular with determinant, $`D2`$. However we find from the left panel of Fig. IV (b) that the region of past nonsingular solutions gets smaller relative to Fig. IV (a). For large anisotropy parameters, the determinant tends to decrease pastwards from asymptotic future with $`D=2`$, thereby resulting in the singularity with $`D=0`$. This is not the case for the nonsingular trajectories with small anisotropy parameters. One example is plotted in Fig. IV (b). Although the determinant decreases from $`t=+\mathrm{}`$ to $`t6`$, it begins to grow toward the asymptotic past before crossing the determinant singularity. From Fig. IV (b) we find that the trajectory comes regularly from asymptotic past and connects to the Friedmann-type branch where the Universe is expanding in all directions. When anisotropy parameters are small and belong to the region (6), the future solution is of the Friedmann-type.

When the average expansion rate $`h`$ gets larger, the region of nonsingular cosmological solutions becomes smaller. We plot in Fig. IV (c) the density plot of nature of the determinant for $`h=0.1`$ and $`\varphi =2`$. The past solutions are dominated by the determinant singularity, whereas the future solutions are nonsingular as is similar to Fig. IV (b). In this case nonsingular trajectories in both past and future exist in only small parameter ranges around the isotropic point.

It is also worth mentioning that the ordinary singularity (31) with asymptotic determinant (38) and (39) for the minus sign of eq. (21) does not appear in the present case. When $`\varphi `$ increases toward the past, $`|f^{}|`$ also grows in the modulus case. Therefore the condition $`(24pqrf^{})^2|2(pq+qr+rp)|`$ is not satisfied, implying that case of (38) does not occur. When $`(24pqrf^{})^2|2(pq+qr+rp)|`$, there are two possibilities for the asymptotic form of $`\dot{\varphi }`$, one of which is $`\dot{\varphi }\pi \delta (e^\varphi e^\varphi )p_0^3/(tt_s)^3`$ and another is $`\dot{\varphi }C(tt_s)`$ where $`C`$ is a constant. In the former case it is easy to show that asymptotic solutions do not exist by integrating $`\dot{\varphi }\pi \delta (e^\varphi e^\varphi )p_0^3/(tt_s)^3`$ \[Note that this is possible in the dilaton-driven case with determinant (39)\]. In the latter case the determinant approaches $`D=0`$ as in eq. (48).

We have done numerical simulations for other cases varying the values of $`h`$ and $`\varphi `$. When $`h`$ is large, $`\varphi `$ is required to be small for the existence of nonsingular cosmological solutions. This property is found in Fig. IV where we plot the regions of nonsingular and singular solutions in the $`(h,\varphi )`$-plane. Nonsingular trajectories come regularly from asymptotic past with $`D=+\mathrm{}`$ where the quadratic curvature term is dominant \[see eq. (43)\], and smoothly connect another nonsingular branch with $`D=2`$ where the Gauss-Bonnet term is negligible \[see eq. (45)\]. For singular solutions nature of singularities is found to be similar to that of the dilaton-driven case discussed in Sec. III A.

## V Summary and discussions

In this paper we have analyzed past and future asymptotic regimes in Bianchi type-I string-inspired cosmological models in the presence of a Gauss-Bonnet curvature invariant. If the loop correction is not taken into account, one has no-go results that the initial big-bang singularity can not be avoided. The Gauss-Bonnet term allows the existence of nonsingular cosmological solutions, depending on the theories we adopt. We investigated two gravity theories, viz., dilaton- and modulus-driven cosmologies. In the former case the dynamics appears to depend significantly on the sign of rhs of eq. (21). Hence, we treated possible signs in eq. (21) separately and constructed three pictures of cosmological evolution:

* Dilaton-driven cosmology with plus sign in eq. (21)
* Dilaton-driven cosmology with minus sign in eq. (21)
* Modulus-driven cosmology with $`\delta <0`$ .

As the quadratic curvature corrections may provide violations of strong and week energy conditions , the Bianchi type-I Universe can recollapse in high-curvature regime. This fact gives us a variety of possible types of trajectories with different past and future asymptotics. Some trajectories can not leave high-curvature regime, some reaches the low-curvature future attractor.

In the dilaton-driven cosmology the 4-dimensional string coupling, $`\lambda `$, is required to be positive. This forbiddens the existence of nonsingular cosmological solutions even when the Gauss-Bonnet term is taken into account. In the case (a) nonsingular trajectories in future correspond to the low-curvature solutions with determinant $`D2`$ where the Gauss-Bonnet term is negligible. These trajectories meet the singularity where the determinant of the system vanishes at finite past. At this determinant singularity, $`\dot{p},\dot{q},\dot{r},\ddot{\varphi }`$ diverge in eq. (12) while $`p,q,r,\dot{\varphi }`$ stay finite. This kind of singularities restricts the presence of nonsingular solutions in an anisotropic background. There exist other kind of singularities (we call ordinary singularities) where at least one expansion rate diverges \[see eqs. (29)-(31)\]. When the singularity (29) or (30) appears at finite past or future, the determinant approaches $`D=\mathrm{}`$. We also find some trajectories which are nonsingular in asymptotic past with determinant, $`D+\mathrm{}`$. These solutions appear for large anisotropy parameters, whose existence, to our knowledge, was not discovered previously. In the case (a) these solutions are found to meet the determinant singularity in future (see Fig. III A). In the case (b) there appears the ordinary singularity of eq. (31) where all expansion rates diverge. In this case the determinant is divergent or approaches a positive constant value, depending on the initial conditions of $`\varphi `$ and $`h`$ \[see (38) and (39)\]. The solutions nonsingular in future typically correspond to this type of singularity or the determinant singularity in past as found by Fig. III B.

In the modulus-driven case the coupling, $`\delta `$, can take either positive or negative value in eq. (42). When $`\delta >0`$ solutions nonsingular in both past and future ar not found as is similar to the dilaton-driven case. For negative values of $`\delta `$, however, there exist nonsingular cosmological solutions where two branches of superinflation and decreasing curvature can be joined to each other. In this case the solutions come regularly from asymptotic past ($`D=+\mathrm{}`$) with determinant being decreased. The determinant approaches the future asymptotic value $`D=2`$ without crossing the determinant singularity. This is the main difference from the dilaton-driven case where the past nonsingular trajectories inevitably meet the determinant singularity in future. When the solutions are singular we find that nature of singularities is similar to the dilaton-driven case discussed in Sec. III A.

Our numerical investigations show that not all combinations of past and future asymptotics are possible. We pay particular attention to nonsingular past and low-curvature nonsingular future regimes (the latter should correspond to our present Universe). For these regimes we get the following results:

* For the dilaton-driven case trajectories with nonsingular past asymptotic meet the determinant singularity in future.
* For cases (a) and (c) trajectories with low-curvature future asymptotic meet the determinant singularity in past, when the past solutions are singular.

The first property tells us that the negative sign of the coupling constant is essentially important for constructing a purely nonsingular string cosmological models. This result is known for FRW Universe  and, hence, is still valid in anisotropic Bianchi I case, despite the existence of a past nonsingular regime in the dilaton-driven cosmology. The second property indicates that the determinant singularity may play a crucial role when we try to trace back in time evolution of the our present Universe in anisotropic background.

Recently string-inspired cosmological models which can avoid the big bang singularity have received much attention  together with the proposal of the Ekpyrotic universe . In those cases the quantum loop corrections or the higher-order derivatives play important roles to determine the dynamics before the graceful exit. It is certainly of interest to extend our analysis to more complicated models such as the multi-field case in the presence of the Gauss-Bonnet term. In fact while there exist nonsingular solutions in the single-field modulus-driven case considered in this work, density perturbations generated by the fluctuation of modulus exhibits blue-spectra with a spectral tilt $`n=10/3`$ . This contradicts with the observational supported flat spectra with $`n1`$. In the multi-field case, however, it may be possible to produce almost scale-invariant spectra if a light scalar field such as axion generates a flat isocurvature perturbations during superinflation  or if the axion is nonminimally coupled to the dilaton with some potential . We leave to future work to construct such nonsingular cosmological models which are consistent with observations.

## ACKNOWLEDGEMENTS

We thank S. Alexeyev and H. Yajima for useful discussions. The work of AT is supported by Russian Foundation for Basic Researches via grants Ns. 00-15-96699 and 02-02-16817. ST is thankful for financial support from the JSPS (No. 04942).


text:
# 1 Introduction

## 1 Introduction

Following the introduction of unification schemes for active galactic nuclei (AGN), it is now widely accepted that the central nucleus of an active galaxy has a Kerr black hole, a relativistic accretion disk, jets, a dusty torus and clouds emitting broad emission lines. Jets and accretion disks are strongly related through a symbiosis, and we postulate that there may also be a large scale symbiosis involving dust components. Since $`\gamma `$-rays are attenuated by photon-photon pair production collisions with lower energy photons from the disk, BLR, or torus, such a symbiosis could have an effect on the production and absorption of $`\gamma `$-rays in AGN as the BLR and torus are powered by the disk, and so may depend on the way in which the jets are fed by the accretion disk.

A $`\gamma `$-ray flare could occur, for example, when a significant amount of energy has been accumulated in the inner part of the accretion disk and is subsequently expelled into the jets. As a consequence, relativistic shocks may form in the jets and could account for the observed rapid and high amplitude $`\gamma `$-ray variability. The intense emission observed during a flare could be also due to an increase in the accretion rate, in the bulk Lorentz factor of the emitting material in the jet, or in the efficiency of particle acceleration. The production of $`\gamma `$-rays from relativistic jets depends also on the density of radiation fields since these photons can be scattered by relativistic electrons to GeV and TeV energies, the external radiation field giving additional target photons to the internal synchrotron radiation assumed in the synchrotron self-Compton (SSC) process . Mannheim  and Protheroe  have shown that the disk radiation can be relevant for $`\gamma `$-ray production if the main constituents of the jet are relativistic protons.

In this paper we shall concentrate on the role played by external radiation fields in the absorption of $`\gamma `$-rays by photon-photon pair production, and we shall calculate the pair production opacity for GeV–TeV $`\gamma `$-rays emitted in the jets of quasars and blazars. We will only consider external radiation fields existing in AGN, and will not discuss pair production by $`\gamma `$-rays interacting with internal radiation fields in the emission region as this has been discussed elsewhere . We shall also not address the problem of absorption by the infrared background radiation during propagation to Earth (e.g. refs. ).

We shall assume that the disk radiation, the BLR radiation, and the torus radiation is anisotropic, but is symmetric about the jet axis. In such a radiation field the reciprocal of the mean interaction length for photon-photon pair production by $`\gamma `$-rays traveling along the jet axis is given by

$$x_{\gamma \gamma }^1(E,z)=_{\epsilon _{\mathrm{min}}}^{\mathrm{}}𝑑\epsilon _{\mathrm{cos}(\theta _{\mathrm{max}})}^{\mathrm{cos}(\theta _{\mathrm{min}})}\frac{dn}{d\mathrm{\Omega }}(\epsilon ,\theta ,z)(1\mathrm{cos}\theta )\sigma _{\gamma \gamma }(s)2\pi d\mathrm{cos}\theta $$
(1)
where $`\theta `$ is the angle between the direction of the soft photons and the jet axis, $`\sigma _{\gamma \gamma }(s)`$ is the photon-photon cross section for pair production at center of momentum frame energy squared $`s=2Eϵ(1\mathrm{cos}\theta )`$, $`dn/d\mathrm{\Omega }`$ is the differential photon number density, and $`\mathrm{cos}(\theta _{\mathrm{max}})`$ and $`\mathrm{cos}(\theta _{\mathrm{min}})`$ are defined by the disk, BLR, or torus geometry involved or, in the case of $`\mathrm{cos}(\theta _{\mathrm{min}})`$, by the pair production threshold condition. We calculate the optical depth for absorption of GeV and TeV photons from the point of $`\gamma `$-ray emission to infinity, in the relevant external radiation fields.

In quasars the torus and the BLR are active components, as can be inferred from the large forest of emission lines observed, and from the thermal infrared spectrum, allowing good estimation of the diffuse radiation within the central parsec. In contrast, blazars, which comprise BL Lacs and Optically Violently Variable (OVV) quasars, are a special category of AGN having the jet aligned closely to the line of sight. The lack of emission lines and UV bump in blazars is usually assumed to be due to a low thermal activity at the centre, implying that the energy density of the external radiation is small in blazars compared to that in quasars. We shall investigate whether or not this is the case, and discuss how the $`\gamma `$-ray output could be different in quasars and blazars depending on the dominant radiation field.

In Section 2 we shall discuss the optical/UV emission from an accretion disk in symbiosis with the jet, and its effect on the BLR emission clouds (ionized by disk photons) and the torus (heating of dust). It is difficult to estimate the disk luminosity in blazars as has been done for some quasars with visible UV bumps . Blazars do not show thermal disk emission, and their spectra from radio to UV frequencies, sometimes to X-ray frequencies, are universally attributed to synchrotron emission from the relativistic jet. The nonthermal boosted radiation from the jet dwarfs any thermal emission, making the direct detection of BLR or torus activity, and hence direct measurements of external photon fields present within the parsec scale of blazars, impossible.

We assume that the accretion disk surrounding the black hole and the relativistic jets are related by a small-scale symbiosis , and that the energy available for the jet is the energy which would be dissipated between the last marginally stable circular orbit (in Kerr geometry), $`R_{\mathrm{ms}}1.23R_g`$ where $`R_g=GM/c^2`$ is the gravitational radius, and the radius of the base of the jet $`R_{\mathrm{jet}}>R_{\mathrm{ms}}`$. The total power of the jet depends on $`R_{\mathrm{jet}}`$. A thick jet base means that more energy is expelled into the jet, and therefore, less energy is dissipated in the disk. A direct consequence of this assumption is that an accretion disk in symbiosis with jets (ADJ) produces a disk luminosity smaller than that of a “standard”  relativistic accretion disk. We suggest that the ADJ model could describe low-luminosity accretion disks in some blazars, and may provide a bridge between the standard accretion disk assumed to exist in quasars and Fanaroff-Riley Class II (FR-II) sources  and the low-power advection dominated accretion disk which may be found in some FR-I sources.

In Section 3, assuming reasonable models for the size and distribution of BLR clouds, we calculate the energy density of the radiation field of the BLR, and the optical depth for absorption of GeV and TeV photons. Broad Line Regions surround accretion disks and may also be important sources of soft photons which could be inverse-Compton scattered by relativistic electrons in the jet to produce $`\gamma `$-rays in “external Compton models” .

In section 4, we discuss the infrared radiation from a dusty torus. Based on the assertion that FR-I and BL Lacs are related , and that there is strong evidence that some FR-I have dusty tori , we assume that blazars must also have dusty tori  and discuss what evidence there is to support this. Since the infrared emission produced by tori could inhibit the escape of TeV $`\gamma `$-rays due to photon-photon pair-production  we study how changes in disk activity can modify the geometry and the heating mechanism of the torus, and how this affects the optical depth for $`\gamma `$-$`\gamma `$ absorption in the torus infrared radiation.

## 2 Radiation from an accretion disk with jets

Relativistic outflows are very nicely described by a model of the small-scale symbiosis between the accretion disk and jets . Donea and Biermann  have used such a model for an accretion disk with jets starting at the inner region of a disk, and were able to reproduce the UV bump in quasars. By fitting observed spectra they derived upper limits for the radius of the base of the jet, $`R_{\mathrm{jet}}`$, and they found that $`R_{\mathrm{jet}}`$ cannot be too far from $`R_{\mathrm{ms}}`$. The jet base, or “footring”, is actually the thin layer between $`R_{\mathrm{ms}}`$ and radius $`R_{\mathrm{jet}}`$, the jet being approximated by a hollow cylinder of inner radius $`R_{\mathrm{ms}}`$ and outer radius $`R_{\mathrm{jet}}`$, such that $`R_{\mathrm{ms}}R_{\mathrm{jet}}R_{\mathrm{out}}`$, where $`R_{\mathrm{out}}`$ is the outer radius of the disk.

The ADJ model assumes that the gravitational potential energy available between $`R_{\mathrm{ms}}`$ and $`R_{\mathrm{jet}}`$ is the energy reservoir of the jets. The total power of the jets is strongly dependent on the accretion mass rate in the disk and on the size of its footring (we not include here the interaction between the black hole and jets which could also result in it putting a non-negligible fraction of its energy into the jets ). The total power of the jets is

$$Q_{\mathrm{jet}}=L_{\mathrm{disk}}L_{\mathrm{disk}}^{\mathrm{jet}}$$
(2)
where $`L_{\mathrm{disk}}^{\mathrm{jet}}`$ is the luminosity of the disk with jets and $`L_{\mathrm{disk}}`$ would be the luminosity of the disk if the conditions necessary to drive outflows were not met, i.e., the case of a standard relativistic disk.

A coupled jet-disk system must obey the laws of conservation of mass and angular momentum. We assume that the jet is fed with mass by the accretion disk, and that the flow of mass into the jet per unit of time $`\dot{M}_{\mathrm{jet}}=q_m\dot{M}`$ is a fraction $`q_m1`$ of the accretion mass rate into the disk $`\dot{M}=2\pi R\mathrm{\Sigma }u^R`$, where $`\mathrm{\Sigma }`$ is the surface density of mass in the disk and $`u^R`$ is the radial velocity of gas at given radius $`R`$. The equation for conservation of mass requires $`\dot{M}=\dot{M}_{\mathrm{disk}}+\dot{M}_{\mathrm{jet}}`$. Derivation of the relevant equations can be found in Donea and Biermann , who followed the standard method of calculating the emission spectrum from an accretion disk . If there is no angular momentum and mass loss into the jet the equations used for mass and angular momentum transport in the disk become the standard equations , with $`R_{\mathrm{jet}}R_{\mathrm{ms}}`$, $`Q_{\mathrm{jet}}0`$ and $`L_{\mathrm{disk}}^{\mathrm{jet}}L_{\mathrm{disk}}`$. In this paper, we shall adopt $`q_m=0.1`$ as used by Falcke et al.  in interpreting the radio-UV correlation in AGN.

The local physics at the inner radius of the disk (radius $`R_{\mathrm{jet}}`$) is directly related to the extraction of angular momentum from the in-falling gas, and so modifies the structure of the relativistic disk . From this one calculates the dissipation energy at radius $`R`$, $`D^{}(R)`$, which must be done numerically in the case of a Kerr black hole, and a detailed discussion of this would serve no useful purpose here. Instead, for the purpose of illustration, and for the sake of simplicity we give here only the relation for the simpler case of a Schwarzschild black hole:

$$D^{}(R)=\frac{3GM\dot{M}}{8\pi R^3}\left[(1q_m)(1q_m)\left(\frac{R_{\mathrm{jet}}}{R}\right)^{1/2}\right]$$
(3)
where $`q_m=\dot{M}_{\mathrm{jet}}/\dot{M}`$ and $`M`$ is the black hole mass; nevertheless, in the present paper we use $`D^{}(R)`$ for the disk of a Kerr black hole. A reasonable outer radius of the footring of the jet would be $`R_{\mathrm{jet}}10R_\mathrm{g}10^4M_8`$ pc where $`M_8=M/10^8M_{}`$.

The disk/jet symbiosis is reflected mainly in a modified photon spectrum from the inner region of the accretion disk where the jet is anchored. The important result is that the spectrum from a Kerr accretion disk is cut off at high frequencies, from extreme UV to soft X-rays. In Fig. 1 we plot disk luminosities for the ADJ model, versus the thickness of the footring of the jet for Kerr black holes with masses $`M=10^8M_{}`$ and $`M=10^9M_{}`$, and for different mass accretion rates given by $`\dot{m}=\dot{M}/\dot{M}_{\mathrm{edd}}`$, where $`M_{\mathrm{edd}}`$ is the Eddington accretion rate. The case $`R_{\mathrm{jet}}=R_{\mathrm{ms}}`$ corresponds to an accretion disk without jets (“standard” accretion disk) and $`Q_{\mathrm{jet}}=0`$. As can be seen, for thicker jet bases more energy is available to power jets (dotted curves show $`Q_{\mathrm{jet}}`$ increasing with $`R_{\mathrm{jet}}R_{\mathrm{ms}}`$) and less energy is radiated by the disk (solid curves show $`L_{\mathrm{disk}}`$ versus thickness). We note that a small variation of the geometry at the coupling between jet and disk at radii $`R_{\mathrm{jet}}`$ less than $`3R_g`$ would induce large variations in the power of the jet, possibly causing flare activity in blazars. This is because, in a Kerr metric, at small radii close to the black hole there is a large amount of gravitational potential energy available for dissipation into the jet.

Because of the jet/disk coupling, the ADJ model predicts a lower disk luminosity and a softer photon spectrum, and some of the implications of this for interactions of accelerated electrons or protons are discussed in a separate paper . As mentioned earlier, a flaring state could arise from enhanced AGN central activity whereby the jets get more energy from the disk. In this case, the disk changes from an ADJ with weak jets into an ADJ with a spectrum cut off at high frequencies but with correspondingly stronger jets. Since, in this case, the average energy of photons from an ADJ with strong jets (flaring) becomes lower than the average energy of photons from an ADJ with weak jets this will change the shape of the resulting $`\gamma `$-ray spectrum in models  where disk photons are inverse-Compton scattered by relativistic electrons in the jet.

Luminosities of disks in blazars are probably $`L_{\mathrm{disk}}<10^{46}`$ erg/s, while quasars typically have $`L_{\mathrm{disk}}=10^{46}`$$`10^{48}`$ erg/s. The ADJ model gives a simplified approach to the symbiosis between the disk and jets, and could describe low-luminosity accretion disks for quasars and some blazars. However, for some blazars the ADJ model could still give UV fluxes much higher than those observed, and in this case an ejection dominated accretion flow (EDAF) model, which is an ADJ model with a wind , may be better. An EDAF model has been applied to the central activity of Sgr A\* where it was shown that a wind plus jets extracts energy from the disk more efficiency. Since the EDAF model does not leave much energy to be dissipated in the accretion disk, the radiation field from the disk would become even less important for GeV and TeV $`\gamma `$-ray absorption. Alternatively, the jets could extract more energy from the accreting gas if the disk with a jets turns into an advection dominated accretion flow (ADAF) .

Fig. 2(a) shows the optical depths from $`z=0.01`$ pc to infinity for absorption of GeV-TeV photons produced in the jet and interacting with photons from an ADJ with $`q_m=0.1`$. We see that at 1 TeV photons are only absorbed for high black hole masses and mass accretion rates close to the Eddington limit. This could apply to quasars where the disk radiation field has a high density around the base of the jet. For low accretion mass rates with $`\dot{M}0.1M_{\mathrm{edd}}`$, and lower black hole masses, $`M10^8M_{}`$, the photon-photon optical depth is $`\tau _{\gamma \gamma }<1`$ at 1 TeV. We note that the disk luminosity corresponding to the case $`\dot{m}=1`$ for $`M_8=1`$ lies roughly at the boundary between blazars and quasars. In Fig. 2(b) we plot the optical depth divided by $`L_{46}L_{\mathrm{UV}}/10^{46}`$ erg s<sup>-1</sup>, and give the likely range for blazars and quasars using $`\dot{m}=1`$ for $`M_8=1`$ as the boundary. In blazars we may have ADJ models with accretion rates as low as $`\dot{m}=0.01`$–0.001 – at lower accretion rates an ADAF or EDAF model including jets may be more realistic. We use $`M_8=1`$ and $`\dot{m}=0.01`$ as the low luminosity boundary of blazars with an ADJ in Fig. 2(b). Similarly, we use a supermassive black hole ($`M_8=10`$) and accretion close to the Eddington limit ($`\dot{m}=1`$) for the upper luminosity boundary of quasars in Fig. 2(b). Of course, these boundaries are by no means rigid and there will be some overlap between the two populations, particularly when we plot $`\tau _{\gamma \gamma }/L_{46}`$.

The ADAF model seems to explain better very low luminosity blazars, but it cannot energize the BLR or heat efficiently dust in the torus; the lack of observed BLR and torus activity for some blazars has been sometimes interpreted as suggesting that blazars do not have tori or BLR. We consider that if blazars are to be included in the unification schemes they should also have BLR and dusty tori, but possibly with lower emission than in quasars, and we shall discuss this possibility later.

## 3 Photons from broad line region clouds

We are interested in obtaining the optical depth for absorption of GeV and TeV photons in the anisotropic radiation field of the BLR. For this we need to know the distribution of the BLR clouds around the jet axis. Broad line regions in AGN are not spatially resolved, but information about their size and the gas distribution can be obtained by reverberation mapping, from which it appears that BLR are stratified, and have typical sizes of the order of light-weeks in Seyfert galaxies and light-months in quasars .

Despite the fact that the majority of BL Lacs do not show any emission lines, there are observations suggesting that there is a weak broad $`H_\alpha `$ emission line in BL Lacertae itself , and even Mrk 421 seems to have a very low BLR luminosity, $`L_{\mathrm{BLR}}=1.5\times 10^{40}`$ erg/s . Interestingly, Celotti et al.  have suggested, based on the rapid TeV variability of Mrk 421, that its accretion rate is lower than $`10^210^3`$ of the Eddington rate. Hence, a less efficient ADJ or an ADAF may then power the BLR for this source. For Mrk 501, the BLR luminosity has been found to be $`L_{\mathrm{BLR}}=4\times 10^{41}`$ erg/s . In addition, Chiaberge et al.  have concluded that BLR and obscuring tori are closely linked in quasars, and both are present only in association with accretion.

In order to model the BLR, and hence its radiation field, we need to know its size and the emission line emissivity as a function of radius. The two-phase (hot/cold) model of the BLR predicts a relatively thin shell-like BLR . Higher ionization lines are found at smaller radii (e.g. for NGC 5548 ), and this may favour a thin shell model for the BLR. This provides a very simple model of the BLR in which the total BLR region luminosity $`L_{\mathrm{BLR}}`$ is emitted from the surface of a sphere of radius $`r_{\mathrm{BLR}}`$, i.e. $`L_{\mathrm{BLR}}/4\pi r_{\mathrm{BLR}}^2`$ per unit area.

Although it is not known if strong winds can exist in blazars, there are models which explain the formation of clouds from bloated stars exposed to winds from accretion disks , and shocks from jets could also contribute to cloud formation . On the other hand, Frometh and Melia  suggest that the unsteady, turbulent accretion flow onto the black hole is subjected to several disturbances capable of producing shocks, which would heat the gas and allow for formation of BLR clouds. They also suggest that circular clouds originate at smaller radii. Results from reverberation mapping  suggest that the BLR in Seyferts is extended with an outer radius at least 10 times larger than the inner radius, and Kaspi and Netzer  use models with a radial dependence of BLR cloud number density and radius to fit their data.

Clearly, there are several possible mechanisms for formation of BLR clouds and it is not known which one operates in blazars. It could be that all these mechanisms play a role in BLR cloud formation. We shall adopt two extreme cases for the distribution of BLR cloud emissivity: (i) a geometrically thin shell at radius $`r_{\mathrm{BLR}}`$, (ii) a geometrically thick shell extending from radius $`r_{\mathrm{in},\mathrm{BLR}}`$ to $`r_{\mathrm{out},\mathrm{BLR}}`$ with the radial dependence of line emissivity based on the best fit parameters of Kaspi and Netzer  used in their modelling of the BLR in NGC 5548 (number density of clouds $`n_{\mathrm{cl}}r^{1.5}`$, cloud radius $`R_{\mathrm{cl}}r^{0.6}`$). In both cases we shall assume that the line emission is isotropic, and that the ionization source (the accretion disk) is located at $`r=0`$ and has luminosity $`L_{\mathrm{UV}}`$. Accretion disk photons excite the BLR material, and emission lines such as $`H_\alpha `$, $`H_\beta `$, $`OIII`$, $`NII`$, etc., are produced. We simplify our problem by making the approximation that the entire BLR luminosity is emitted in the $`H_\alpha `$ (6563 Å) line.

The BLR is assumed to be optically thin with optical depth $`\tau _{\mathrm{BLR}}`$ such that the total BLR luminosity is $`L_{\mathrm{BLR}}\tau _{\mathrm{BLR}}L_{\mathrm{UV}}`$. Celotti et al.  analyzed a sample of quasars and found an average value of $`L_{\mathrm{BLR}}4\times 10^{45}`$ erg/s, and for a sample of 12 blazars (including Mrk 421 and Mkr 501) found an average value of $`L_{\mathrm{BLR}}2\times 10^{43}`$ erg/s. However, flat spectrum radio quasars and OVV quasars could have higher $`L_{\mathrm{BLR}}`$ as required by Blazejowski et al.  in order to explain the $`\gamma `$-ray spectra for these objects. In our calculations we use $`L_{\mathrm{UV}}=10^{46}`$ erg/s and $`\tau _{\mathrm{BLR}}=0.01`$, which are typical values used for fitting the broad emission lines with photoionization models (for $`\tau =0.01`$, $`L_{\mathrm{BLR}}=10^{44}`$ erg/s). For $`\gamma `$-ray attenuation by photon-photon pair production we need the photon number density per solid angle of BLR photons along the jet axis, $`dn/d\mathrm{\Omega }(\theta ,z)=I(\theta ,z)/c`$ where $`I(\theta ,z)`$ is the intensity. Calculation of the intensity is straightforward for distribution (i) where the BLR is a geometrically thin shell. In Fig. 3 we show the full angular dependence of the BLR $`H_\alpha `$ intensity at different distances $`z`$ along the jet for the case where the BLR is a thin shell at $`r_{\mathrm{BLR}}=0.01`$ pc.

For distribution (ii) we shall calculate the intensity by integrating the emissivity along lines of sight through the BLR sketched in Fig. 4. A fraction $`d\tau _{\mathrm{BLR}}=n_{\mathrm{cl}}\sigma _{\mathrm{cl}}dr`$ of the central luminosity would then be re-processed in $`r(r+dr)`$. Since this spherical shell has volume $`4\pi r^2dr`$ the emissivity (in erg s<sup>-1</sup> sr<sup>-1</sup> cm<sup>-3</sup>) of re-processed radiation inside the BLR at radius $`r`$ is

$$j(r)=\frac{Ld\tau _{\mathrm{BLR}}}{16\pi ^2r^2dr}=\frac{Ln_{\mathrm{cl}}\sigma _{\mathrm{cl}}}{16\pi ^2r^2}.$$
(4)
For $`n_{\mathrm{cl}}=n_0(r/r_{in,\mathrm{BLR}})^\alpha `$ and $`R_{\mathrm{cl}}=R_0(r/r_{in,\mathrm{BLR}})^\beta `$, where $`n_0`$ and $`R_0`$ is the number density and the radius of clouds at $`r_{\mathrm{in},\mathrm{BLR}}`$, respectively, the optical depth of the entire BLR to central UV radiation is

$$\tau _{\mathrm{BLR}}=_{r_{\mathrm{in},\mathrm{BLR}}}^{r_{\mathrm{out},\mathrm{BLR}}}n_0\pi R_0^2\left(\frac{r}{r_{in,\mathrm{BLR}}}\right)^{\alpha +2\beta }𝑑r$$
(5)
from which we obtain the constant $`n_0R_0^2`$ for any assumed value of $`\tau _{\mathrm{BLR}}`$. Note that we shall use the power-law exponents preferred by Kaspi and Netzer , $`\alpha =1`$ and $`\beta =1.5`$. We treat separately the three possible emission region locations shown in Fig. 4: (a) blob 1 inside the cavity inside the BLR, (b) blob 2 located within the BLR region, and (c) blob 3 is outside the BLR. The intensity of radiation at angle $`\theta `$ to the jet axis at a position $`z`$ within the BLR region (e.g. blob 2 in Fig. 4) is given by

$$I(z,\theta )=_0^{l_{min1}}j(r)𝑑l+_{l_{min2}}^{l_{\mathrm{max}}}j(r)𝑑l,$$
(6)
where

$$r^2=z^22zl\mathrm{cos}\theta +l^2,$$
(7)
$$l_{\mathrm{max}}=z\mathrm{cos}\theta +z\left(\mathrm{cos}^2\theta +\left(\frac{r_{out,BLR}}{z}\right)^21\right)^{1/2},$$
(8)
and

$$l_{\mathrm{min1},2}=z\mathrm{cos}\theta \pm z\left(\mathrm{cos}^2\theta +\left(\frac{r_{in,BLR}}{z}\right)^21\right)^{1/2}.$$
(9)
The total (angle-integrated) energy density of BLR is:

$$u(z)=\frac{I(z,\theta )}{c}𝑑\mathrm{\Omega },$$
(10)
and is plotted in Fig. 5 as a function of distance along the jet. The energy density of BLR photons clearly depends also on the distribution of BLR clouds. At $`zr_{\mathrm{BLR}}`$ (distribution i) and $`zr_{\mathrm{out},\mathrm{BLR}}`$ (distribution ii) the BLR photon energy density approaches the limit $`u_{\mathrm{BLR}}L_{\mathrm{BLR}}/4\pi cz^2`$ (dashed line in Fig. 5). For distribution (i), where we have a geometrically thin spherical shell for the BLR, the energy density peaks sharply at $`r_{\mathrm{BLR}}`$ as found also by Ghisellini and Madau . In the case of distribution (ii) the energy density peaks near $`r_{\mathrm{in},\mathrm{BLR}}`$, but the curves are smoother as a result of the finite thickness of the BLR. Both cases are idealizations and there is probably a spread in cloud size at any given radius. From Fig. 5 we see that it is very important to know the right distribution of BLR clouds since there can be a difference of more than two orders of magnitude between the energy density of photons calculated for the same BLR luminosity but different BLR cloud distributions.

In Fig. 6(a) we show the optical depth from various emission region heights, $`z_0`$, to infinity for absorption of 100 GeV–100 TeV $`\gamma `$-rays in the radiation field of the BLR for distribution (i) where the BLR clouds are located in a geometrically thin shell at radius 0.01 pc. The optical depth is plotted for a BLR luminosity of $`10^{44}`$ erg/s, but can be linearly scaled for other BLR luminosities. The anisotropy of the BLR radiation field is reflected in the shapes of the curves as one goes to higher $`\gamma `$-ray energies. Fig. 6(b) gives the optical depth for cloud distribution (ii), with the inner radius at 0.01 pc (the same as $`r_{\mathrm{BLR}}`$ in part a) and outer radius 0.4 pc. In this case, while at $`z=0.01`$ pc the optical depth is less than for distribution (i), $`\tau _{\gamma \gamma }`$ drops off more slowly with $`z_0`$, such that at only slightly larger distances the optical depth becomes orders of magnitude higher than for distribution (i) until $`r_{\mathrm{out},\mathrm{BLR}}^b`$ is reached. This is in part because the possibility of near head-on collisions of $`\gamma `$-rays with BLR photons for $`z<r_{\mathrm{out},\mathrm{BLR}}^b`$ for distribution (ii) – large angles between directions of interacting photons means more collisions will be above the pair production threshold.

Fig. 7 shows the $`\gamma `$-$`\gamma `$ opacity for 10 TeV photons from $`z=z_0`$ to infinity vs. $`L_{\mathrm{BLR}}`$ for the two cases shown in Fig. 6. For blazars with a very low BLR luminosity, 10 TeV photons can probably escape through the BLR. However, in quasars BLR typically have $`L_{\mathrm{BLR}}10^{44}`$ erg/s, and so absorption of 10 TeV $`\gamma `$-rays can be quite important if the emission region is inside the BLR (distribution i) or is close to $`r_{\mathrm{in},\mathrm{BLR}}`$ (distribution ii).

We comment here on a very intriguing problem, namely the absence of a Ly absorption edge in quasars . It was expected that BLR clouds situated along the line of sight should produce a sharp Ly absorption edge in more than 30% of quasars . However, Maiolino et al.  have noted that none of the 101 quasars in the sample analyzed by Zheng et al.  show this feature. There are two possible explanations for the missing edge: (a) BLR clouds have a sufficiently large component of their random velocity along our line of sight to result in a broadening of the Ly edge and thereby making it difficult to observe, or (b) the BLR does not have a spherically symmetric distribution around the ionizing accretion disk, such that there are no BLR clouds intersecting our line of sight to the disk, and Maiolino et al.  suggest a flattened geometry could do this. Such a flattening in the cloud distribution could produce weaker densities of BLR photons and so would impede less the escape of TeV photons from the central region. However, we adopt a spherical distribution of clouds, and instead interpret the absence of an observed Ly absorption edge in 101 quasars to suggest that the BLR clouds must have a covering factor $`f_{\mathrm{cov}}0.03`$ or less (95% upper limit; for a covering factor 0.01 the probability of a BLR cloud being in the line of sight to the disk in none of 101 observations is 36%, quite consistent with the observations). For a given covering factor, $`f_{\mathrm{cov}}`$, the fraction of the disk luminosity intercepted by the BLR must then be less than $`f_{\mathrm{cov}}`$, and justifies our use of $`\tau _{\mathrm{BLR}}=0.01`$. Nevertheless, $`\tau _{\mathrm{BLR}}0.1`$ have also been widely used, and would lead to photon-photon pair production opacities a factor of 10 higher than those presented here.

## 4 Radiation from a dusty torus

Chiaberge et al.  have suggested that tori do not exist in FR-I sources, or they must have a very thin geometry. This conclusion results from the correlation found between optical and radio cores of some FR-I sources. However, the sample used in this analysis represents only just under half of the FR-I galaxies tabled by Zirbel and Baum , and one needs additional observations of the polarized optical spectrum to reach a firm conclusion about this . Nevertheless, these results could be interpreted as suggesting that perhaps there is a correlation between a flat infrared torus and a flat distribution of BLR clouds as suggested by Maiolino . We do not know what mechanism, other than accretion, could cause such a flattening of the BLR and the torus. However, such a flattening might be understood in terms of a non-homogeneous spherical shell of BLR or dust in which the polar caps above and below the disk plane are much more diluted than the equatorial belt. Zier and Biermann  have shown that the caps of a spherical layer of dust could be diluted so much, that even a dense torus would have a doughnut shape.

We have already discussed some aspects related to the infrared radiation from ubiquitous dusty tori in quasars. Unification schemes for AGN generally involve the following components: Kerr black hole, a relativistic accretion disk, jets, broad line clouds and a dusty torus. For blazars the situation appears at first sight to be a little different, and we wonder why that might be. There are two issues which should be considered in connection with dusty tori in blazars. Firstly, since there is no direct evidence of any thermal emission from tori we should consider what indirect evidence (at any wavelength) there is about the existence of dust in the centres of host galaxies. Secondly, we should consider what kind of torus geometry would fit observations made at different wavelengths best.

Since unification schemes propose that FR-I and BL Lacs are related, FR-I sources being the mis-oriented counterparts of BL Lacs , we have searched the literature for information about dust structures associated with blazars/BL Lacs/FR-I. We note that Blazejowski et al.  required an external IR photon field from a torus scattered by relativistic jet electrons, in addition to those in the SSC model, in order to fit $`\gamma `$-ray spectra of OVV quasars observed by EGRET. A key to the existence of tori in blazars could come from establishing a link between the torus and the BLR. For example, when there is no detection of direct optical emission, the interpretation of spectro-polarimetric data on FR-I objects with strong evidence for infrared obscuration ( and references therein) could suggest the existence of a BLR hidden by a thick torus. In this case, the observed polarized broad lines would be the result of scattering into the line of sight by free electrons in zones whose shape and orientation is determined by the torus’ inner geometry. Hence, the detection of free electron scattering regions should be a good diagnostic of torus geometry. Falcke et al.  proposed that the opening angle of the torus might play a critical role in the FR-I and FR-II dichotomy – a closed torus covering a large fraction of $`4\pi `$ steradians, as seen from the black hole, would obscure the internal activity of FR-I objects.

de Koff et al.  have discussed how the properties of the dust depend upon the radio properties of the object, and found that FR-II, which are powerful radio galaxies, have a rather clumpy dust distribution, and they suggest that this translates into a large opening angle of the torus. A less powerful radio jet (such as those in FR-I) would only slowly disperse the dust, and this type of obscuration might contribute to the deceleration of FR-I jets by entraining the material from the torus.

We consider first the FR-I galaxy Centaurus A. Alexander et al.  explain the IR spectrum of Cen A using a combined model of infrared emission from startbursts, cirrus clouds and an AGN-type torus. A compact optically thick torus with $`r_{\mathrm{out},\mathrm{torus}}=3.6`$ pc and an opening angle of $`\varphi =30^{}`$ would account for the observed flux. Centaurus A appears as a concentrated IR source in the H band  within a scale of tens of pc. No optical emission has been observed from the central source and, as Falcke et al.  suggested, this source may have a closed torus hiding the BLR which has not been detected even in polarized flux. Israel et al. , Rydbeck et al.  and Turner et al.  concluded that there is a dense nuclear torus with diameter $`<230`$ pc. It is also interesting that Bryant et al.  show a clear picture of a warped dust lane obscuring the nuclear source – the radio jet being perpendicular to a large-scale torus which seems to be inclined with respect to the plane of the dust-lane.

Turning now to other FR-I radio galaxies, 3C218 (Hydra A) also shows evidence for a nuclear obscuration . 3C270 is another FR-I radio galaxy for which a broad $`H_\alpha `$ emission line has been detected . Ferrarese et al.  proposed that its nucleus is surrounded by a dust torus with a diameter of 120 pc and optical depth $`\tau 1`$. An inclination angle of $`20^{}`$ between the axis of the torus and the line of sight would be in accordance with the detection of un-obscured nuclear optical emission .

Possible evidence against tori in blazars comes from the FR-I radio galaxy M87, which seems to have a very low IR flux which cannot be explained with a standard torus model . The torus in this source could be very diluted, and heated less as a result of the extremely low accretion activity of M87  such that its emission is dwarfed by the jet emission.

At kiloparsec scales, starburst activity dominates the infrared emission, but as one goes deeper into the nucleus of the galaxy, de Koff et al.  found that well organized dust structures were present in FR-I galaxies in the 3CR catalog, and this could indicate that the flow of matter towards the nucleus is rather steady in FR-I, allowing the formation of distinct torus features. These structures tended to be sharply defined small-scales disks with radii less than $`2.5`$ kpc. The detection of large-scale dusty features (sometimes shaped as bars or dust-lanes) suggests that there could be an association between a small-scale infrared torus, such as those in Seyfert galaxies, and a large-scale torus sometimes identified with the dust structure of the galaxy. The structure of the nuclear torus could then depend on the large-scale dust distribution – well-organized kpc dust structures could extend inward towards smaller scales (small-scale torus surrounded by large-scale torus). We define this as a symbiosis between the large and small scale dusty features .

The above discussion leads us to the conclusion that dust could be present everywhere, even in FR-I objects where the torus cannot easily be detected. Since FR-I and BL Lacs are thought to be similar (apart from orientation) and there is some evidence of BLR and torus activity in FR-I, we postulate that the linkage between BLR and tori in quasars discussed by Chiaberge et al.  may apply also in blazars. We shall calculate next the optical depth for $`\gamma `$-$`\gamma `$ absorption in the infrared radiation field of tori.

We know that the IR emission from the torus is strongly related to the activity of the central object. The inner radius depends on $`L_{\mathrm{disk}}`$, and is given by the sublimation radius of the dust :

$$R_{\mathrm{in},\mathrm{torus}}T_{1500}^{2.8}L_{\mathrm{disk},46}^{1/2}pc$$
(11)
where $`L_{\mathrm{disk}}=10^{46}L_{\mathrm{disk},46}`$ erg/s, and $`T_{1500}`$ is the dust temperature in units of 1500 K (the dust sublimation temperature is taken to be 1500 K). The inner torus radius is plotted in Fig. 8 versus disk luminosity for three dust temperatures.

We start with a torus centred on the black hole, symmetric about the jet axis, and having a rectangular cross section with full height $`h`$, and inner and outer radii $`R_{\mathrm{in},\mathrm{torus}}`$ and $`R_{\mathrm{out},\mathrm{torus}}`$. We shall discuss how $`\tau _{\gamma \gamma }`$ changes for different scales of tori, e.g., simulating an open torus by having either a rectangular cross section or a rectangular cross section with the inner edges cut away at angle $`\varphi `$ (see Fig. 9).

As discussed earlier, in the context of the symbiosis between jets and accretion disks, an increase in jet power could be at the expense of disk luminosity. A lower disk luminosity at UV frequencies would reduce the heating of the inner surface of the dusty torus, causing the inner radius of the torus, given by Eq. 11, to be small (Fig. 8). This implies that fat dusty tori ($`R_{\mathrm{in},\mathrm{torus}}h`$) could exist in AGN with very low central activity such as blazars, and this would be in accordance with the closed torus model discussed by Falcke et al. . On the other hand, for a “cold torus” with $`T<1000`$ K, the inner radius would be far from the nucleus for a given luminosity (see Fig. 8). Such inefficient heating of the dust would probably describe best the torus in M87 which seems to be a peculiar object in that it has an extremely low mass accretion rate.

Due to its size, the energy density of the infrared photons remains fairly uniform inside the torus. A flat torus (i.e., $`Rh`$) could approximate the models proposed by Chiaberge et al. ; see the discussion from the previous section and comments by Maiolino et al.  about a flattened BLR distributions. Although flat tori have not been detected, and we are uncertain about the stability of such a configuration of dust, we believe that it is an interesting possibility. For simplicity we model in this section a torus with a single black body temperature such that the IR intensity is constant within the solid angle subtended by the torus. We shall extend the work of Protheroe & Biermann  for different torus geometries. Fig. 9 illustrates the geometry for interaction, at points A and B, of $`\gamma `$-rays with IR photons emitted from the surfaces of the torus.

Accretion rates are smaller in blazars than in quasars. Since the thickness of the torus is related to the accretion inflow and disk luminosity (through radiation pressure effects)  one expects tori in FR-I to be thinner, diluted, and/or less efficiently heated – as is probably the case in M87. Therefore, we have varied $`h`$ between $`1`$ pc and $`3`$ pc, $`R_{\mathrm{in},\mathrm{torus}}`$ between $`0.1`$ pc and 1 pc, and $`R_{\mathrm{out},\mathrm{torus}}`$ between 2 pc and 10 pc. The optical depth from $`z=0`$ to infinity, for $`\gamma `$-rays traveling along the jet axis, is shown in Fig. 10 where we have compared the case of a torus with $`R_{\mathrm{in},\mathrm{torus}}=0.1`$ pc (solid curves) with a torus having $`R_{\mathrm{in},\mathrm{torus}}=1`$ pc (dashed curves) for three different torus heights $`h`$. We find that more $`\gamma `$-rays with energies around $`10^{2.53.5}`$ GeV would be absorbed by $`\gamma `$-$`\gamma `$ interactions when the torus comes closer to the central source. We also find that an open torus geometry ($`\varphi >0`$) can significantly modify $`\tau _{\gamma \gamma }`$ as shown in Fig. 11(a) where it is again seen to have the greatest effect near the pair production threshold. If the torus has a large outer radius, then more IR photons from the upper surface of the torus are available for interaction with $`\gamma `$-rays along the jet. For example, for an extended torus with $`R_{\mathrm{out},\mathrm{torus}}=10`$ pc the absorption of photons with energy above 2 TeV is larger (compare dotted line with upper solid curve in Fig. 11a). As noted previously by Protheroe and Biermann , we can see that for photon energies above about several hundreds of GeV the opacity is large enough such that no TeV $`\gamma `$-rays can emerge if the source is near the centre of the torus.

Risaliti et al.  have shown that a torus is stable if the mass of the dust in the torus does not exceed the dynamical mass flowing toward the center. This constrains the outer torus radius, $`R_{\mathrm{out},\mathrm{torus}}`$, to have values less than 10 pc for a column density of dust along the line of sight that intersects the obscuring medium of $`N10^{24}`$ cm<sup>-2</sup>. FR-I objects could have diluted tori with column densities of $`N10^{22}`$ cm<sup>-2</sup>, or lower, but still be optically thick, and this would impose a upper limit of $`100`$ pc for the torus’ outer radius. However, Maiolino et al.  interpret their data as suggesting that the inner parts of tori are much denser than the outer parts, and find a gradient in the torus’ covering factor that would give $`R_{\mathrm{out},\mathrm{torus}}<20`$ pc for Seyfert galaxies. The covering factor plays a crucial role in these models, as well as the gas mass enclosed within the torus . For simplicity, we take a value of 10 pc as the maximum radius of the dust torus which could be relevant to our analysis.

In Fig. 11(b) we show how the optical depth depends on the distance of $`\gamma `$-ray emission region along the jet. So far, we have assumed that the torus radiates as a perfect black body. However, if the torus is patchy or not optically thick at IR wavelengths we may have diluted black body radiation, $`u_{\mathrm{IR}}\eta _{\mathrm{IR}}aT^4`$ erg/cm<sup>3</sup>, where $`a`$ is the radiation constant. All our curves would be multiplied by $`1/\eta _{\mathrm{IR}}`$ in this case. Alternatively, for a poor heating mechanism related to a low luminosity of the accretion disk the dust could be heated to lower temperatures, e.g. $`T100`$ K, and this results in reducing the GeV to TeV opacity considerably, with the pair production threshold $`\gamma `$-ray energy increasing as $`1/T`$ as illustrated in Fig. 11(b) by the dotted curves. We note however, that Blazejowski et al.  require $`\eta _{IR}0.3`$ to fit $`\gamma `$-ray spectra of OVV quasars with an external inverse Compton model, and this still allows for significant TeV $`\gamma `$-ray absorption in this type of blazar. Extremely low values of $`\eta _{IR}`$ may apply in objects such as M87.

## 5 Constraints on the position of the $`\gamma `$-ray emitting region and conclusions

The ADJ model gives a simplified approach to the symbiosis between disks and jets. It could explain the low central activity in some blazars, and for very low mass accretion rates an ADJ could transform into an EDAF or an ADAF. The strong variability observed in blazars from radio to $`\gamma `$-ray frequencies could be a consequence of the feeding mechanism of the jet at the boundary layer. A local perturbation produced in the boundary layer, such as a rapid variation of $`R_{\mathrm{jet}}`$ or $`\dot{M}`$, would propagate through the nozzle and along the jet causing transient phenomena. As the jet extracts more energy from the disk, this could leave less energy to be dissipated in the disk as radiation. The jet/disk symbiosis then mainly modifies the energetics in the central parsec of AGN: increasing the power of the jet could weaken the disk which therefore would not efficiently ionize the BLR clouds or efficiently heat the dust in the torus, and so would affect important target photon fields for interactions of $`\gamma `$-rays (and also accelerated electrons and protons) along the jet.

Fig. 12 summarizes our results in a way that makes it easy to see the relative importance of the three radiation fields to the absorption of high energy $`\gamma `$-rays. We show the optical depth as a function of height of the emission region above the plane of the disk for three $`\gamma `$-ray energies spanning the range currently explored with optical Cherenkov telescopes: (a) 0.5 TeV which is close to the threshold of present telescopes, (b) 2.5 TeV, (c) 12.5 TeV typically observed during flaring activity in nearby blazars. The shaded areas give the range in optical depth for the models discussed earlier, but normalized for an accretion disk luminosity $`L_{\mathrm{disk}}=10^{46}`$ erg s<sup>-1</sup>, a broad line region luminosity $`L_{\mathrm{BLR}}=10^{44}`$ erg s<sup>-1</sup>, and torus IR dilution factor $`\eta _{\mathrm{IR}}=1`$. We see that at all three energies, all three radiation fields can play a very important role. We shall consider each radiation field in turn, first for quasars and then for blazars.

In quasars, with typically $`L_{\mathrm{disk}}10^{46}`$ erg s<sup>-1</sup> the accretion disk radiation will surely cut off $`\gamma `$-rays above 0.5 TeV if the emission region is below $`0.005`$–0.03 pc (depending on $`\dot{m}`$). At 2.5 TeV and 12.5 TeV the corresponding minimum emission region heights to avoid absorption in quasars are $`0.02`$–0.08 pc and $`0.03`$–0.2 pc. The BLR is active in quasars, and is likely to have $`L_{\mathrm{BLR}}10^{44}`$ erg s<sup>-1</sup> such that the emission region would have to be above 0.03–0.3 pc (depending on BLR cloud distribution) at all three energies. For a luminous disk as in quasars, one would expect an active torus to cut off TeV energy $`\gamma `$-rays unless the emission region were well above the torus: $`z_0>0`$–2 pc (0.5 TeV, note $`\tau _{\gamma \gamma }^{\mathrm{IR}}`$ can be small at the energy), $`z_0>2`$–10 pc (2.5 TeV), and $`z_0>4`$–30 pc (12.5 TeV), depending on torus geometry. All known quasars are at a distance such that multi-TeV $`\gamma `$-rays are expected to be absorbed in the infrared background radiation during their propagation to Earth. However, observations at energies below 1 TeV are possible with current atmospheric Cherenkov telescopes, and multi-GeV observations will be made with GLAST  within a few years. Future measurements of the shape of the $`\gamma `$-ray spectrum of quasars at the cut-off might have something to say about which radiation field is responsible for the cut-off, and hence the location of the emission region, although this is likely to be very difficult (compare the optical depth vs. energy in the disk radiation, BLR radiation and torus radiation in Figs. 2, 6 and 11).

The situation is somewhat different in blazars. Here, we can have significantly less powerful accretion disks than in quasars. For example, to avoid attenuation of $`\gamma `$-rays in the disk radiation for a disk luminosity of $`10^{44}`$$`10^{46}`$ erg s<sup>-1</sup> one would need to be above $`0.001`$–0.005 pc (0.5 TeV), and above $`0.005`$–0.03 pc ($`>2.5`$ TeV), depending on $`\dot{m}`$. Note that at all energies the optical depth in the BLR field becomes less than unity if the BLR luminosity is lower than $`10^{42}`$ erg s<sup>-1</sup>. However, for $`L_{\mathrm{BLR}}=10^{42}`$$`10^{44}`$ erg s<sup>-1</sup> the $`\gamma `$-ray emission region would need to be above 0.03–0.3 pc (depending on BLR cloud distribution), as for quasars, to avoid attenuation. In the case of the torus radiation, if blazars have active tori, they could present a serious problem for the escape of TeV $`\gamma `$-rays from emission regions which are not well above the torus (as in the case of quasars). However, because of the low disk luminosity, the torus’ inner radius would probably be small ($`0.1`$ pc for $`L_{\mathrm{disk}}10^{44}`$ erg s<sup>-1</sup>, see Eq. 11), and the torus would be likely to be small, generally, and the optical depth would probably be at or below the lower boundary of that shown. Inefficient heating could also result in dilution of the torus radiation ($`\eta _{\mathrm{IR}}<1`$) or in a lower temperature, and consequently lower optical depths (see e.g. Fig. 11b).

We have demonstrated that the external radiation fields present within the parsec scale of AGN may create problems for the escape of high energy $`\gamma `$-rays. The effect is dramatic for quasars where the densities of external radiation fields are much higher than in blazars. The symbiosis between the accretion disk and jets could result in lower radiation fields around the base of the jet during flaring activity if the jet takes a higher fraction of the matter entering the disk, and this could minimize the dramatic effect of the $`\gamma `$-ray absorption. However, although heated by the disk, the torus IR field would not react rapidly to changes in disk luminosity because of its larger scale, and so the escape of TeV $`\gamma `$-rays would require the emission region to be located at parsec distances from the black hole, far above the torus unless the torus is very inefficiently heated although even in this case one can expect some attenuation. This is in accordance with conclusion of Protheroe and Biermann , and also the recent work of Jorstad et al.  who concluded that the temporal connection between the radio and $`\gamma `$-ray variations suggests that the $`\gamma `$-rays appear to originate far away from the radio core.

In conclusion, the origin of the 100 GeV-1000 GeV $`\gamma `$-ray photons must be above the inner part of BLR region for any quasar. Blazars usually have low central luminosities, and this translates into a having a low BLR luminosity, and possibly a closed torus. Mrk 421 seems to have a BLR with an extremely low luminosity $`L_{\mathrm{BLR}}1.5\times 10^{40}`$ erg/s  and so could be transparent to TeV $`\gamma `$-rays. However, TeV $`\gamma `$-rays would need to be produced far above the torus unless the torus is very inefficiently heated. Farther along the jet external radiation fields become less important. In this case the radiation produced inside the emission region in the jet (e.g. synchrotron emission) becomes the dominant radiation field for photon-photon pair production by $`\gamma `$-rays produced in the jet.

## Acknowledgments

We thank Peter Biermann and Geoff Bicknell for useful discussions. This work was supported by an Australian Research Council grant to RJP.