format: text

text:
# Untitled Document

US-FT-6/02 hep-th/0212139

Singularity Resolution in Gauged Supergravity and Conifold Unification

J. D. Edelstein$`^{\mathrm{}}`$ 1 jedels@math.ist.utl.pt, A. Paredes$`^{}`$ 2 angel@fpaxp1.usc.es and A.V. Ramallo$`^{}`$ 3 alfonso@fpaxp1.usc.es

$`^{\mathrm{}}`$ Departamento de Matemática, Instituto Superior Tecnico

Av. Rovisco Pais, 1049–001, Lisboa, Portugal

$`^{}`$ Instituto de Física de La Plata – Conicet, Universidad Nacional de La Plata

C.C. 67, (1900) La Plata, Argentina

$`^{}`$Departamento de Física de Partículas, Universidad de Santiago de Compostela

E-15782 Santiago de Compostela, Spain

We obtain a unified picture for the conifold singularity resolution. We propose that gauged supergravity, through a novel prescription for the twisting, provides an appropriate framework to smooth out singularities in the context of gravity duals of supersymmetric gauge theories.

11 December 2002

1. Introduction

String theory compactifications in a Calabi–Yau threefold have been the focus of a countless number of papers due to the fact that they provide effective four dimensional vacua with $`𝒩=1`$ supersymmetry. Particular emphasis has been given to the study of singularities on these manifolds –in particular, conical singularities– as long as non-trivial phenomena take place on them such as gauge symmetry enhancement or the appearance of new massless particles. The archetype of these is the well–known conifold. It is a complex three–manifold which is a cone over the homogeneous space $`T^{1,1}=\frac{SU(2)\times SU(2)}{U(1)}`$. The conical singularity can be resolved in two different ways according to whether an $`S^2`$ or an $`S^3`$ is blown up at the singular point. The former is known as the resolved (or Kähler deformed) conifold, while the latter is the (complex) deformed conifold. Both regular manifolds depend on a single parameter (namely, the resolution $`a`$ and the deformation $`\mu `$), are non compact and asymptotically behave as the singular conifold. That is, the three solutions display the same UV behavior for the associated gauge theories. Supersymmetry and matching holonomy conditions in the context of String theory  ensure that there must exist manifolds with $`G_2`$ holonomy metrics whose Gromov–Hausdorff limits are precisely the Ricci flat Kähler metrics on the resolved and deformed conifolds. These were explicitly found in . It was shown afterwards that these $`G_2`$ manifolds arise as solutions of the same system of first–order equations, this providing a nicely unified picture of the resolved and deformed conifolds from the perspective of M–theory .

In a different approach based on lower dimensional gauged supergravity, it was recently shown that the resolved conifold comes out when studying the gravity dual of $`D6`$–branes wrapping an holomorphic $`S^2`$ in a $`K3`$ manifold . The low–energy dynamics is governed, when the size of the cycle is taken to zero, by a five dimensional supersymmetric gauge theory with eight supercharges. If the manifold is large enough and smooth, the dual description is given in terms of a purely gravitational configuration of eleven dimensional supergravity which is the direct product of Minkowski five dimensional spacetime and the resolved conifold. The general solution to this system was later shown  to be given by the generalized resolved conifold . This is a one-parameter (say, $`b`$) generalization of the resolved conifold. There is an analogous extension metric both for the deformed conifold –though it is not regular–, and the singular conifold. We will call the latter regularized conifold, following , because $`b`$ smoothens the curvature singularity and the metric is regular upon imposing a $`\text{ZZ}_2`$ identification of the $`U(1)`$ fiber. It is an ALE space that asymptotically approaches $`T^{1,1}/\text{ZZ}_2`$.

Lower dimensional gauged supergravities provide an explicit arena to impose the twisting conditions required to wrap a $`D`$–brane in a supersymmetric cycle . (See  for a recent review.) Loosely speaking, in the conventional twist, the gauge connection has to be identified with the spin connection. This notion of the twist can be generalized, as shown in , in a way that involves non-trivially the scalar fields that arise in lower dimensional gauged supergravity from the external components of the metric. The solutions obtained by these means usually correspond to the near horizon limit of wrapped D–branes . However, in a recent paper , we have shown that the twist can be further generalized so that it encompasses much more general solutions either corresponding to wrapped D–branes or to special holonomy manifolds with certain RR fluxes turned on. On the one hand, the new twisting condition can be thought of as a non-trivial embedding of the world–volume in the lower dimensional spacetime. More interestingly, as we will show in this paper through an archetypical example, it can also be understood as a singularity resolution mechanism 4 Notice that, in a sense, it is natural to expect that lower dimensional gauged supergravity degrees of freedom cure singularities. For example, even when using the standard twisting conditions, the resolution of the conifold singularity has been shown to be driven by turning on a scalar field in gauged supergravity .: the ordinary twisting imposes the value of the gauge fields at infinity, while the lower dimensional gauged supergravity governs the nontrivial dynamics in the bulk. This mechanism resembles that used in the Maldacena–Núñez solution , where the singularity is solved by turning on a non–Abelian gauge field that asymptotically approaches the Abelian one that twists the gauge theory.

In this letter we present a unified scenario for conifold singularity resolutions. In a sense, we are providing the unified picture of the resolved and deformed conifolds from the perspective of M–theory advocated in . The main difference being that we deal with conifolds in eleven dimensions instead of $`G_2`$ manifolds. A unique system encompass the generalized resolution and deformation of the conifold singularity, each of them emerging as the only two possible solutions of an algebraic constraint. Notice the difference with the $`G_2`$ case studied in , where the algebraic constraints are involved enough so as to admit several well distinct solutions. Here there are only two. Furthermore, we show that it is possible to impose at the same time both solutions of the algebraic constraint, this leading to the regularized conifold metric, which describes a complex line bundle over $`S^2\times S^2`$.

2. D6–branes wrapped on $`𝐒^\mathrm{𝟐}`$ revisited

The Lagrangian describing the dynamics of the sector of Salam and Sezgin’s eight dimensional gauged supergravity  on which we would like to focus (entirely coming from the eleven dimensional metric), reads

$$e^1=\frac{1}{4}R\frac{1}{4}e^{2\varphi }(F_{\mu \nu }^i)^2\frac{1}{4}(P_{\mu ij})^2\frac{1}{2}(_\mu \varphi )^2\frac{1}{32}e^{2\varphi }(e^{8\lambda }4e^{2\lambda }),$$

where $`\varphi `$ is the dilaton, $`\lambda `$ is a scalar in the coset $`SL(3,\mathrm{IR})/SO(3)`$, $`A^i`$ is an $`SU(2)`$ gauge potential, $`e`$ is the determinant of the vierbein $`e^a`$, $`F^i`$ is the Yang–Mills field strength and $`P_{ij}`$ is a symmetric and traceless 1–form defined by

$$P_{ij}+Q_{ij}=\left(\begin{array}{ccccc}d\lambda & & A^3& & A^2e^{3\lambda }\\ \\ A^3& & d\lambda & & A^1e^{3\lambda }\\ \\ A^2e^{3\lambda }& & A^1e^{3\lambda }& & 2d\lambda \end{array}\right),$$

$`Q_{ij}`$ being the antisymmetric counterpart. The supersymmetry transformations for the fermions are given by

$$\delta \psi _\gamma =𝒟_\gamma ϵ+\frac{1}{24}e^\varphi F_{\mu \nu }^i\widehat{\mathrm{\Gamma }}_i(\mathrm{\Gamma }_\gamma ^{\mu \nu }10\delta _\gamma ^\mu \mathrm{\Gamma }^\nu )ϵ\frac{1}{288}e^\varphi ϵ_{ijk}\widehat{\mathrm{\Gamma }}^{ijk}\mathrm{\Gamma }_\gamma Tϵ,$$

$$\delta \chi _i=\frac{1}{2}(P_{\mu ij}+\frac{2}{3}\delta _{ij}_\mu \varphi )\widehat{\mathrm{\Gamma }}^j\mathrm{\Gamma }^\mu ϵ\frac{1}{4}e^\varphi F_{\mu \nu i}\mathrm{\Gamma }^{\mu \nu }ϵ\frac{1}{8}e^\varphi (T_{ij}\frac{1}{2}\delta _{ij}T)ϵ^{jkl}\widehat{\mathrm{\Gamma }}_{kl}ϵ,$$

where $`T_{ij}=\mathrm{diag}(e^{2\lambda },e^{2\lambda },e^{4\lambda })`$, $`T=\delta ^{ij}T_{ij}=2e^{2\lambda }+e^{4\lambda }`$, and the covariant derivative is

$$𝒟ϵ=\left(+\frac{1}{4}\omega ^{ab}\mathrm{\Gamma }_{ab}+\frac{1}{4}Q_{ij}\widehat{\mathrm{\Gamma }}^{ij}\right)ϵ.$$

We shall use the following representation for the Dirac matrices:

$$\mathrm{\Gamma }^\mu =\gamma ^\mu \mathrm{II},\widehat{\mathrm{\Gamma }}^i=\gamma _9\tau ^i,$$

where $`\gamma ^\mu `$ are eight dimensional Dirac matrices, $`\gamma _9=i\gamma ^0\gamma ^1\mathrm{}\gamma ^7`$ ($`\gamma _9^2=\mathrm{\hspace{0.17em}1}`$), $`\tau ^i`$ are Pauli matrices and $`\widehat{\mathrm{\Gamma }}^i`$ are the Dirac matrices along the $`SU(2)`$ group manifold, whereas $`\mathrm{\Gamma }_7\mathrm{\Gamma }_r`$ corresponds to the radial direction.

We shall consider the following ansatz for the eight dimensional metric

$$ds_8^2=e^{2f}dx_{1,4}^2+e^{2h}d\mathrm{\Omega }_2^2+dr^2,$$

where $`d\mathrm{\Omega }_2^2=d\theta ^2+\mathrm{sin}^2\theta d\phi ^2`$ is the metric of the unit $`S^2`$. The ansatz for the gauge field is better presented in terms of the triplet of Maurer–Cartan 1-forms on $`S^2`$

$$\sigma ^1=d\theta ,\sigma ^2=\mathrm{sin}\theta d\phi ,\sigma ^3=\mathrm{cos}\theta d\phi ,$$

that obey the conditions $`d\sigma ^i=\frac{1}{2}ϵ_{ijk}\sigma ^j\sigma ^j`$. It is:

$$A^1=g(r)\sigma ^1,A^2=g(r)\sigma ^2,A^3=\sigma ^3.$$

Notice that the twisting in  corresponds to $`g(r)=0`$. We will check at the end that $`g(r)0`$ asymptotically so, from the point of view of the dual gauge theory, the twisting is not modified. The field strength, $`F^i=dA^i+\frac{1}{2}ϵ_{ijk}A^jA^k`$, reads:

$$F^1=g^{}dr\sigma ^1,F^2=g^{}dr\sigma ^2,F^3=(g^21)\sigma ^1\sigma ^2.$$

When uplifted to eleven dimensions, the unwrapped part of the metric should correspond to flat five dimensional Minkowski spacetime. This condition determines the relation $`f=\varphi /3`$ that we impose from now on. Actually, it is not difficult to write down the form of the eleven dimensional metric for the ansatz we are adopting. Let $`w^i`$ for $`i=1,2,3`$ be a set of SU(2) left invariant one forms of the external three sphere satisfying $`dw^i=\frac{1}{2}ϵ_{ijk}w^jw^k`$. Then, the uplifted eleven dimensional metric is:

$$\begin{array}{cc}\hfill ds_{11}^2=& dx_{1,4}^2+e^{2h\frac{2\varphi }{3}}d\mathrm{\Omega }_2^2+e^{\frac{2\varphi }{3}}dr^2+\mathrm{\hspace{0.17em}4}e^{\frac{4\varphi }{3}+2\lambda }\left(w^1+g\sigma ^1\right)^2+\hfill \\ & +\mathrm{\hspace{0.17em}4}e^{\frac{4\varphi }{3}+2\lambda }\left(w^2+g\sigma ^2\right)^2+\mathrm{\hspace{0.17em}4}e^{\frac{4\varphi }{3}4\lambda }\left(w^3+\sigma ^3\right)^2.\hfill \end{array}$$

In order to seek for supersymmetric solutions to the system, we start by subjecting the spinor to the following angular projection

$$\mathrm{\Gamma }_{\theta \phi }ϵ=\widehat{\mathrm{\Gamma }}_{12}ϵ,$$

which is imposed by the Kähler structure of the ambient $`K3`$ manifold in which the two–cycle lives. The equations $`\delta \chi _1=\delta \chi _2=0`$ give:

$$\left(\lambda ^{}+\frac{2}{3}\varphi ^{}\right)ϵ=ge^h\mathrm{sinh}3\lambda \widehat{\mathrm{\Gamma }}_1\mathrm{\Gamma }_\theta \mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123}ϵe^{\varphi +\lambda h}g^{}\widehat{\mathrm{\Gamma }}_1\mathrm{\Gamma }_\theta ϵ\frac{1}{4}e^{\varphi 4\lambda }\mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123}ϵ,$$

while $`\delta \chi _3=0`$ reads:

$$\begin{array}{cc}\hfill \left(\mathrm{\hspace{0.17em}2}\lambda ^{}\frac{2}{3}\varphi ^{}\right)ϵ=& \left[e^{\varphi 2\lambda 2h}(g^21)\frac{1}{4}e^\varphi \left(e^{4\lambda }\mathrm{\hspace{0.17em}2}e^{2\lambda }\right)\right]\mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123}ϵ\hfill \\ & +2ge^h\mathrm{sinh}3\lambda \widehat{\mathrm{\Gamma }}_1\mathrm{\Gamma }_\theta \mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123}ϵ.\hfill \end{array}$$

One can combine these two equations to eliminate $`\lambda ^{}`$:

$$\varphi ^{}ϵ+e^{\varphi +\lambda h}g^{}\widehat{\mathrm{\Gamma }}_1\mathrm{\Gamma }_\theta ϵ+\left[\frac{1}{2}e^{\varphi 2\lambda 2h}(g^21)+\frac{1}{8}e^\varphi \left(e^{4\lambda }+\mathrm{\hspace{0.17em}2}e^{2\lambda }\right)\right]\mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123}ϵ=0.$$

From this last equation, it is clear that the supersymmetric parameter must satisfy a projection of the sort:

$$\mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123}ϵ=\left(\beta +\stackrel{~}{\beta }\widehat{\mathrm{\Gamma }}_1\mathrm{\Gamma }_\theta \right)ϵ,$$

where $`\beta `$ and $`\stackrel{~}{\beta }`$ are (functions of the radial coordinate) proportional to the first derivatives of $`\varphi ^{}`$ and $`g^{}`$:

$$\varphi ^{}=\left[\frac{1}{2}e^{\varphi 2\lambda 2h}(g^21)+\frac{1}{8}e^\varphi \left(e^{4\lambda }+\mathrm{\hspace{0.17em}2}e^{2\lambda }\right)\right]\beta ,$$

$$e^{\varphi +\lambda h}g^{}=\left[\frac{1}{2}e^{\varphi 2\lambda 2h}(g^21)+\frac{1}{8}e^\varphi \left(e^{4\lambda }+\mathrm{\hspace{0.17em}2}e^{2\lambda }\right)\right]\stackrel{~}{\beta }.$$

This radial projection provides the generalized twist introduced in . It encodes a non-trivial fibering of the two sphere with the external three sphere as will become clear below. Since $`(\mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123})^2ϵ=ϵ`$ and $`\{\mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123},\widehat{\mathrm{\Gamma }}_1\mathrm{\Gamma }_\theta \}=0`$, one must have $`\beta ^2+\stackrel{~}{\beta }^2=1`$. Thus, we can represent $`\beta `$ and $`\stackrel{~}{\beta }`$ as

$$\beta =\mathrm{cos}\alpha ,\stackrel{~}{\beta }=\mathrm{sin}\alpha .$$

Also, it is clear that both independent projections (2.1) and (2.1) leave unbroken eight supercharges as expected. Inserting the radial projection (2.1), as well as (2.1), in (2.1), we get an equation determining $`\lambda ^{}`$:

$$\lambda ^{}=ge^h\mathrm{sinh}3\lambda \stackrel{~}{\beta }\left[\frac{1}{3}e^{\varphi 2\lambda 2h}(g^21)\frac{1}{6}e^\varphi \left(e^{4\lambda }e^{2\lambda }\right)\right]\beta ,$$

together with an algebraic constraint:

$$ge^h\mathrm{sinh}3\lambda \beta +\left[\frac{1}{2}e^{\varphi 2\lambda 2h}(g^21)\frac{1}{8}e^\varphi \left(e^{4\lambda }\mathrm{\hspace{0.17em}2}e^{2\lambda }\right)\right]\stackrel{~}{\beta }=\mathrm{\hspace{0.17em}0}.$$

Let us now consider the equations obtained from the supersymmetric variation of the gravitino. From the components along the unwrapped directions one does not get anything new, while from the angular components we get:

$$\begin{array}{cc}\hfill h^{}ϵ=& ge^h\mathrm{cosh}3\lambda \widehat{\mathrm{\Gamma }}_1\mathrm{\Gamma }_\theta \mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123}ϵ+\frac{2}{3}e^{\varphi +\lambda h}g^{}\widehat{\mathrm{\Gamma }}_1\mathrm{\Gamma }_\theta ϵ\hfill \\ & \frac{1}{6}\left[5e^{\varphi 2\lambda 2h}(g^21)+\frac{1}{4}e^\varphi \left(\mathrm{\hspace{0.17em}2}e^{2\lambda }+e^{4\lambda }\right)\right]\mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123}ϵ.\hfill \end{array}$$

By using the projection (2.1) we obtain an equation for $`h^{}`$:

$$h^{}=ge^h\mathrm{cosh}3\lambda \stackrel{~}{\beta }+\frac{1}{6}\left[5e^{\varphi 2\lambda 2h}(g^21)+\frac{1}{4}e^\varphi \left(\mathrm{\hspace{0.17em}2}e^{2\lambda }+e^{4\lambda }\right)\right]\beta ,$$

together with a second algebraic constraint:

$$ge^h\mathrm{cosh}3\lambda \beta +\left[\frac{1}{2}e^{\varphi 2\lambda 2h}(g^21)\frac{1}{8}e^\varphi \left(\mathrm{\hspace{0.17em}2}e^{2\lambda }+e^{4\lambda }\right)\right]\stackrel{~}{\beta }=\mathrm{\hspace{0.17em}0}.$$

Finally, from the radial component of the gravitino we get the functional dependence of the supersymmetric parameter $`ϵ`$:

$$_rϵ=\frac{5}{6}e^{\varphi +\lambda h}g^{}\widehat{\mathrm{\Gamma }}_1\mathrm{\Gamma }_\theta ϵ\frac{1}{12}\left[e^{\varphi 2\lambda 2h}(g^21)+\frac{1}{4}\left(\mathrm{\hspace{0.17em}2}e^{2\lambda }+e^{4\lambda }\right)\right]\mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123}ϵ.$$

The projection (2.1) gives the generalized twisting conditions first studied in . Its interpretation follows from a similar reasoning as the one used in that reference: using the trigonometric parametrization (2.1), the generalized projection can be written as:

$$\mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123}ϵ=e^{\alpha \widehat{\mathrm{\Gamma }}_1\mathrm{\Gamma }_\theta }ϵ,$$

which can be solved as:

$$ϵ=e^{\frac{1}{2}\alpha \widehat{\mathrm{\Gamma }}_1\mathrm{\Gamma }_\theta }ϵ_0,\mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123}ϵ_0=ϵ_0.$$

We can determine $`ϵ`$ by plugging (2.1) into the equation for the radial component of the gravitino (2.1). Using (2.1), we get two equations. The first one gives the characteristic radial dependence of $`ϵ_0`$ in terms of the eight dimensional dilaton, namely:

$$_rϵ_0=\frac{\varphi ^{}}{6}ϵ_0ϵ_0=e^{\frac{\varphi }{6}}\eta ,$$

with $`\eta `$ being a constant spinor. The other equation determines the radial dependence of the phase $`\alpha `$:

$$\alpha ^{}=2e^{\varphi +\lambda h}g^{}.$$

Thus, the spinor $`ϵ`$ can be written as:

$$ϵ=e^{\frac{\varphi }{6}}e^{\frac{1}{2}\alpha \widehat{\mathrm{\Gamma }}_1\mathrm{\Gamma }_\theta }\eta ,\mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123}\eta =\eta ,\mathrm{\Gamma }_{\theta \phi }\widehat{\mathrm{\Gamma }}_{12}\eta =\eta .$$

The meaning of the phase $`\alpha `$ can be better understood by using the following $`\mathrm{\Gamma }`$–matrices identity $`\mathrm{\Gamma }_{x^0\mathrm{}x^4}\mathrm{\Gamma }_{\theta \phi }\mathrm{\Gamma }_r\widehat{\mathrm{\Gamma }}_{123}=1`$, so that

$$\mathrm{\Gamma }_{x^0\mathrm{}x^4}\left(\mathrm{cos}\alpha \mathrm{\Gamma }_{\theta \phi }\mathrm{sin}\alpha \mathrm{\Gamma }_\theta \widehat{\mathrm{\Gamma }}_2\right)ϵ=ϵ,$$

which shows that the D6–brane is wrapping a two–cycle which is non-trivially embedded in the $`K3`$ manifold as seen from the uplifted perspective that is implied in (2.1). Indeed, the case $`\alpha =0`$ corresponds to the D6–brane wrapping a two–sphere that is fully contained in the eight dimensional spacetime where supergravity lives, studied in .

3. Solution of the algebraic constraints

In the previous section we derived two algebraic constraints (2.1) and (2.1) that the functions of our ansatz must obey. Let us presently solve them. By adding and subtracting the two equations, we get:

$$\mathrm{tan}\alpha \frac{\stackrel{~}{\beta }}{\beta }=2ge^{\varphi +\lambda h}=\frac{ge^{3\lambda h}}{e^{\varphi 2\lambda 2h}(g^21)\frac{1}{4}e^{\varphi 4\lambda }}.$$

Whereas the first part of this equation allows us to write $`\alpha `$ in terms of the remaining functions, the last equality provides an algebraic constraint that restricts our ansatz. It is not hard to arrive at the following simple equation:

$$g\left[g^21+\frac{1}{4}e^{2\varphi 2\lambda +2h}\right]=\mathrm{\hspace{0.17em}0}.$$

There are obviously two solutions. The first one is clearly $`g=0`$, which corresponds to $`\stackrel{~}{\beta }=0`$, $`\beta =1`$, or $`\alpha =0`$. One can check that this is a consistent truncation of the system of equations that actually reduce to the case studied in , whose integral is the generalized resolved conifold (see also ). Indeed, the resulting eleven dimensional metric can be written as $`ds_{11}^2=dx_{1,4}^2+ds_6^2`$, where the six dimensional metric $`ds_6^2`$ is:

$$\begin{array}{cc}\hfill ds_6^2& =\left[\kappa (\rho )\right]^1d\rho ^2+\frac{\rho ^2}{9}\kappa (\rho )(d\psi +\underset{a=1,2}{}\mathrm{cos}\theta _ad\varphi _a)^2\hfill \\ & +\frac{1}{6}(\rho ^2+\mathrm{\hspace{0.17em}6}a^2)(d\theta _1^2+\mathrm{sin}^2\theta _1d\varphi _1^2)+\frac{1}{6}\rho ^2(d\theta _2^2+\mathrm{sin}^2\theta _2d\varphi _2^2),\hfill \end{array}$$

with $`\kappa (\rho )`$ being:

$$\kappa (\rho )=\frac{\rho ^6+\mathrm{\hspace{0.17em}9}a^2\rho ^4b^6}{\rho ^6+\mathrm{\hspace{0.17em}6}a^2\rho ^4}.$$

In eq. (3.1), $`\rho `$ is a new radial variable, $`\theta _1\theta `$, $`\varphi _1\phi `$ and $`(\theta _2,\varphi _2,\psi )`$ parametrize the $`w^i`$’s. The constants $`a`$ and $`b`$ provide the generalized resolution of the conifold singularity . In the context of gauged supergravity, even when this solution corresponds to the conventional twist, $`a`$ and $`b`$ are non-zero when certain scalar fields are excited . The case $`a=0`$, $`b0`$ corresponds to the above mentioned regularized conifold .

The other solution to eq. (3.1) gives a non trivial relation between $`g`$ and the remaining functions of the ansatz, namely:

$$g^2=\mathrm{\hspace{0.17em}1}\frac{1}{4}e^{2\varphi 2\lambda +2h}.$$

It is not difficult to find the values of $`\beta `$ and $`\stackrel{~}{\beta }`$ for this solution of the constraint:

$$\beta =\frac{1}{2}e^{\varphi \lambda +h},\stackrel{~}{\beta }=g.$$

Notice that they satisfy $`\beta ^2+\stackrel{~}{\beta }^2=1`$ as a consequence of the relation (3.1). Moreover, one can verify that (3.1) is consistent with the first-order equations. Indeed, by differentiating eq. (3.1) and using the first-order equations for $`\varphi `$, $`\lambda `$ and $`h`$ (eqs. (2.1), (2.1) and (2.1)), we arrive precisely at the first-order equation for $`g`$ written in (2.1). It can be also checked that eq. (2.1) is identically satisfied for this solution of the constraint. Thus, one can eliminate $`g`$ from the first-order equations arriving at the following system of equations for $`\varphi `$, $`\lambda `$ and $`h`$:

$$\begin{array}{cc}\hfill \varphi ^{}=& \frac{1}{8}e^{2\varphi +\lambda +h},\hfill \\ \hfill \lambda ^{}=& \frac{1}{24}e^{2\varphi +\lambda +h}\frac{1}{2}e^{3\lambda h}+\frac{1}{2}e^{3\lambda h},\hfill \\ \hfill h^{}=& \frac{1}{12}e^{2\varphi +\lambda +h}+\frac{1}{2}e^{3\lambda h}+\frac{1}{2}e^{3\lambda h}.\hfill \end{array}$$

4. The generalized deformed conifold

In order to integrate the system (3.1), let us define the function $`z=\varphi +\lambda h`$ and a new radial coordinate $`\tau `$, $`dr=\mathrm{\hspace{0.17em}2}e^{\varphi 2\lambda }d\tau `$. Then, if the dot denotes the derivative with respect to $`\tau `$, it follows from (3.1) that $`z`$ satisfies the equation:

$$\dot{z}=\frac{1}{2}e^z2e^z.$$

This equation can be immediately integrated:

$$e^z=\frac{1}{2}\frac{\mathrm{cosh}(\tau +\tau _0)}{\mathrm{sinh}(\tau +\tau _0)},$$

where $`\tau _0`$ is an integration constant, which from now on we will absorb in a redefinition of the origin of $`\tau `$. We can obtain $`\varphi `$ by noticing that it satisfies the equation:

$$\dot{\varphi }=\frac{1}{4}e^z.$$

Since we know $`z(\tau )`$ explicitly, we can obtain immediately $`\varphi (\tau )`$, namely:

$$e^\varphi =\widehat{\mu }\left(\mathrm{cosh}\tau \right)^{\frac{1}{2}},$$

where $`\widehat{\mu }`$ is a constant of integration. Finally, $`h`$ satisfies the following differential equation:

$$\dot{h}=\frac{1}{6}e^z+e^z+e^{6\varphi 5z6h}.$$

If we define, $`y=e^{6h}`$ and use the expressions of $`z`$ and $`\varphi `$ as functions of $`\tau `$, we get:

$$\dot{y}=\frac{\mathrm{cosh}^2\tau +2}{\mathrm{cosh}\tau \mathrm{sinh}\tau }y+192\widehat{\mu }^6\frac{(\mathrm{sinh}\tau )^5}{(\mathrm{cosh}\tau )^2},$$

which is also easily integrated by the method of variation of constants. In order to express the corresponding result, let us define the function:

$$K(\tau )\frac{\left(\mathrm{sinh}2\tau \mathrm{\hspace{0.17em}2}\tau +C\right)^{\frac{1}{3}}}{2^{\frac{1}{3}}\mathrm{sinh}\tau },$$

where $`C`$ is a new constant of integration. Then, $`h`$ is given by:

$$e^h=\mathrm{\hspace{0.17em}3}^{\frac{1}{6}}\mathrm{\hspace{0.17em}2}^{\frac{5}{6}}\widehat{\mu }\frac{\mathrm{sinh}\tau }{(\mathrm{cosh}\tau )^{\frac{1}{3}}}\left[K(\tau )\right]^{\frac{1}{2}}.$$

As we know $`z`$, $`h`$ and $`\varphi `$, we can obtain $`\lambda `$. The result is:

$$e^\lambda =\left(\frac{3}{2}\right)^{\frac{1}{6}}(\mathrm{cosh}\tau )^{\frac{1}{6}}\left[K(\tau )\right]^{\frac{1}{2}}.$$

Finally, we can get $`g`$ from the solution of the constraint (eq. (3.1)), namely:

$$g=\frac{1}{\mathrm{cosh}\tau }.$$

It follows immediately from (4.1) that $`g0`$ as $`\tau \mathrm{}`$, as claimed above. Moreover, by using the explicit form of this solution we can find the value of the phase $`\alpha `$:

$$\mathrm{cos}\alpha =\frac{\mathrm{sinh}\tau }{\mathrm{cosh}\tau },\mathrm{sin}\alpha =\frac{1}{\mathrm{cosh}\tau },$$

Notice that $`\alpha \pi /2`$ when $`\tau 0`$, whereas $`\alpha 0`$ for $`\tau \mathrm{}`$. In order to express neatly the form of the corresponding eleven dimensional metric, let us define the following set of one-forms:

$$\begin{array}{cc}\hfill g^1=& \frac{1}{\sqrt{2}}\left[\sigma ^2w^2\right],g^2=\frac{1}{\sqrt{2}}\left[\sigma ^1w^1\right],\hfill \\ \hfill g^3=& \frac{1}{\sqrt{2}}\left[\sigma ^2+w^2\right],g^4=\frac{1}{\sqrt{2}}\left[\sigma ^1+w^1\right],\hfill \\ \hfill g^5=& \frac{1}{\sqrt{2}}\left[\sigma ^3+w^3\right],\hfill \end{array}$$

and a new constant $`\mu `$, related to $`\widehat{\mu }`$ as $`\mu =2^{\frac{11}{4}}3^{\frac{1}{4}}\widehat{\mu }`$. Then, by using the uplifting formula (2.1), the resulting eleven dimensional metric $`ds_{11}^2`$ can again be written as $`ds_{11}^2=dx_{1,4}^2+ds_6^2`$, where now the six dimensional metric is:

$$\begin{array}{cc}\hfill ds_6^2=\frac{1}{2}\mu ^{\frac{4}{3}}K(\tau )[& \frac{1}{3K(\tau )^3}(d\tau ^2+(g^5)^2)+\mathrm{cosh}^2\left(\frac{\tau }{2}\right)((g^3)^2+(g^4)^2)+\hfill \\ & +\mathrm{sinh}^2\left(\frac{\tau }{2}\right)((g^1)^2+(g^2)^2)],\hfill \end{array}$$

which, for $`C=0`$ is nothing but the standard metric of the deformed conifold, with $`\mu `$ being the corresponding deformation parameter.

The metric (4.1) for $`C0`$ was studied in ref. , where it was shown to display a curvature singularity when $`\mu 0`$. On the contrary, for $`\mu =0`$ and $`C0`$ this metric is regular, after a $`\text{ZZ}_2`$ identification of the $`U(1)`$ fiber, and reduces to the one written in (3.1) for $`a=0`$ and $`b0`$ (the regularized conifold), the parameter $`b`$ being related to the constant $`C`$ . It is not difficult to reobtain this result within our formalism. Notice, first of all, that both solutions of the constraint (3.1) are not incompatible, i.e. one can take $`g=0`$ in eq. (3.1) if $`z=\varphi +\lambda h`$ is fixed to the particular constant value $`e^z=1/2`$. Notice that this is consistent with eq. (4.1). Actually, this value of $`z`$ can be obtained by taking $`\tau _0\mathrm{}`$ in the general solution (4.1). Moreover, the differential equation (4.1) for $`\varphi `$ in this $`g=0`$ case reduces to $`\dot{\varphi }=1/2`$, which can be immediately integrated to give $`e^\varphi =Ae^{\tau /2}`$, with $`A`$ being a non-zero constant. Again, this solution can be obtained from the general expression (4.1) by first reintroducing the $`\tau _0`$ parameter (i.e. by changing $`\tau \tau +\tau _0`$) and then by taking $`\tau _0\mathrm{}`$ and $`\widehat{\mu }0`$ in such a way that $`\widehat{\mu }e^{\tau _0/2}=\sqrt{2}A`$. Notice that this corresponds to taking $`\mu =0`$, as claimed. It follows from this discussion that the regularized conifold is a boundary in the moduli space separating the regions that correspond to the generalized deformed and resolved conifolds, as depicted in the figure. Notice that we cannot continously connect the deformed and resolved conifolds through a supersymmetric trajectory of non-singular metrics.

Fig. 1. Representation of the moduli space of generalized resolutions of the conifold singularity. The two regions depicted correspond to the two solutions of our constraint. The generalized deformed conifold metric is singular. A point on each of the three lines represents, from left to right, the resolved, regularized and deformed conifold. They meet at a single point, the singular conifold.

5. Summary and discussion

In this letter, we have shown that lower dimensional gauged supergravities are an appropriate framework to resolve singularities in the study of gravity duals of supersymmetric gauge theories arising in D–branes that wrap a supersymmetric cycle. The key ingredient is provided by the novel twist prescription recently introduced in . The value of the gauge fields at infinity implied by the conventional twisting is preserved, the lower dimensional gauged supergravity governing the non-trivial dynamics in the bulk. The singularity resolution takes place by switching on the appropriate fields of the gauged supergravity which correspond to the generalized twisting.

We have presented a unified scenario for conifold singularity resolutions from the perspective of M–theory: a single system encompassing both the generalized resolution and deformation of the conifold singularity, each of them emerging as the only two possible solutions of an algebraic constraint. It might be possible to relate this constraint to those appearing in the study of $`G_2`$ manifolds carried out in  by a reduction mechanism of the sort discussed in .

It would be interesting to understand the meaning of $`b`$ on the dual five dimensional gauge theory. In the regularized conifold it plays the rôle of a mass scale: If a stack of D3–branes and fractional branes is at the tip of the conifold , $`b0`$ breaks the otherwise conformal invariance associated to the AdS factor for small $`\rho `$. See, for example, .

The mechanism presented in this letter must be useful in studying other singularity resolutions. It would be also interesting to understand the appearance of cascading solutions with chiral symmetry breaking occuring in the IR –through the resolution of naked singularities– , in the framework of gauged supergravity. We hope to report on some of these issues in a near future.

Acknowledgements:

It is a pleasure for us to acknowledge valuable comments from Rafael Hernandez, Carlos Núñez, Leopoldo Pando Zayas and Konstadinos Sfetsos. This work has been supported in part by MCyT and FEDER under grant BFM2002-03881, by Xunta de Galicia, by Fundación Antorchas and by Fundacão para a Ciência e a Tecnologia under grants POCTI/1999/MAT/33943 and SFRH/BPD/7185/2001.

References

relax J. Gomis, “D-Branes, Holonomy and M-Theory”, hep-th/0103115. relax M. Cvetic, G. W. Gibbons, H. Lu and C. N. Pope, “M-theory conifolds”, Phys. Rev. Lett. 88 (2002) 121602, hep-th/0112098. relax M. Cvetic, G. W. Gibbons, H. Lu and C. N. Pope, “A $`G_2`$ unification of the deformed and resolved conifolds”, Phys. Lett. B 534 (2002) 172, hep-th/0112138. relax J.D. Edelstein and C. Núñez, “D6 branes and M-theory geometrical transitions from gauged supergravity”, JHEP 0104, (2001) 028, hep-th/0103167. relax J. D. Edelstein, A. Paredes and A. V. Ramallo, “Wrapped branes with fluxes in 8d gauged supergravity”, hep-th/0207127. relax P. Candelas and X. C. de la Ossa, “Comments On Conifolds”, Nucl. Phys. B 342 (1990) 246. relax G. Papadopoulos and A. A. Tseytlin, “Complex geometry of conifolds and 5-brane wrapped on 2-sphere”, Class. Quant. Grav. 18, 1333 (2001), hep-th/0012034. relax L. A. Pando Zayas and A. A. Tseytlin, “3-branes on resolved conifold”, JHEP0011 (2000) 028, hep-th/0010088. relax L. A. Pando Zayas and A. A. Tseytlin, “3-branes on spaces with $`\mathrm{IR}\times S^2\times S^3`$ topology”, Phys. Rev. D63 (2001) 086006, hep-th/0101043. relax M. Bershadsky, C. Vafa and V. Sadov, “D-Branes and Topological Field Theories”, Nucl. Phys. B 463 (1996) 420, hep-th/9511222. relax J. M. Maldacena and C. Núñez, “Supergravity description of field theories on curved manifolds and a no go theorem”, hep-th/0007018. relax J. D. Edelstein, “Large N dualities from wrapped D-branes”, hep-th/0211204. relax R. Hernandez and K. Sfetsos, “An eight-dimensional approach to G(2) manifolds”, Phys. Lett. B 536 (2002) 294, hep-th/0202135. relax H. J. Boonstra, K. Skenderis and P. K. Townsend, “The domain wall/QFT correspondence”, JHEP9901 (1999) 003 \[hep-th/9807137\]. relax J. D. Edelstein, A. Paredes and A. V. Ramallo, “Let’s twist again: General metrics of G(2) holonomy from gauged supergravity”, hep-th/0211203. relax J. M. Maldacena and C. Núñez, “Towards the large N limit of pure $`𝒩=1`$ super Yang Mills”, Phys. Rev. Lett. 86 (2001) 588, hep-th/0008001. relax A. Salam and E. Sezgin, “d=8 Supergravity”, Nucl. Phys. B 258 (1985) 284. relax R. Hernandez and K. Sfetsos, “Holonomy from wrapped branes”, hep-th/0211130. relax I. R. Klebanov and M. J. Strassler, “Supergravity and a confining gauge theory: Duality cascades and (chi)SB-resolution of naked singularities”, JHEP0008 (2000) 052, hep-th/0007191. relax H. Lü and J. F. Vazquez-Poritz, “$`S^1`$-wrapped D3-branes on conifolds,” Nucl. Phys. B 633, 114 (2002), hep-th/0202175.


text:
# A generalised lattice Boltzmann algorithm for the flow of a nematic liquid crystal with variable order parameter.

## 1 Introduction

The lattice Boltzmann (LB) method has been extensively studied as a mesoscopic method of simulating isotropic fluids (eg ). The strengths of the method lie in modelling flow in complex geometries (eg ) or in multi-component flow (eg ). Recently, a number of LB schemes have been developed to represent the flow of anisotropic fluids such as liquid crystals .

Materials which exhibit liquid crystal phases have anisometric molecules  and nematic liquid crystals are of particular interest because of their application in display devices. There is also increasing interest in modelling liquid crystal colloids in which colloidal particles are embedded in a nematic phase . The colloidal particles interact through the distortions and defects which they generate in the nematic elastic field. The forces between the particles and the dynamics of their motion is accessible to experiment (cf Poulin et al ). As a consequence of these additional colloidal interactions the particles rearrange themselves into new structures and the LB method provides a particularly effective way of dealing with the complex boundary conditions in such problems. The method can recover both the nematostatics and nematodynamics of these phases and hence the dynamics of the colloidal phase can be captured and this has been achieved with an earlier LB method developed by the authors . Additionally, the colloidal particles may be an isotropic fluid and work is currently in progress to extend the two phase LB algorithms which have been developed for isotropic fluids to describe a mixture in which one of the phases is a nematic liquid crystal. This current work is based on the algorithms presented in this paper extended to embody the isotropic-nematic interface following the macroscopic description of Rey . The hydrodynamics associated with pair annihilation of line defects is also a matter of current interest and results have been reported using both lattice Boltzmann  and conventional solvers  of the equations for nematodynamics in the presence of a variable order parameters.

The orientational ordering of the molecules in the nematic phase is characterised by a director field, $`n_\alpha (𝐱,t)`$, a unit vector which essentially defines the ‘average orientation’ of the molecules. However, the nematic ordering is more fully characterised by a traceless and symmetric order tensor, $`Q_{\alpha \beta }`$. In this work, for simplicity, we assume that the director is confined to a two dimensional plane, and hence that $`Q_{\alpha \beta }`$ may be written in the form

$$Q_{\alpha \beta }(𝐱,t)=S(2n_\alpha n_\beta \delta _{\alpha \beta })$$
(1)
The principal eigenvector of $`Q_{\alpha \beta }`$ is the director and the principal eigenvalue, $`S(𝐱,t)`$, is the scalar order parameter.

The scheme proposed by Care et al  to recover the Ericksen, Leslie and Parodi (ELP) equations  used two coupled LB equations. These equations govern the continuum fluid dynamics of an incompressible nematic with an order parameter, $`S`$, which is both position and time independent. In the scheme , one of the LB equations carried the momentum and the second carried a vector density corresponding to the director field. The Denniston et al  scheme is also based on two coupled lattice Boltzmann schemes one of which carries a momentum density and the second of which carries an tensor density from which the macroscopic order tensor can be recovered. This latter scheme recovers the Beris–Edwards equations  for the flow of a nematic liquid crystal with variable order parameter.

In the work presented in this paper, a third scheme is presented which recovers the Qian–Sheng  equations for the flow of a nematic liquid crystal with a variable scalar order parameter. The scheme is based on a single LB equation which governs the evolution of a tensor density and from which both the macroscopic order and momentum evolution equations are recovered. It is assumed that the dynamics generated by the Qian–Sheng equations are essentially equivalent to those generated by the Beris-Edwards equations  although this equivalence has not been demonstrated explicitly. Recent work by Sonnet, Maffettone and Virga  provides the basis upon which the variety of schemes with a variable order parameter may be compared. The Qian–Sheng equations, like the Beris–Edwards equations, reduces to the ELP formalism in the limit that the order parameter becomes independent of time and position. There are a number of important differences between the scheme proposed here and that of Denniston. In the scheme described here, the target equations are those of Qian–Sheng rather than Beris–Edwards, there is a single Boltzmann equation, the equilibrium distribution function is isotropic (as is expected on physical grounds), and the scheme is able to recover the full tensorial coupling of the order tensor to the velocity gradient tensor, from a single lattice distribution function.

The target macroscopic equations of the Qian–Sheng scheme are summarised in Section 2, the new scheme is described in Section 3 and a Chapman Enskog analysis of the scheme is presented in Section 4. A number of analytical results are developed from the Qian–Sheng equations in Section 5, which provide the basis for validation of the LB algorithm in Section 6. The conclusions are presented in Section 7.

## 2 Qian-Sheng formalism

In this section we summarise the target macroscopic equations for the LB method. The two governing equations of the Qian scheme are the momentum evolution equation

$$\rho D_tu_\beta =_\beta (P\delta _{\alpha \beta }+\sigma _{\alpha \beta }^d+\sigma _{\alpha \beta }^f+\sigma _{\alpha \beta }^{})$$
(2)
and the order tensor evolution equation

$$J\ddot{Q}_{\alpha \beta }=h_{\alpha \beta }+h_{\alpha \beta }^{}\lambda \delta _{\alpha \beta }\epsilon _{\alpha \beta \gamma }\lambda _\gamma $$
(3)
It is shown by Qian that in the limit of constant order parameter, the solutions of these equations is identical to those obtained from the ELP equations. Throughout this work we use the repeated index notation for summations over Cartesian indices. In the above equations $`D_t=_t+u_\mu _\mu `$ is the convective derivative, $`P`$ is the pressure, $`\sigma _{\alpha \beta }^d`$ is the distortion stress tensor given by

$$\sigma _{\alpha \beta }^d=\frac{F_{LdG}}{(_{\alpha Q_{\mu \nu }})}_\beta Q_{\mu \nu }$$
(4)
and $`\sigma _{\alpha \beta }^f`$ is the stress tensor associated with an externally applied field. In the current work we only consider an externally applied magnetic field for which the stress tensor is given by Landau and Lifshitz

$$\sigma _{\alpha \beta }^f=\frac{1}{4\pi }(H_\alpha H_\beta \frac{1}{2}H^2\delta _{\alpha \beta })$$
(5)
The viscous stress tensor $`\sigma _{\alpha \beta }^{}`$ is given by

$`\sigma _{\alpha \beta }^{}`$ $`=`$ $`\beta _1Q_{\alpha \beta }Q_{\mu \nu }A_{\mu \nu }+\beta _4A_{\alpha \beta }+\beta _5Q_{\alpha \mu }A_{\mu \beta }+\beta _6Q_{\beta \mu }A_{\mu \alpha }+`$ (6)
$`{\displaystyle \frac{1}{2}}\mu _2N_{\alpha \beta }\mu _1Q_{\alpha \mu }N_{\mu \beta }+\mu _1Q_{\beta \mu }N_{\mu \alpha }`$
where $`A_{\alpha \beta }=\frac{1}{2}(_\alpha u_\beta +_\beta u_\alpha )`$ is the symmetric velocity gradient tensor. The elastic molecular field is given by

$$h_{\alpha \beta }=\frac{F_{LdG}}{Q_{\alpha \beta }}+_\mu \frac{F_{LdG}}{(_\mu Q_{\alpha \beta })}$$
(7)
where, in the current work, the Landau-deGennes free energy is assumed to be of the form

$`F_{LdG}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\alpha _FQ_{\mu \nu }^2+L_1Q_{\mu \nu ,\tau }^2+L_2Q_{\mu \nu ,\nu }Q_{\mu \tau ,\tau })`$ (8)
$`\beta _FQ_{\mu \nu }Q_{\nu \tau }Q_{\tau \mu }+\gamma _F(Q_{\mu \nu }^2)^2`$
where $`Q_{\alpha \beta ,\gamma }_\gamma (Q_{\alpha \beta })`$. Using this form of the free energy and the definition of the molecular field, equation (7), we find

$$h_{\alpha \beta }=L_1_\mu ^2Q_{\alpha \beta }+L_2_\beta _\mu Q_{\alpha \mu }\alpha _FQ_{\alpha \beta }+3\beta _FQ_{\alpha \mu }Q_{\beta \mu }4\gamma _FQ_{\alpha \beta }Q_{\mu \nu }^2$$
(9)
The quantities $`\lambda `$ and $`\lambda _\alpha `$ in equation (3) are Lagrange multipliers which impose the constraints on the elastic molecular field which arise because the order tensor, $`Q_{\alpha \beta }`$, is symmetric and traceless. For a three dimensional system they have the values

$$\begin{array}{ccccc}\lambda =\frac{1}{3}(h_{\mu \mu }\frac{1}{2}\mu _2A_{\mu \mu })& & \lambda _\alpha & =& \frac{1}{2}\epsilon _{\alpha \mu \nu }h_{\mu \nu }\end{array}$$
(10)
The term $`A_{\mu \mu }`$ in the first of equations (10) does not appear in the Qian-Sheng equations, but is necessary in order to correct for the slight compressibility of the LB fluid. In the presence of an external magnetic field the free energy is assumed to be augmented by a term of the form

$$F_H=\frac{1}{2}(\chi _{}\chi _{})Q_{\alpha \beta }H_\alpha H_\beta $$
(11)
which gives an additional term in the molecular field of the form

$$h_{\alpha \beta }^H=\frac{1}{2}\chi _aH_\alpha H_\beta $$
(12)
where $`\chi _a=\chi _{}\chi _{}`$ is the anisotropy in the susceptibility. The viscous molecular field, $`h_{\alpha \beta }^{}`$, is given by

$$h_{\alpha \beta }^{}=\frac{1}{2}\mu _2A_{\alpha \beta }+\mu _1N_{\alpha \beta }$$
(13)
where $`N_{\alpha \beta }`$ is the co-rotational derivative defined by

$$N_{\alpha \beta }=_tQ_{\alpha \beta }+u_\mu _\mu Q_{\alpha \beta }\epsilon _{\alpha \mu \nu }\omega _\mu Q_{\nu \beta }\epsilon _{\beta \mu \nu }\omega _\mu Q_{\nu \alpha }$$
(14)
where the vorticity, $`\underset{¯}{\omega }=\frac{1}{2}(\times \underset{¯}{u})`$.

In order to develop a lattice Boltzmann scheme we re-arrange equation (6) by substituting from equation (14) for the co-rotational derivative, $`N_{\alpha \beta }`$, and obtain the form

$`\sigma _{\alpha \beta }^{}`$ $`=`$ $`\beta _1Q_{\alpha \beta }Q_{\mu \nu }A_{\mu \nu }+(\beta _4{\displaystyle \frac{(\mu _2)^2}{2\mu _1}})A_{\alpha \beta }+(\beta _5+{\displaystyle \frac{1}{2}}\mu _2)Q_{\alpha \mu }A_{\mu \beta }+`$ (15)
$`(\beta _6{\displaystyle \frac{1}{2}}\mu _2)Q_{\beta \mu }A_{\mu \alpha }+{\displaystyle \frac{\mu _2}{2\mu _1}}(h_{\alpha \beta }\delta _{\alpha \beta }\lambda \epsilon _{\alpha \beta \mu }\lambda _\mu )+`$
$`Q_{\beta \mu }h_{\mu \alpha }Q_{\alpha \mu }h_{\mu \beta }\epsilon _{\mu \beta \nu }Q_{\alpha \mu }\lambda _\nu \epsilon _{\mu \alpha \nu }Q_{\beta \mu }\lambda _\nu `$
Equation (3) for the order tensor evolution can be cast in the form

$`D_tQ_{\alpha \beta }`$ $`=`$ $`{\displaystyle \frac{\mu _2}{2\mu _1}}A_{\alpha \beta }+{\displaystyle \frac{1}{\mu _1}}h_{\alpha \beta }+ϵ_{\alpha \mu \nu }\omega _\mu Q_{\nu \beta }+ϵ_{\beta \mu \nu }\omega _\mu Q_{\nu \alpha }`$ (16)
$`{\displaystyle \frac{1}{\mu _1}}(\delta _{\alpha \beta }\lambda +\epsilon _{\alpha \beta \mu }\lambda _\mu )`$
where we have set the inertial density of the fluid, $`J`$, to zero. Equations (15) and (16) are the expressions which are recovered by the LB scheme proposed below. After each time step the macroscopic density, momentum and the order tensor are recovered from the LB tensor density. The gradients of the order tensor are then used to construct the molecular field through equation (9) and this modifies the dynamics through appropriate forcing terms in the LB equation.

## 3 The lattice Boltzmann algorithm

In order to recover both the momentum evolution equation (2) and the order tensor evolution equation (3) within a lattice Boltzmann scheme, the scalar density of a standard lattice Boltzmann scheme, $`f_i(\underset{¯}{r},t)`$, is replaced by a tensor density $`g_{i\alpha \beta }(\underset{¯}{r},t)`$ where $`i`$ is the normal velocity index and $`\alpha `$ and $`\beta `$ label either a two, or three, dimensional Cartesian basis. One can think of the tensor density as carrying information about the ordering of that population of the fluid ‘element’ on the velocity link $`i`$, associated with a particular position and time. Hence the densities of a standard LB scheme have been generalised to carry information about the order associated with the fluid in addition to the density and momentum.

The density $`g_{i\alpha \beta }(\underset{¯}{r},t)`$ is assumed to have the following moments

$`\rho `$ $`=`$ $`{\displaystyle \underset{i}{}}g_{i\mu \mu }`$ (17)
$`\rho u_\alpha `$ $`=`$ $`{\displaystyle \underset{i}{}}c_{i\alpha }g_{i\mu \mu }`$ (18)
$`\rho S_{\alpha \beta }`$ $`=`$ $`{\displaystyle \underset{i}{}}g_{i\alpha \beta }`$ (19)
where $`S_{\alpha \beta }`$ is an order matrix with unit trace which is related to $`Q_{\alpha \beta }`$ by

$$S_{\alpha \beta }=\frac{1}{2}(Q_{\alpha \beta }+\delta _{\alpha \beta })$$
(20)
where $`d`$ is the dimensionality of the space which the Cartesian indices $`\alpha `$ and $`\beta `$ span. The LB algorithm governing the evolution of $`g_{i\alpha \beta }`$ is taken to be of the form

$$g_{i\alpha \beta }(\underset{¯}{r}+\delta \underset{¯}{c}_i,t+\delta )=g_{i\alpha \beta }(\underset{¯}{r},t)+\underset{j}{}g_{j\mu \nu }^{(neq)}(\underset{¯}{r},t)M_{i\alpha \beta j\mu \nu }+\varphi _{i\alpha \beta }+\chi _{i\alpha \beta }$$
(21)
where $`\varphi _{i\alpha \beta }`$ and $`\chi _{i\alpha \beta }`$ are the forcing terms for the momentum and order respectively and there is an implicit summation over repeated Cartesian indices. The non-equilibrium distribution function is given by

$$g_{i\alpha \beta }^{(neq)}=g_{i\alpha \beta }g_{i\alpha \beta }^{(0)}$$
(22)
where the equilibrium distribution function, $`g_{i\alpha \beta }^{(0)}`$, is taken to be of the form

$$g_{i\alpha \beta }^{(0)}=S_{\alpha \beta }f_i^{(0)}.$$
(23)
The distribution function, $`f_i^{(0)}`$, is assumed to be second order in the velocity $`u_\alpha `$ and is determined in the usual way, by the requirements

$$\rho =\underset{i}{}f_i^{(0)}\rho u_\alpha =\underset{i}{}c_{i\alpha }f_i^{(0)}$$
(24)
and Galilean invariance in the form

$$\underset{i}{}c_{i\alpha }c_{i\beta }g_{i\alpha \beta }^{(0)}=c_s^2\rho \delta _{\alpha \beta }+\rho u_\alpha u_\beta $$
(25)
Hence, we have

$$\underset{i}{}g_{i\alpha \beta }^{(0)}=\rho S_{\alpha \beta }\underset{i}{}c_{i\alpha }g_{i\alpha \beta }^{(0)}=\rho S_{\alpha \beta }u_\alpha $$
(26)
In Section 4, it is shown that in order to recover the required tensor coupling of the order tensor to the velocity gradient tensor it is necessary to have sixth order isotropy of the velocity tensors. In this work, the required isotropy is achieved by using a $`D2Q13`$ three speed hexagonal lattice with velocity vectors, $`𝐜_p`$, given by

$`𝐜_\mathrm{𝟎}`$ $`=`$ $`\{0,0\}`$
$`𝐜_\mathrm{𝟏}`$ $`=`$ $`c\{\pm 1,0\},c\{\pm 1/2,\pm \sqrt{3/2}\}`$
$`𝐜_\mathrm{𝟐}`$ $`=`$ $`c\{0,\pm \sqrt{3}\},c\{\pm 3/2,\pm \sqrt{3/2}\}`$ (27)
the subscripts 0, 1 and 2 being associated with particles with velocity 0, $`c`$ and $`\sqrt{3c}`$ respectively. Imposing the conditions (24) and (25) leads to the result

$$g_{i\alpha \beta }^{(0)}=\rho S_{\alpha \beta }t_i\left[1+\frac{1}{c_s^2}u_\alpha c_{i\alpha }+\frac{1}{2c_s^2}u_\alpha u_\beta \left(\frac{c_{i\alpha }c_{i\beta }}{c_s^2}\delta _{\alpha \beta }\right)\right]$$
(28)
where, in $`D2Q13`$, $`t_0=11/25,t_1=9/100`$ and $`t_2=1/300`$ and the velocity of sound $`c_s=(3/10)^{1/2}`$. The collision operator is taken to be of the form

$$M_{i\alpha \beta j\mu \nu }=\delta _{\alpha \mu }\delta _{\beta \nu }\left(\frac{\delta _{i,j}}{\tau _j}+\frac{\delta _{i0}}{\tau _j}+\frac{t_ic_{i\epsilon }c_{j\epsilon }}{c_s^2}(\frac{1}{\tau _j}2)\right)$$
(29)
where the direction dependent relaxation parameter is defined to be

$$\tau _i=\tau _0(1+\eta _0c_{i\mu }c_{i\nu }\delta _{\mu \nu }+\eta _2c_{i\mu }c_{i\nu }Q_{\mu \nu })$$
(30)
This collision operator may be seen more transparently after performing the contraction with the non-equilibrium distribution function as in equation (21)

$$\underset{j\mu \nu }{}M_{i\alpha \beta j\mu \nu }g_{j\mu \nu }^{(neq)}=\frac{g_{i\alpha \beta }^{(neq)}}{\tau _i}+\delta m_{i\alpha \beta }+\delta p_{i\alpha \beta }$$
(31)
where

$`\delta m_{i\alpha \beta }`$ $`=`$ $`\delta _{i0}{\displaystyle \underset{j}{}}{\displaystyle \frac{g_{j\alpha \beta }^{(neq)}}{\tau _j}}`$ (32)
is a correction to conserve mass. The term

$`\delta p_{i\alpha \beta }`$ $`=`$ $`{\displaystyle \frac{t_ic_{i\epsilon }}{c_s^2}}{\displaystyle \underset{j}{}}c_{j\epsilon }g_{j\alpha \beta }^{(neq)}\left({\displaystyle \frac{1}{\tau _j}}2\right)`$ (34)
includes a term in $`1/\tau _j`$ which is a correction to conserve momentum and a term associated with a factor $`2`$ which is explained at the end of Section 4.1. The collision operator is essentially an LBGK collision operator (eg ) with a relaxation time which is direction dependent. When written in the matrix form, the collision operator is seen to lie between a conventional LBGK operator and a linearised lattice Boltzmann scheme (eg ) although the usual circulant properties of the matrix associated with an isotropic fluid are destroyed by the direction dependent scattering in the method presented here.

The momentum forcing term is given by

$$\varphi _{i\alpha \beta }=t_iS_{\alpha \beta }c_{i\mu }_\nu F_{\nu \mu }$$
(35)
where

$`F_{\alpha \beta }={\displaystyle \frac{1}{c_s^2}}({\displaystyle \frac{\mu _2(h_{\alpha \beta }\epsilon _{\alpha \beta \mu }\lambda _\mu )}{2\mu _1}}`$ (36)
$`+Q_{\alpha \mu }(h_{\mu \beta }\epsilon _{\mu \beta \nu }\lambda _\nu )Q_{\beta \mu }(h_{\mu \alpha }\epsilon _{\mu \alpha \nu }\lambda _\nu )\sigma _{\alpha \beta }^d\sigma _{\alpha \beta }^f)`$
and the angular forcing term is given by

$`\chi _{i\alpha \beta }=t_i\rho \left[(\epsilon _{\alpha \nu \mu }Q_{\beta \mu }+\epsilon _{\beta \nu \mu }Q_{\alpha \mu })\omega _\nu \right]`$ (37)
$`{\displaystyle \frac{t_i\rho }{2\mu _1}}\left[\mu _2A_{\alpha \beta }2h_{\alpha \beta }+2\delta _{\alpha \beta }\lambda +2\epsilon _{\alpha \beta \mu }\lambda _\mu \right]`$
## 4 Analysis of the algorithm

### 4.1 Chapman Enskog expansion

In this section the key results of the Chapman Enskog analysis are given which demonstrate how the scheme described in Section 2 can be recovered from the LB algorithm described in Section 3. We follow a standard Chapman Enskog analysis () and expand the density and the time derivatives in the form

$`g_{i\alpha \beta }`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\delta ^ng_{i\alpha \beta }^{(n)}`$ (38)
$`_t`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\delta ^n_{t,n}`$ (39)
We assume for the purposes of the following analysis that the forcing terms, $`\varphi _{i\alpha \beta }`$ and $`\chi _{i\alpha \beta }`$, can both be introduced at order $`O(\delta ^2)`$ since both include gradient terms in both the velocity and director field (cf). The choice is supported by the agreement of the Chapman Enskog analysis with the measured simulation data; however the assumption relies upon assumptions implicit in LB hydrodynamics in general. Accordingly, it requires more careful analysis which will be undertaken in later work.

The Chapman Enskog expansion gives to $`O(\delta )`$

$$_{t,0}g_{i\alpha \beta }^{(0)}+c_{i\mu }_\mu g_{i\alpha \beta }^{(0)}=\underset{j}{}g_{j\mu \nu }^{(1)}M_{i\alpha \beta j\mu \nu }$$
(40)
and to $`O(\delta ^2)`$

$`_{t,1}g_{i\alpha \beta }^{(0)}+(_{t,0}+c_{i\mu }_\mu )\left\{g_{i\alpha \beta }^{(1)}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{j}{}}g_{j\mu \nu }^{(1)}M_{i\alpha \beta j\mu \nu }\right\}=`$ (41)
$`{\displaystyle \underset{j}{}}g_{j\mu \nu }^{(2)}M_{i\alpha \beta j\mu \nu }+\varphi _{i\alpha \beta }+\chi _{i\alpha \beta }`$
We now take moments of equations (40) and (41) in order to recover the macroscopic equations to which the algorithm is equivalent. If we sum the two equations over the index $`i`$, take the trace of the distribution density and sum the resulting equations, we recover

$$(_{t,0}+_{t,1})\rho +_\mu (\rho u_\mu )=0$$
(42)
which is the continuity equation to second order. If we multiply equations (40) and (38) by $`c_{i\nu }`$, sum over the index $`i`$, take the trace of the distribution density and sum the resulting equations we recover

$$(_{t,0}+_{t,1})(\rho u_\alpha )+_\mu (\mathrm{\Pi }_{\alpha \mu }^{(0)}+\mathrm{\Pi }_{\alpha \mu }^{(1)}+\frac{1}{2}\mathrm{\Omega }_{\alpha \mu }^{(1)})=c_s^2_\mu F_{\mu \alpha }$$
(43)
where

$$\mathrm{\Pi }_{\alpha \beta }^{(n)}=\underset{i}{}c_{i\alpha }c_{i\beta }g_{i\mu \mu }^{(n)}$$
(44)
and

$$\mathrm{\Omega }_{\alpha \beta }^{(n)}=\underset{ij}{}g_{j\mu \nu }^{(n)}c_{i\alpha }c_{i\beta }M_{i\kappa \kappa j\mu \nu }$$
(45)
Further progress can be made by noting that, to $`O(u)`$, the lowest order equation (40) can be written in the form

$$t_i\left(\frac{c_{i\mu }c_{i\nu }}{c_s^2}\delta _{\mu \nu }\right)_\mu (\rho S_{\alpha \beta }u_\nu )=\underset{j}{}g_{j\mu \nu }^{(1)}M_{i\alpha \beta j\mu \nu }$$
(46)
Substituting into equation (45), this yields the result

$$\mathrm{\Omega }_{\alpha \beta }^{(1)}=\underset{i}{}t_ic_{i\alpha }c_{i\beta }\left(\frac{c_{i\mu }c_{i\nu }}{c_s^2}\delta _{\mu \nu }\right)_\mu (\rho u_\nu )$$
(47)
where the unit trace property of $`S_{\alpha \beta }`$ has been used. In order to obtain an expression for $`\mathrm{\Pi }_{\alpha \beta }^{(1)}`$ we use equations (31) and (46) to give, to first order in the velocity,

$$g_{i\alpha \beta }^{(1)}=\tau _i\left[t_i\left(\frac{c_{i\mu }c_{i\nu }}{c_s^2}\delta _{\mu \nu }\right)_\mu (\rho S_{\alpha \beta }u_\nu )\delta m_{i\alpha \beta }\delta p_{i\alpha \beta }\right]$$
(48)
where we have made the approximation, correct to $`O(\delta ^2)`$, that $`g_{i\alpha \beta }^{(neq)}g_{i\alpha \beta }^{(1)}`$. The term $`\delta m_{i\alpha }`$ does not contribute to $`\mathrm{\Pi }_{\alpha \beta }^{(1)}`$, as can be seen from equation (44), since it is included in the rest mass. It can further be shown that the contribution from the term $`\delta p_{i\alpha \beta }`$ is essentially zero. We therefore find

$`\mathrm{\Pi }_{\alpha \beta }^{(1)}=\tau _0{\displaystyle \underset{i}{}}`$ (49)
$`t_ic_{i\alpha }c_{i\beta }(1+\eta _0c_{i\mu }c_{i\nu }\delta _{\mu \nu }+\eta _2c_{i\mu }c_{i\nu }Q_{\mu \nu })\left[\left({\displaystyle \frac{c_{i\tau }c_{i\epsilon }}{c_s^2}}\delta _{\tau \epsilon }\right)_\tau (\rho u_\epsilon )\right]`$
where we have substituted for the anisotropic relaxation parameter, $`\tau _i`$. Assuming that the velocity tensors

$$E_{\alpha _1\mathrm{}\alpha _n}^{(n)}=\underset{i}{}t_ic_{i\alpha _1}\mathrm{}\mathrm{}c_{i\alpha _n}$$
(50)
are isotropic up to sixth order, we recover the required form of the tensor coupling between the order tensor and the velocity gradient tensor.

Using these results it can be shown that, apart from the term associated with $`\beta _1`$, the required momentum evolution equation (2) is recovered when equation (36) is used as the momentum forcing term. The term in $`\beta _1`$ could be recovered by including a fourth order velocity product in the anisotropic relaxation parameter but this would require the velocity tensors to be isotropic to eighth order. However, for simplicity, the term in $`\beta _1`$ is omitted in the results presented in this paper.

The order tensor evolution equation is recovered from the two equations (40) and (41) by summing over the index $`i`$ and adding the two resulting equations to give

$$_{t,0}(\rho Q_{\alpha \beta })+_{t,1}(\rho Q_{\alpha \beta })+_\gamma (\rho Q_{\alpha \beta }u_\gamma )=\underset{i}{}\chi _{i\alpha \beta }$$
(51)
Equation (51) combined with the forcing term (37) gives the required evolution equation (3). However, in order to obtain the result (51) it is necessary to suppress a contribution of the form

$$_\gamma (\underset{i}{}c_{i\gamma }g_{i\alpha \beta }^{(1)})$$
(52)
at second order in the Chapman Enskog analysis and this is achieved using the term associated with the factor $`2`$ in equation (34).

The term (52) arises because conservation of mass and momentum are alone insufficient to constrain all the first order elements of the distribution function. This arises because in order to recover the required macroscopic equations, the form of the equilibrium distribution function must satisfy equations (26) . However, the limitations on the definitions of the moments of the distribution function lead to the result

$$\underset{i}{}c_{i\gamma }g_{i\alpha \beta }\rho S_{\alpha \beta }u_\gamma $$
(53)
and as a consequence the argument of the derivative in (52) is non-zero.

### 4.2 Choice of LB Parameters

The following correspondence is found between the parameters of the LB algorithm and the parameters of the target scheme

$`\beta _4`$ $`=`$ $`{\displaystyle \frac{\mu _2^2}{4\mu _1}}\rho c_s^2\left(12\tau _0\zeta \right)`$
$`\beta _5`$ $`=`$ $`{\displaystyle \frac{\mu _2}{2}}+\eta _2\rho \tau _0c_s^2`$
$`\beta _6`$ $`=`$ $`{\displaystyle \frac{\mu _2}{2}}+\eta _2\rho \tau _0c_s^2`$ (54)
where the parameter $`\zeta =1+\eta _0(1+d/4)`$ with $`d`$ being the dimensionality of the LB scheme which is taken to be 2 for the results presented in Section (6). The viscosity set, $`\{\beta _4,\beta _5,\mu _1,\mu _2\}`$, of the Qian scheme is recovered from the parameter set, $`\{\tau _0,\eta _0,\eta _2,\mu _1,\mu _2\}`$, of the LB algorithm. We note from Qian() that $`\beta _6\beta _5=\mu _2`$ and this is seen to be consistent with the last two of equations (4.2). There is one free parameter in the LB scheme since we require to recover only four Qian parameters from five LB parameters. The free parameter is taken to be $`\eta _0`$ which is adjusted to place the LB scheme in a stable region of its parameter space, as is explained below.

The equations (4.2) can be inverted to give

$`\eta _2`$ $`=`$ $`{\displaystyle \frac{4\mu _1\zeta (2\beta _5+\mu _2)}{\mu _2^24\mu _1(\beta _4+\rho c_S^2)}}`$
$`\tau _0`$ $`=`$ $`{\displaystyle \frac{4\beta _4+4\rho c_s^2\frac{\mu _2^2}{\mu _1}}{8\rho c_s^2\zeta }}`$ (55)
In the absence of any anisotropic terms, the last expression in equation (4.2) becomes $`\beta _4/\rho =c_s^2(12\tau _0)`$, which is equivalent to the standard result for the kinematic viscosity in an LBGK scheme. The relationship between the Qian parameters $`\{\beta _4,\beta _5,\mu _1,\mu _2\}`$ and the standard Leslie Coefficients is given in Section 5.

We now establish an approximate criterion for the stability of the LB algorithm by considering a system which is evolving towards a uniform density and velocity distribution at all points. The algorithm involves repeated application of the collision operator, $`M_{i\alpha \beta j\mu \nu }\underset{¯}{\underset{¯}{M}}`$, and ignoring forcing terms, the non-equilibrium part of the distribution function is transformed at each time step according to the symbolic equation

$$\underset{¯}{g}^{(neq)}(\underset{¯}{\underset{¯}{1}}+\underset{¯}{\underset{¯}{M}})\underset{¯}{g}^{(neq)}$$
(56)
In order for the non-equilibrium part to decay to zero under successive applications of equation (56), it is necessary that

$$1<(\underset{¯}{\underset{¯}{1}}+\underset{¯}{\underset{¯}{M}})_{i\alpha \beta j\mu \nu }<1$$
(57)
for all choices of indices. This equation may be recast as the stability criteria

$$2<M_{(Diag)}<0|M_{(OffDiag)}|<1$$
(58)
However, it should be remembered that the matrix $`\underset{¯}{\underset{¯}{M}}`$ is a function of position, as a consequence of its dependence on the order tensor through the anisotropic relaxation parameter, $`\tau _i`$. Hence, in a given simulation, the conditions (58) cannot be guaranteed under all flow conditions. The criteria have been achieved in the implementation of the algorithm described in Section 6 by explicit evaluation of all the non-zero elements of the collision matrix prior to running the algorithm and after making appropriate assumptions about the ordering in the system. The parameter $`\eta _0`$ is then chosen to ensure that the conditions (58) are satisfied for all possible values of the director and scalar order parameter.

## 5 Analytical results from the Qian formalism

In this section we quote without detailed proof some analytical results which can be derived from the Qian formalism and which will be used in Section 6 to validate the algorithm. Further discussion of the results and their significance is being prepared for a future publication. The results quoted in Section 6 have been derived for a system in which the director is confined to lie in a two dimensional plane and the order tensor for such a system may be written in the form

$$Q_{\alpha \beta }=S(2n_\alpha n_\beta \delta _{\alpha \beta })$$
(59)
where the time and position dependent quantities $`S`$ and $`n_\alpha `$ are the order parameter and director. If the form (59) is substituted into the free energy expression (8) and all the time and spatial derivatives are removed, we find the equilibrium Landau-deGennes free energy in the presence of a uniform magnetic field is given by

$$F_{LdG}^{(equ)}=\alpha _FS_H^2+4\gamma S_H^4\frac{1}{2}S_H\chi _aH_\alpha ^2$$
(60)
where it is assumed that the director is parallel to the magnetic field. The equilibrium order parameter, which minimises equation (60), has the form

$$S_H=S_++S_{}$$
(61)
with

$$S_\pm =\frac{1}{12}\left(\frac{27\chi _aH^2}{\gamma }\pm 3\sqrt{3}\sqrt{\frac{8\alpha _F^3+27\gamma \chi _a^2H^4}{\gamma ^3}}\right)^{1/3}$$
(62)
where the free energy parameter, $`\alpha `$, will be negative in the nematic phase and we choose the real solutions of $`S_0`$. In the limit the $`H0`$, we have

$$S_0=\frac{1}{2}\sqrt{\frac{\alpha _F}{2\gamma }}$$
(63)
which is the equilibrium order parameter which effectively minimises the Landau-deGennes free energy in the absence of director gradients. Using equation (59) with $`S`$ set to $`S_0`$ we may follow the arguments given in the Appendix (B) of Qian to determine the relationship between the Leslie coefficients of a material and the viscosity coefficients used in the Qian scheme. It is found for a system with the order tensor given by (59)

$$\begin{array}{ccccccc}\hfill \beta _1& =& \alpha _1/(4S_0^2)\hfill & & \hfill \beta _4& =& (1/2)\left(2\alpha _4+\alpha _5+\alpha _6\right)\hfill \\ \hfill \beta _5& =& \alpha _5/(2S_0)\hfill & & \hfill \beta _6& =& \alpha _6/(2S_0)\hfill \\ \hfill \mu _1& =& (\alpha _3\alpha _2)/(8S_0^2)\hfill & & \hfill \mu _2& =& (\alpha _2+\alpha _3)/(2S_0)\hfill \end{array}$$
(64)
where $`S_0`$ is the order parameter defined by equation (63). In simple shear flow between two parallel plates and in the absence of any external field, equation (3) for the evolution of the order tensor can be solved to find the steady state value of the angle of orientation and the order parameter. In the centre of the flow, remote from the walls, the gradients in the order tensor may be ignored and the director is found to lie at an angle $`\theta `$ with respect to the flow velocity where

$$\mathrm{cos}(2\theta )=\frac{4S\mu _1}{\mu _2}$$
(65)
and this is directly equivalent to the result found for the ELP equations, $`\mathrm{cos}(2\theta )=\gamma _1/\gamma _2`$ . The order parameter is also modified in shear flow and is given by the solution of the cubic equation

$$ax^3+bx^2+cx+d=0$$
(66)
where $`x=S^2`$ and

$$\begin{array}{ccccccc}a\hfill & =\hfill & 1024\gamma ^2\hfill & & b\hfill & =\hfill & 256\alpha \gamma \hfill \\ c\hfill & =\hfill & 16(\alpha ^2+(\mu _1_yu_x)^2)\hfill & & d\hfill & =\hfill & (\mu _2_yu_x)^2\hfill \end{array}$$
(67)
This equation may be solved and to first order in the velocity gradient the order parameter is found to be

$$S=S_0+\frac{\sqrt{\mu _2^216S_0^2\mu _1^2}}{8|\alpha _F|}|_yu_x|$$
(68)
and hence the order parameter increases as a consequence of flow alignment. Finally we consider the Miesowicz viscosities in a system with a variable order parameter. The stress tensor evolution equation can be solved for a simple shear flow in the presence of a strongly aligning external field. If the director is given by $`\{n_x,n_y\}`$, the associated Miesowicz viscosity, $`\eta _{eff}`$, is given by

$`\eta _{eff}`$ $`=`$ $`{\displaystyle \frac{1}{4S_0^2}}\{S_0^2[2\alpha _4+\alpha _5+\alpha _6]+S^2[\alpha _3\alpha _2+4\alpha _1n_x^2n_y^2]`$ (69)
$`+SS_0[\alpha _6(2n_x^21)+(\alpha _2+\alpha _3)(n_x^2n_y^2)+\alpha _5(2n_y^21)]\}`$
where $`S`$ is the order parameter in the flow and $`S_0`$ is the equilibrium value of the order parameter defined in equation (63). In the limit that the order parameter is fixed (ie $`S=S_0`$), equation (69) reduces to the standard results

$`\eta _{par}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\alpha _4+\alpha _3+\alpha _6)`$
$`\eta _{perp}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\alpha _4\alpha _2+\alpha _5)`$ (70)
where $`\eta _{par}`$ and $`\eta _{perp}`$ correspond to the director being aligned parallel, $`𝐧=\{1,0\}`$ and perpendicular, $`𝐧=\{0,1\}`$, to the flow. In the presence of both a velocity gradient and a strong aligning field the change in the order parameter is dominated by the aligning field and is hence given to a good approximation by equations (61) and (62).

## 6 Results

A number of simulations were undertaken in order to validate the method described in the previous section. All the results were obtained using a $`D2Q13`$ lattice.

In the absence of flow and in the presence of periodic boundary conditions, the order parameter is found to be consistent with the value predicted by equation (63) to machine accuracy. If an external magnetic field is applied through the introduction of a molecular field term of the form (12), it is found that the director becomes fully aligned with the external field and that the order parameter follows equation (61) with an accuracy of better than $`0.1\%`$. These two results confirm that the free energy correctly controls the dynamical equations through the associated molecular field. In principal, therefore, the effects of temperature could be introduced into the model through the coefficients in the Landau-deGennes free energy.

The angle of alignment of the director in a shear flow in the absence of an external magnetic field is given by equation (65). In order to test that the technique recovered this result, simulations were run with the free energy parameters, $`(\alpha _F=0.0512,\beta _F=0,\gamma _F=0.01)`$, chosen to give an equilibrium order parameter of $`0.8`$ and the parameters, ($`L_1=0.001,L_2=0`$), selected to give nematic elastic behaviour equivalent to the one constant approximation. The LB parameters were chosen to be, $`(\tau _0=1.1000,\eta _0=0.3036,\eta _2=0.1565,\mu _1=0.3823,s\mu _2=1.261,\rho =1.8)`$. If the order parameter $`S_0`$ is assumed to be 0.8, these values recover the viscosity ratios, $`(\alpha _2/\alpha _4=0.9556,\alpha _3/\alpha _4=0.0144,\alpha _5/\alpha _4=0.5565,\alpha _6/\alpha _4=0.4135)`$, which are consistent with MBBA at 25<sup>o</sup> C

With this choice of parameters, the angle of rotation of the director is predicted to be $`7^o`$ and this is recovered to an accuracy of less than $`1\%`$. The parameters $`\mu _1`$ and $`\mu _2`$ were then adjusted to give a range of values of angle of alignment and the results are shown in Figure (1) together with the expected curve, equation (65), which is shown as a continuous curve with the value of the order parameter taken from the simulation data. The change in the order parameter predicted by equation (68) is also recovered to within less than $`0.3\%`$.

As a final test of the algorithm, the Miesowicz viscosities were measured for simulations in which the director orientation was controlled by a strong magnetic field. Using the parameters appropriate for MBBA the expected ratio of the viscosity with the director aligned parallel and perpendicular to the direction of flow is 5.03 and this value was recovered with an error of $`0.1\%`$. The expected form of the viscosity is given by equation (69) and it can be shown that this implies that the viscosity should be a linear function of $`1/\tau _0`$. Figure (2) shows that the expected linearity is observed within the simulations; the slope to intercept ratio for each curve is in error by less than $`0.2\%`$

## 7 Conclusions

In this paper we have presented a generalised lattice Boltzmann scheme which recovers the tensor order parameter equations for a nematic liquid crystal proposed by Qian and Sheng . The generalised method is based on a single tensor density from which all the macroscopic quantities can be recovered. A Chapman Enskog analysis demonstrates that the algorithm recovers the target macroscopic equations and test simulations demonstrate that the method correctly recovers the evolution of the director, the order parameter and the velocity gradients in the presence of shear and magnetic fields.

The method will generalise straightforwardly to three dimensions and work is currently in progress to extend the scheme to model a mixture of a isotropic and nematic fluids.

Acknowledgements  We thank Dr D J Cleaver for valuable conversations throughout the developments described in this paper.


text:
# A comparative study of the X-ray afterglow properties of optically bright and dark GRBs

## 1 Introduction

About 50% of well-localized GRBs show optical transients (OTs) successive to the prompt gamma-ray emission, whereas an X-ray counterpart is present in 90% of cases. It is possible that late and shallow observations could not detect the OTs in some cases; several authors argue that dim and/or rapid decaying transients could bias the determination of the fraction of truly obscure GRBs (Fynbo et al., 2001; Berger et al., 2002). However, recent reanalysis of optical observations (Reichart & Yost, 2001; Ghisellini et al., 2000; Lazzati et al., 2002) has shown that GRBs without OT detection (usually dark GRBs, FOAs Failed Optical Afterglows, or GHOSTs, Gamma ray burst Hiding an Optical Source Transient) have had on average weaker optical counterparts, at least 2 magnitudes in the R band, than GRBs with OTs. Therefore, they appear to constitute a different class of objects, albeit there could be a fraction undetected for bad imaging.

Two hypothesis have been put forward to explain the behaviour of GHOSTs. First, they are similar to the other bright GRBs, except for the fact that their lines of sight pass through large and dusty molecular clouds, that cause high absorption. Second, they are more distant than GRBs with OT, at $`z5`$ (Fruchter, 1999), so that the Lyman break is redshifted into the optical band. These GRBs might be associated with the explosion of ancient Population III, high mass stars. Nevertheless, the distances of a few dark GRBs have been determined and they do not imply high redshifts (Djorgovski et al., 2001b; Antonelli et al., 2000; Piro et al., 2002)

Goal of this paper is an analysis of a complete sample of BeppoSAX X-ray afterglows in order to distinguish between these various scenarios, including all x-ray fast observation from the launch to February 2001. In §2 and §3 we present the data analysis of the afterglows and we show the results, whose implications are discussed in §4. Finally, we summarize our conclusions in §5.

## 2  Data Analysis

We have analyzed all the 31 fast BeppoSAX observations of GRB X-ray afterglows taken by the Low Energy (0.1 - 10 keV) and Medium Energy (1.6 - 10 keV) Concentrator Spectrometer (LECS and MECS respectively, see Parmar et al. 1997, Boella et al. 1997) up to GRB010222. We excluded only GRB960720 for the late follow up, GRB990705 due to its high contamination of a nearby X-ray source, and GRB980425 due to its peculiarity. X-ray follow up observations usually start $`910`$ hours after the high energy event and the typical observation time is $`2\times 10^5`$ seconds for MECS and $`5\times 10^4`$ for LECS. The exposure - or integration - lasts $`1/3`$ of the observation.

In order to find out the GRB X-ray afterglows, we first have built up the images of each GRB with the MECS and selected sources with $`3\sigma `$ significance within the WFC error box. Successively, we have built the light curves of these sources to recognize afterglows through their typical fading emission. The counts have been collected within a circle centered on source with radius $`r=4`$ arcmin. Then we have subtracted the background collected in annuli around the extraction area and 5 times more extended. Local background have been used in order to take into account possible time fluctuations. TOOs successive to the first one (typically $`2`$ day after) have been used, if available. We have fitted the light curves with simple power law $`N_{cts}t^\delta `$ (where $`N_{cts}`$ is counts per second) and 26 sources with decaying index $`\delta >0`$ (at $`90\%`$ Confidence Level) have been recognized as GRB afterglows. In the case of GRB970111, GRB991106 and GRB000615 we have detected 1 source within WFC error box that does not show a significant fading behaviour. We will refer to them as ”candidate” afterglows111In the case of GRB991106, the source in the WFC error box could be a type-I Galactic X-ray burst (Cornelisse et al., 2002).. We have calculated the probability to have serendipitous sources with their flux within the WFC error box, adopting the Log N - Log S distribution for BeppoSAX released by Giommi et al. (2000). The probability is $`0.027`$ for each one, while the probability that all of them are not afterglows is $`P10^5`$.

The MECS integration time for GRB990907 was only 1070 seconds, so that the presence of a fading flux could not be verified. The X-ray source detected was recognized as the GRB afterglow because the probability to have a serendipitous source with flux $`10^{12}`$ erg cm-2 sec-1 (see further) in the WFC error box was $`0.007`$.

Finally, in the case of GRB990217 and GRB010220, we have not detected any source with $`3\sigma `$ significance.

To obtain flux, we have produced spectra for the afterglows from LECS and MECS first TOO data. For absorption and spectral index, we have selected those with more than 150 photons in the MECS ( background subtracted ). 5 GHOSTs and 9 OTGRBs passed this criterion.

We have generally taken the LECS data between 0.1 and 4.0 keV and the MECS data between 1.6 and 10 keV. The backgrounds we have used are the library ones because they have a very good signal-to-noise ratio, due to long exposition222 In the case of GRB970111, 970402 and 991014 the use of local background enabled us to gather better results.. However, we have taken the minimum energy for LECS to be 0.4 keV if the Galactic column density was $`N_H5\times 10^{20}`$ cm-2 because in this case the low energy backgrounds differ from the library ones, which have been taken at high Galactic latitudes and lower column densities (Stratta et al., 2002). If we had not adopted this criterion, our analysis would have led to overestimate the true absorption at the source.

The standard spectrum model to fit the data consists of a constant, Galactic absorption, extragalactic absorption (i.e. in situ) and a power law. The constant has been included because LECS and MECS observe a decaying source at different times. Its value is allowed to vary within a range, obtained in each case by fitting LECS and MECS data in the 1.6 - 4 keV interval (to avoid absorption effects) with a simple power law model. The redshift in our fits has been forced to be 1 for all bursts. This value corresponds roughly to the average redshift of OTGRBs. We have adopted this ”working hypothesis” to obtain a homogeneous set that allows us to compare the absorption properties of dark GRBs in the assumption that they are at the same distance.

We have calculated the 1.6 - 10 keV flux of dark and bright GRBs 11 hours after the burst trigger. We have chosen this time to avoid effects of changes in decaying slope. The average count rate in the MECS has been associated with the flux given by the spectrum. Successively, we have taken the count rate at 11 hours, which is given by light curves, to compute the flux at that time. In most cases, observations include it. In a few cases (e.g. GRB000926) the flux has been extrapolated.

For GRB990907, the counts collected were very few and we have not been able to do any spectral analysis. We have estimated the flux assuming a spectral index $`\alpha =1.05`$. For the two non-detections, we have calculated the $`3\sigma `$ upper limits on counts and converted them to flux adopting again $`\alpha =1.05`$. In all successive analysis, upper limits have been included as true afterglows as well as candidate afterglows.

As a first assessment of our study, we can say that X-ray afterglows follow the prompt gamma emission in 26 of 31 cases, which constitute 84% of the sample. If all doubtful sources are considered as afterglows, then the fraction of X-ray afterglows increases up to $`94\%`$. Instead, optical afterglows are 11 and constitute only $`37\%`$333 GRB980515 has not been included in this calculation, see further.  of the sample. We note that all these fractions are in agreement with published data.

We do not know any optical study on GRB980515. We have calculated its X-ray flux but this burst has not been included in our successive analysis.

## 3  The X-ray Spectral and Flux Properties

The data we have obtained are the result of the convolution of the intrinsic distribution with the measurement error distribution. Under the assumption that both are gaussian, it is possible to deconvolve the two distribution. We have followed a maximum likelihood method (Maccacaro et al., 1988) to gather jointly the best estimates of parent distribution mean and standard deviation. We have used these best estimates (hereafter, indicated with index m) for successive analysis, but we have calculated and shown also the weighted mean and standard deviation of our data. The complete set of fit parameters is given in Table 1 and plotted in Figures 1 and 2.

For GHOSTs, the weighted mean and the standard deviation of the measured energy indexes are $`\alpha =1.3\pm 0.17`$ (hereafter errors are at $`1\sigma `$, unless otherwise indicated) and $`\sigma =0.3`$ respectively. The best estimates for the parent population are $`\alpha ^m=1.3_{0.26}^{+0.27}`$, $`\sigma ^m=0^{+0.4}`$ 444In a few cases, the best estimates of the standard deviation in the parent population are equal to or compatible with zero. This suggests that measurements are dominated by experimental errors.. In the case of OTGRBs, $`\alpha =1.04\pm 0.03`$, $`\sigma =0.44`$ and $`\alpha ^m=1.05_{0.06}^{+0.11}`$, $`\sigma ^m=0.05_{0.05}^{+0.13}`$ for the observed and the parent distribution respectively. Energy indexes of dark and optically bright burst are compatible at $`1\sigma `$ level.

The mean value and the standard (linear) deviation of the measured absorption (hereafter, in units of $`10^{22}\mathrm{cm}^2`$) are, respectively, $`N_H=0.13_{0.13}^{+0.42}`$, $`\sigma =3`$ for dark GRBs, and $`N_H=0.13\pm 0.06`$, $`\sigma =1.7`$ for OTGRBs. The best estimates for the parent population are $`N_H^m=0.14_{0.14}^{+1.46}`$, $`\sigma ^m=0^{+1.58}`$ for dark GRBs, and $`N_H^m=0.13_{0.075}^{+0.13}`$, $`\sigma ^m=0^{+0.35}`$ (see also Stratta et al. 2002) for OTGRBs. The amount of absorption does not appear statistically different for optically bright and dark GRBs in the assumption that they lie at the same average $`z`$.

The logarithmic weighted means and the standard deviations of the observed X-ray fluxes (c.g.s. units) are $`<`$ log $`F`$ $`>`$ $`=12.38\pm 0.02`$ , $`\sigma =0.34`$ for dark GRBs and
$`<`$log $`F`$ $`>`$ $`=11.45\pm 0.01`$ , $`\sigma =0.65`$ for OTGRBs. Best estimates for the parent population are $`<`$log $`F`$ $`>^m`$ = $`12.53_{0.09}^{+0.11}`$, $`\sigma ^m=0.23_{0.05}^{+0.09}`$ for dark GRBs and $`<`$ log $`F>^m`$ = $`11.85_{0.23}^{+0.22}`$, $`\sigma ^m=0.47_{0.12}^{+0.2}`$ for GRBs with OTs. The GHOST mean flux is likely overestimated, because we have considered upper limits as detections.

The logarithm of ratio between the mean fluxes of the two parent populations is $`0.68\pm 0.25`$, which corresponds to 4.8 in linear units. A K-S test performed on the flux distributions shows that the probability that optically bright and dark GRBs derive from the same population is $`P=2\times 10^3`$. This is a conservative result, because it has been obtained by including the upper limits and the non-fading sources as true afterglows in the set of dark GRBs. If we were to substitute the the non-fading source fluxes with the $`3\sigma `$ upper limits of their WFC error boxes, then the distributions of dark and optically bright GRBs would be even more different because limits are lower.

## 4 Discussion

Our analysis shows that dark GRBs have on average weaker X-ray flux than bright GRBs. Then, we could simply explain why we miss their optical detection by assuming that dark bursts are weaker than OT GRBs in the optical band by the same ratio. Dark bursts should have had OTs at least 2 magnitudes fainter than OT GRBs; the 4.8 flux ratio that we have found corresponds to $`1.7`$ magnitudes.

In order to check the viability of this hypothesis, we have calculated the optical flux density in the R band and hence the optical-to-X flux ratios (hereafter $`f_{OX}`$) of each OTGRB and GHOST 11 hours after the burst (Lazzati et al. 2002, Fynbo et al. 2001 and reference therein; Vreeswijck et al. 1999; Masetti et al. 2001). Upper limits on optical fluxes of GHOSTs have been extrapolated from the tightest constraint available and adopting an optical flux decaying index $`\delta =1.15`$. Results are shown in Table 1 and plotted in Figures 3, 4 and 5. We note that the optical and X ray fluxes of OTGRBs are correlated: the higher is the the X-ray flux, the more luminous is the optical counterpart. The probability that it occurs by chance is only $`1.5\%`$. The logarithmic standard deviation of $`f_{OX}`$’s is $`\sigma _{f_{OX}}=0.42_{0.12}^{+0.2}`$, which corresponds to a multiplicative factor of 2.6, while the logarithmic mean is $`<`$ log $`f_{OX}^{OT}`$ $`>^m`$ = $`0.3\pm 0.22`$ if the X-ray and optical fluxes are expressed in $`10^{13}`$ erg cm-2 sec-1 and $`\mu `$Jy respectively. We have fitted the distribution of X-ray and optical fluxes with the function log $`F_{optical}`$ = K + $`A`$ log(F1.6-10keV). The best fit values are $`A=0.81`$, K = $`0.41`$. We have also calculated the average optical-to-x spectral index, $`\overline{\alpha }_{OX}`$, as function of $`<`$ log $`\overline{f}_{OX}^{OT}`$ $`>^m`$, by adopting the X-ray and optical density flux at 2 keV and R band respectively and X-ray spectral index $`\alpha _X=1.05`$. Our result is $`\overline{\alpha }_{OX}^{OT}=0.79\pm 0.054`$.

If we exclude GRB980519, which seems to be the only afterglow explained by interaction of a jet outflow with a star wind medium (Jaunsen et al., 2000), the correlation is strengthened: $`\sigma _{f_{OX}}=0.28_{0.08}^{+0.14}`$, which corresponds to a multiplicative factor of 1.9; the probability of a chance occurrence is $`<0.001`$; $`<`$ log $`f_{OX}^{OT}>^m`$ = $`0.18_{0.14}^{+0.16}`$. The best fit values are $`A=0.91`$, $`K=0.24`$.

We can immediately recognize that $`75\%`$ of the GHOSTs of our sample (14 of 19) have optical flux upper limits consistent with OT detections (see Fig. 3), so that they may not be actually ”dark”. Optical follow-ups conducted for these bursts would not have been deep enough to detect the faintest OTs in our set. A similar fraction has been found out by Fynbo et al. (2001) and Berger et al. (2002), comparing sets of non-detections with the light curves of the dim afterglows of GRB000630 and GRB020124. It is worth noting that the $`f_{OX}`$ upper limits of these 14 objects are quite similar to OTGRBs ones: a K-S test performed shows that the probability they belong to the same population is not marginal555GRB010220 has limits on optical and X-ray flux, so its $`f_{OX}`$ is not constrained. However, wherever it is, it would not affect much the result.. This result gives support to the fact that they could be faint sources with optical properties assimilable to X-ray ones.

The remaining 5 objects, which constitute $`25\%`$ of GHOSTs, have optical flux upper limits lower than all OTGRBs in our set (see Fig. 2). Furthermore, their optical emission must have been even 2-3 times fainter than the dim afterglow of GRB000630 and 4-6 fainter than GRB020124 (Berger et al., 2002), which had an R band flux $`F_{11h}^R3.4\mu `$ Jy and $`F_{11h}^R7.9\mu `$ Jy respectively 666 Data extrapolated with best fit values given by the authors and corrected for Galactic extinction.. These two optical afterglows, however, are not the weakest ever occurred. In our set, the OT of GRB980613 is even dimmer (see table 1) and establishes a more stringent test.

We wonder if the hypothesis of weaker flux at all wavelenghts can hold for these 5 GHOSTs (hereafter, we will refer to them as ”the darkest” or the ”most obscure”, etc. for simplicity). If so, we would expect that their X-ray fluxes were proportionally very weak like the optical fluxes, so that their $`f_{OX}`$’s should be not very different from OTGRBs. In 4 cases of 5 the $`f_{OX}`$’s are lower than all OTGRBs and $`410`$ times lower than the average optical-to-x flux ratio of OTGRBs. The exception is GRB990217, which has upper limits both in optical and in X-rays flux. If we use both of them, then we get log $`f_{OX}=0.2`$ which is much more similar to $`<`$ log $`f_{OX}^{OT}`$ $`>^m`$. We have performed a K-S test on $`f_{OX}`$ between all these darkest bursts and all OTGRBs. The probability that they are drawn from the same distribution is $`P0.01`$ ($`2.6\sigma `$ confidence level). The average optical-to-X spectral index of these objects $`\overline{\alpha }_{OX}0.62`$, well below that of OTGRBs. Therefore, we have a strong indication that for these bursts the spectrum is depleted in the optical band, by $`2`$ magnitude on average.

The x-ray mean flux of these 5 GHOSTs is $`12.48\pm 0.16`$. The logarithmic ratio between the OTGRB X-ray flux mean and this mean is log $`r=0.63\pm 0.28`$, which corresponds to a 4.3 factor. A hypothesis for the absence of OT and fainter X-ray flux is that of very high redshift. The GRB prompt emission and X-ray afterglow of the strongest bursts (e.g. GRB990123 and GRB990510) could be detectable even they occur at $`z>10`$ (Lamb & Reichart, 2000a). However, if GHOSTs were at $`z5`$, then extragalactic hydrogen clouds would entirely wash out optical emission (Piro, 2002; Fruchter, 1999; Becker et al., 2001).

To estimate the average redshift of the most obscure GRBs, we shall use the formula (Lamb & Reichart, 2000a):

$$F(\nu ,t)=\frac{L_\nu (\nu ,t)}{4\pi D^2(z)(1+z)^{1\alpha +\delta }}$$
(1)
where $`\alpha `$ is the spectral index, $`\delta `$ is the decaying (temporal) index and $`D(z)`$ is the comoving distance. We assume the cosmological parameters values $`H_0=65\mathrm{k}\mathrm{m}\mathrm{sec}^1\mathrm{Mpc}^1`$, $`\mathrm{\Omega }_M=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. The average of the known redshifts of OTGRBs in our set is $`\overline{z}_{OT}1.5`$.

In the simplest model of GRB afterglows, $`\delta =4/3,\alpha =2\delta /3`$, so $`1\alpha +\delta =5/9`$. For such parameters, the average redshift of the darkest GRBs should be $`2.6\overline{z}_D8.7`$ under the assumption that the lower mean flux were only due to their larger distances and not to an intrinsic difference in their luminosity. Using the best estimate of $`\alpha =1.05`$ calculated for OTGRBs and the average $`\delta =1.33`$ of the strongest bursts of our sample, we obtain $`2.3\overline{z}_D7.8`$. We should also expect a distribution of burst redshifts around $`\overline{z}_D`$. These facts make the high redshift scenario for the most obscure GHOSTs still plausible.

Adopting the hypothesis GRBs are the final result of very massive star evolution, an interesting issue to address is what might be the progenitors of GRBs at very high redshifts.

Currently, we observe only old and low-mass Population II stars, but even high mass stars could have formed. Theories suggest that the first stars of the Universe - the so-called Population III - might have very large mass, so they could possibly be good candidates. Recent calculations suggest (Lamb & Reichart, 2000b; Valageas & Silk, 1999; Gnedin & Ostriker, 1997; Ostriker et al., 1996) that the star formation rate has two peaks. The first one, at $`20z16`$ is due to Population III stars. The second one, due to Population II, is higher and much broader and it is at a redshift in the range $`12z2`$. Also the number of stars (i.e. the star formation rate time-dilated and weighted by the comoving volume of the universe) shows two peaks at $`z8`$ and $`z2`$.

In a few cases, however, the redshifts of some dark GRBs have been almost securely found, e.g. GRB970828 at $`z=1`$ (Djorgovski et al., 2001b) , GRB 000210 at $`z=0.85`$ (Piro et al., 2002), GRB 000214 at $`z=0.44`$ (Antonelli et al., 2000), while GRB981226 is also likely to have not occurred at very high redshift (Frail et al., 1999), because the candidate host galaxy is still detected in the R band. With present statistics, at least $`15\%`$ of the examined dark GRBs are not at very large redshift. It should be noticed that two of them are included in the list of most obscure objects in our set.

An hypothesis to explain the lack of the optical emission, alternative to the very high redshift scenario, may be strong absorption (Djorgovski et al., 2001a). So far, we have collected many indications that GRBs take place in dense environments, like the Giant Molecular Clouds (hereafter GMCs) (Piro, 2002). GMCs are very rich in dust, which extinguishes very efficiently the optical and UV light. Piro et al. (2002) argues that in the case of GRB000210 the lower limit on amount of obscuration is 1.6 magnitude in the R band. This value has been obtained extrapolating a power law spectrum, described by the fireball model, from the X-ray band to the optical band and comparing the expected flux with the upper limits. We find a similar result through our model-independent analysis of optical-to-x flux ratios. The $`f_{OX}`$ upper limits of the burst is 3.8 times lower than $`\overline{f}_{OX}^{OT}`$, which corresponds to $`1.5`$ magnitude depletion. The measurements of Chandra X-ray Observatory showed that the amount of local absorbtion is able to explain this obscuration, under the assumption that the dust-to-gas ratio of the intervening medium is the same of the Galaxy. However, we note that in the case of OTGRBs the dust-to-gas ratio seems not to be consistent with the Galactic one (Stratta et al., 2002; Galama & Wijers, 2001). Similarly, Djorgovski et al. (2001b) derived extinction in the case of GHOST GRB970828 (Yoshida et al., 2001), for which a significative amount of X-ray absorption was detected.

If the most obscure GHOSTs were similar to GRBs with OT except for higher absorption, we would expect to see differences in values of $`N_H`$. From our results, we cannot affirm that $`N_H`$ in these bursts shows this tendency, also due to considerable errors (see Table 1 and Figure 1). For those with good statistics, we do not find any absorption value $`3\sigma `$ higher than the Galactic value but marginal evidence ($`2\sigma `$). On the other hand, we cannot rule out the hypothesis of obscuring GMCs altogether. The upper limits on $`N_H`$, a few $`\times 10^{22}`$, are in fact the typical column densities of GMCs. The optical absorption, however, does not imply that most obscure GHOSTs have X-ray flux weaker than OTGRBs, as we have found in our analysis, because X-ray absorption is almost negligible at energy larger than 1.6 keV. Reichart & Yost (2001) try to reconcile this fact with the hypothesis of dusty birthplaces for GRBs and, in particular, they considered the effect of variously beamed GRB fireballs on their dusty environments. The energetics of GRBs are more or less the same for all events (Frail et al., 2001), but the beaming angles differ, being narrower for stronger bursts. The larger the beaming angle is, the more difficult it is for the diluited prompt UV and X-ray emission to destroy dust along the line of sight (Waxman & Draine 2000; Fruchter et al. 2001; Draine & Hao 2002; Perna & Lazzati 2002), so that we see a weak GRBs without OT. With a narrower beaming angle, the prompt emission will destroy a larger fraction of dust and the GRBs will appear strong and with OT. If this hypothesis is correct, on the basis of our results we have to assume that the average beaming angle of the darkest GHOSTs is $`2`$ times wider than the OTGRB one. According to Frail et al., the average beaming angle of BeppoSAX OT GRBs is $`\theta 0.1\mathrm{rad}`$, so that the average beaming angle of the darkest GRBs should be $`\theta 0.2\mathrm{rad}`$. This prediction is important, because it can be experimentally tested by observing and timing the presence of achromatic breaks in the light curves.

Another consequence of dark GRB occurrence in high density environments should be the detection of semi-ionized absorber in the low energy X-ray spectrum. So far, this kind of feature has not been found. The ionization front, however, should be rather sharp (see e.g. Draine & Hao 2002), and therefore it would be hard to detect signatures of semi-ionized species in the X-ray spectra of the bursts.

## 5 Conclusions

We have discussed the issue of GRBs with X-ray but no optical afterglows. We have performed a standard temporal and spectral analysis of a complete sample of 31 GRB X-ray follow up observations of BeppoSAX, i.e. all the fast observations from the launch until February 2001. We have found that X-ray afterglows follow the prompt gamma emission in $`84\%94\%`$ of the cases.

We have obtained the 1.6 - 10 keV fluxes 11 hours after the trigger for each GRB and the values of NH at $`z=1`$ to compare the absorption properties for strong X-ray afterglows. While absorption of optically bright and dark GRBs does not appear to be significatively different, the fluxes of GRBs with OT are on average about 5 times stronger than GHOST ones. The probability that GHOSTs and optically bright GRBs belong to the same population in fluxes is $`0.002`$.

From the very fact that X-ray flux of dark GRBs is 5 times lower than OTGRBs, the optical flux could be $`2`$ magnitude lower under the assumption that the shape of the optical-to-X spectrum is the same of OTGRBs. This difference could explain the non-detection of the optical transient. In order to test this hypothesis, we have calculated the optical-to-X flux ratios of OTGRBs and upper limits for GHOSTs. OTGRBs show a tight correlation of optical and X-ray fluxes. The mean for OTGRBs is $`<`$ log $`f_{OX}^{OT}`$ $`>^m`$ = $`0.3\pm 0.22`$ and $`\sigma =0.4`$ ; the probability of a chance correlation is a marginal $`1\%`$. We find that $`75\%`$ of GHOSTs have $`f_{OX}`$ upper limits similar to OTGRB ones; however, the remaining $`25\%`$ of dark bursts are fainter in optical than in X-rays, being their average optical-to-x flux ratio $`<`$ log $`f_{OX}`$ $`>^m`$ $`0.4`$. Thus, we have a strong indication that for these bursts the spectra are different from OTGRBs. This result is significant at $`2.6\sigma `$ level.

Two different interpretations for this effect can be given: 1) location at $`z>5`$ ; 2) higher absorption. In the very high redshift scenario, the optical flux of the sample is extinguished by the intervening Ly-$`\alpha `$ systems, while the X-ray flux lower than OTGRBs is understood in terms of a higher distance.

However, given the fact that some GHOSTs of the sample almost certainly do not lie at very high redshift, we have considered the alternative possibility of occurrence in dusty and dense environments like GMCs. We have not found that these bursts have a higher absorption than optically bright GRBs, but we note that upper limits on $`N_H`$ are consistent with those of giant clouds. In the case of GRB000210 our model-independent analysis has shown a depletion in the optical which is compatible with the X-ray absorption measured by Chandra, assuming a dust-to-gas ratio similar to that of our Galaxy.

In the near future, a key role will be played by fast and deep follow up X-ray and optical observations of GRBs, which will allow us to constrain better their spectral properties. In particular, observations in the IR band are a really important tool because they are less sensitive to dust and to Ly-$`\alpha `$ extinction. They will enable us to investigate dark GRB properties like distance, that is a crucial piece of information to disclose the nature of these objects.

M.D.P. is grateful for support from University ”La Sapienza” grant. BeppoSAX was a joint program of Italian Space Agency (ASI) and Netherlands Agency for Aerospace Programs (NIVR). We warmly thank all members of BeppoSAX Scientific Operation Center and Operation Control Center.