format: text

text:
# Maxwellian gas undergoing a stationary Poiseuille flow in a pipe

## 1 Introduction

The Poiseuille flow, first studied by Hagen and Poiseuille towards the half of the 19th century, is still one of the classical examples in fluid dynamics . It describes the steady flow in a slab or in a pipe under the action of a longitudinal pressure gradient, what produces a longitudinal macroscopic velocity with a typical parabolic profile in the transverse directions. Essentially the same effect is generated when the longitudinal pressure difference is replaced by a uniform external force $`m𝐠`$ (e.g. gravity) directed longitudinally. This latter mechanism for driving the Poiseuille flow does not produce longitudinal gradients and so is more convenient than the former in computer simulations as well as from the theoretical point of view, especially to assess the reliability of the continuum description . The first study of the Poiseuille flow driven by an external force we are aware of was carried out by Kadanoff et al. , who simulated the flow with the FHP lattice gas automaton  to confirm the validity of a hydrodynamic description for lattice gas automata. In the context of a dilute gas, Esposito et al.  studied the Boltzmann equation and found that if the force is sufficiently weak there is a solution which converges, in the hydrodynamic limit, to the local equilibrium distribution with parameters given by the stationary solution of the Navier–Stokes (NS) equations. A generalized Navier–Stokes theory was seen to give a reasonable account of a fluid composed of molecules with spin when compared with molecular dynamics simulations .

The first kinetic theory analysis of the Poiseuille flow clearly exhibiting non-Newtonian behavior was carried out by Alaoui and Santos , who found an exact solution of the Bhatnagar–Gross–Krook (BGK) model kinetic equation for a particular value of $`g`$. The general solution of the BGK model under the form of an expansion in powers of $`g`$ through fifth order was obtained by Tij and Santos . The most interesting outcome of the solution, as first recognized by Malek Mansour et al. , was the presence of a positive quadratic term in the temperature profile to second order in $`g`$, in addition to the negative quartic term predicted by the NS description. As a consequence of this new term, the temperature profile exhibits a bimodal shape, namely a local minimum at the middle of the channel surrounded by two symmetric maxima at a distance of a few mean free paths. In contrast, the continuum hydrodynamic equations predict a temperature profile with a (flat) maximum at the middle. The Fourier law is dramatically violated since in the slab enclosed by the two maxima the transverse component of the heat flux is parallel (rather than anti-parallel) to the thermal gradient. This phenomenon is not an artifact of the BGK model since the same results, except for the numerical values of the coefficients, were derived from an exact solution of the Boltzmann equation for Maxwell molecules through second order , as well as from approximate solutions of the Boltzmann equation for hard spheres by Grad’s moment method . It is interesting to note that this correction to the NS temperature profile is not captured by the next hydrodynamic description, namely the Burnett equations . Therefore, in order to describe a non-monotonic temperature field \[which is an $`𝒪(g^2)`$-order effect\] from a hydrodynamic description one needs to go at least to the super-Burnett equations . A systematic asymptotic analysis of the BGK equation with diffuse boundary conditions to second order in the Knudsen number  validates the results of Ref.  as corresponding to the Taylor series expansion of the normal (or Hilbert) solution around the center of the gap. The theoretical predictions of a bimodal temperature profile have been confirmed at a qualitative and semi-quantitative level by numerical Monte Carlo simulations of the Boltzmann equation  and by molecular dynamics simulations . In the case of a dense gas, even though the non-monotonic behavior of the temperature profile may disappear, the existence of a quadratic term coexisting with the classical quartic term is supported by molecular dynamics simulations . It is worth mentioning that when the Poiseuille flow is driven by a longitudinal pressure gradient rather than by an external force, the NS description is in much better agreement with Monte Carlo simulations of the Boltzmann equation .

Practically all the studies about the Poiseuille flow driven by an external force have considered the planar geometry , i.e. the fluid is assumed to be enclosed between two parallel, infinite plates orthogonal to the $`y`$-direction, the force acting along the $`z`$-direction. In the stationary laminar flow the physical quantities have a unidirectional dependence on the $`y`$ coordinate only. While the planar geometry is the simplest one to study the Poiseuille flow, it is much more realistic to consider that the fluid is enclosed in a cylindrical pipe, the external force being directed along the symmetry axis $`z`$. In that case, the quantities depend on the distance $`r=\sqrt{x^2+y^2}`$ from the axis and the fluxes have in general both radial and tangential components. Two of us have recently solved the BGK model for the cylindrical Poiseuille flow through fourth order in $`g`$ . The results showed that the structure of the hydrodynamic and flux profiles in the pipe geometry is similar to that of the slab geometry. In particular, the temperature exhibits a non-monotonic behavior as one moves apart from the pipe axis. On the other hand, the results are quantitatively sensitive to the geometry of the problem. For instance, the normal stress along the gradient direction is uniform in the planar case ($`P_{yy}=\text{const}`$), while it has a non-trivial spatial dependence in the cylindrical case ($`P_{rr}\text{const}`$). Also, comparison with the BGK solution for the planar geometry  shows that the differences between the kinetic theory results and the NS predictions are in general more pronounced in the planar case than in the cylindrical case.

As is well known, in the BGK model kinetic equation the complicated nonlinear structure of the Boltzmann collision operator is replaced by a single relaxation-time term to the local equilibrium distribution . This is usually sufficient to account for many of the nonequilibrium properties of the underlying Boltzmann equation, at least at a qualitative level . However, the existence of only one effective collision frequency does not allow the BGK model to describe quantitatively those states where momentum and heat fluxes coexist and are inextricably intertwined, as happens in the Poiseuille flow. This limitation of the BGK model is responsible for an incorrect value $`\text{Pr}=1`$ of the Prandtl number , in contrast to the Boltzmann value  $`\text{Pr}\frac{2}{3}`$. In the planar Poiseuille flow, it is possible to assess the validity of the BGK solution  by comparing it with the results derived from the Boltzmann equation for Maxwell molecules  and hard spheres . The aim of this paper is to fill the gap caused by the absence of results from the Boltzmann equation for the cylindrical Poiseuille flow. We consider a gas of Maxwell molecules and solve the infinite hierarchy of moment equations stemming from the Boltzmann equation through third order in the external force. The results confirm the correctness of the functional dependence of the hydrodynamic and flux profiles obtained from the BGK model . On the other hand, the numerical values of the coefficients are in general different in both kinetic theories, as expected. The organization of the paper is the following. Section 2 is devoted to the description of the flow under study and its solution in a hydrodynamic description to NS order. The kinetic theory description is presented in Sec. 3. In order to solve the moment hierarchy in a systematic and recursive way, we carry out a perturbation expansion in powers of the external field in Sec. 4, the technical details being relegated to Appendix B. The results are extensively discussed and compared with those of the NS and BGK descriptions in Sec. 5. Finally, the main conclusions of the paper are briefly presented in Sec. 6.

## 2 Navier–Stokes description of the cylindrical Poiseuille flow

The cylindrical Poiseuille flow studied in this paper refers to a monatomic dilute gas enclosed in a long pipe of radius $`R`$. Let us take the $`z`$-axis as the symmetry axis of the cylinder. The gas particles are subjected to the action of a constant external force per unit mass $`𝐠=g\widehat{𝐳}`$ (e.g. gravity) parallel to $`z`$. After a certain transient stage, the system reaches a steady laminar flow. This nonequilibrium state is characterized by gradients of the hydrodynamic variables along the directions $`x`$ and $`y`$ orthogonal to the cylinder axis. By symmetry, the hydrodynamic fields are expected to depend on $`x`$ and $`y`$ through the distance $`r=\sqrt{x^2+y^2}`$ from the $`z`$-axis. In our analysis, we are interested in the bulk region of the flow. This means that the radius $`R`$ of the cylinder is assumed to be large enough (as compared with the mean free path) to allow for the existence of such a region.

The steady-state balance equations expressing the conservation of momentum and energy are

$$𝖯=\rho 𝐠,$$
(1)
$$𝐪+𝖯:𝐮=0,$$
(2)
where $`\rho `$ is the mass density, $`𝐮=u_z\widehat{𝐳}`$ is the flow velocity, $`𝖯`$ is the pressure tensor and $`𝐪`$ is the heat flux vector. The geometry of the problem suggests the use of cylindrical coordinates (see Appendix A). Using Eqs. (86)–(90), Eqs. (1) and (2) yield

$$\frac{}{r}\left(rP_{rr}\right)=P_{\varphi \varphi },$$
(3)
$$\frac{1}{r}\frac{}{r}\left(rP_{rz}\right)=\rho g,$$
(4)
$$P_{rz}r\frac{u_z}{r}+\frac{}{r}\left(rq_r\right)=0,$$
(5)
Equations (3)–(5) are exact, but they do not constitute a closed set.

In the Navier–Stokes (NS) description, the momentum and heat fluxes are assumed to be linear functions of the hydrodynamic gradients . In the geometry of the cylindrical Poiseuille flow the NS constitutive equations are written as

$$P_{rr}=P_{\varphi \varphi }=P_{zz}=p,$$
(6)
$$P_{rz}=\eta \frac{u_z}{r},$$
(7)
$$q_r=\kappa \frac{T}{r},$$
(8)
$$q_z=0.$$
(9)
In Eq. (6), $`p=\frac{1}{3}\text{tr}𝖯`$ is the hydrostatic pressure. Equation (7) is Newton’s friction law, $`\eta `$ being the shear viscosity. In addition, Eq. (8) is Fourier’s law, where $`\kappa `$ is the thermal conductivity and $`T`$ is the temperature. The latter is related to the number density $`n=\rho /m`$ (where $`m`$ is the mass of a particle) and the pressure $`p`$ by the (local) equilibrium equation of state. In particular, for a dilute gas, $`p=nk_BT`$, $`k_B`$ being the Boltzmann constant. Equation (6) implies the absence of normal stress differences in a sheared Newtonian fluid, while Eq. (9) means that the heat flux is parallel to the thermal gradient in a fluid obeying Fourier’s law. The combination of Eqs. (3)–(8) gives the following closed set of hydrodynamic equations:

$$\frac{p}{r}=0,$$
(10)
$$r^1\frac{}{r}\left(r\eta \frac{u_z}{r}\right)=\rho g,$$
(11)
$$\frac{}{r}\left(r\kappa \frac{T}{r}\right)=\eta r\left(\frac{u_z}{r}\right)^2.$$
(12)
In principle, the NS hydrodynamic equations (10)–(12) are too complicated to obtain its explicit solution for arbitrary $`g`$. On the other hand, the solution can be found as a series expansion in powers of $`g`$. To third order the result is

$$p(r)=p_0,$$
(13)
$$u_z(r)=u_0\frac{\rho _0g}{4\eta _0}r^2\left(1+\frac{1}{144}\frac{\rho _0^2g^2}{\eta _0\kappa _0T_0}r^4\right)+𝒪(g^5),$$
(14)
$$T(r)=T_0\frac{1}{64}\frac{\rho _0^2g^2}{\eta _0\kappa _0}r^4+𝒪(g^4),$$
(15)
where we have particularized to a dilute gas of Maxwell molecules, in which case $`\eta T`$ and $`\kappa T`$ . The subscript 0 in Eqs. (13)–(15) denotes quantities evaluated at $`r=0`$, i.e. $`u_0=u_z(0)`$, $`T_0=T(0)`$, $`\eta _0=\eta (T_0)`$, …. The corresponding expressions for the fluxes are

$$P_{rz}(r)=\frac{\rho _0g}{2}r\left(1+\frac{1}{192}\frac{\rho _0^2g^2}{\eta _0\kappa _0T_0}r^4\right)+𝒪(g^5),$$
(16)
$$q_r(r)=\frac{1}{16}\frac{\rho _0^2g^2}{\eta _0}r^3+𝒪(g^4),$$
(17)
complemented with Eqs. (6) and (9). It is interesting to note that the spatial variable $`r`$ can be eliminated between Eqs. (14) and (15) to obtain the following nonequilibrium “equation of state”:

$$T=T_0\frac{\eta _0}{4\kappa _0}(u_zu_0)^2+𝒪(g^4).$$
(18)
By equation of state we mean in this context a relationship holding locally among the hydrodynamic fields ($`p`$, $`u_z`$ and $`T`$) and that does not contain $`g`$ explicitly (at least up to a certain order). Since the pressure is uniform in the NS description, it does not enter into Eq. (18). Equation (18) shows that in the Poiseuille flow the velocity and temperature profiles are closely related.

From Eqs. (13)–(15) we can obtain the mass rate of flow in the NS description. Let us consider a circular section of radius $`a`$ centered on the cylinder axis and orthogonal to it. The mass of fluid flowing across this surface per unit time is defined, in the reference frame moving with the fluid at $`r=a`$, as

$$\dot{M}(a)=2\pi _0^a𝑑rr\rho (r)\left[u_z(r)u_z(a)\right].$$
(19)
Using Eqs. (13)–(15) and taking into account that $`\rho =mp/k_BT`$, one gets

$$\dot{M}(a)=\frac{\pi \rho _0^2ga^4}{8\eta _0}\left(1+\frac{5}{384}\frac{\rho _0^2g^2}{\eta _0\kappa _0T_0}a^4\right)+𝒪(g^5).$$
(20)
Before closing this Section, let us make a few remarks about our use of the term “Navier–Stokes solution”. By that we mean the solution to the coupled set of hydrodynamic equations (10)–(12). In fluid mechanics, on the other hand, it is usual to adopt a more restrictive point of view in which by Navier–Stokes equation one means the momentum equation with a constant viscosity $`\eta `$ and a constant density $`\rho `$ (“incompressible fluid”). In that case, the solution of Eq. (11) would simply be

$$u_z(r)=u_0\frac{\rho g}{4\eta }r^2,$$
(21)
regardless of the value of $`g`$. Accordingly, the mass rate of flow would be

$$\dot{M}(a)=\frac{\pi \rho ^2ga^4}{8\eta }.$$
(22)
Equation (21) gives the typical parabolic velocity profile under Poiseuille flow, while Eq. (22) is known as Poiseuille’s law . Nevertheless, the shear viscosity $`\eta `$ depends in general on the temperature $`T`$ and so the assumption of a constant viscosity implicitly implies a constant temperature. The same conclusion follows from Eq. (13) and the incompressibility assumption $`\rho =\text{const}`$. The energy equation (12), however, shows that a velocity profile is compatible with a constant temperature only if the thermal conductivity $`\kappa `$ is allowed to be infinite. As a matter of fact, Eqs. (14) and (20) reduce to Eqs. (21) and (22), respectively, if we formally set $`\kappa _0\mathrm{}`$. Since in a dilute gas the relative strength of the thermal conductivity is comparable to that of the shear viscosity (as measured by the Prandtl number $`\text{Pr}5k_B\eta /2m\kappa \frac{2}{3}`$), the energy equation (12) cannot be neglected and hence $`T\text{const}`$.

In summary, in the remainder of this paper we will refer to Eqs. (10)–(12) as Navier–Stokes equations and will compare their solution for a dilute gas of Maxwell molecules \[Eqs. (13)–(15)\], along with the momentum and energy fluxes \[Eqs. (6), 9), (16) and (17)\], with the results derived from the Boltzmann equation.

## 3 Kinetic theory description

In kinetic theory, the relevant information about the nonequilibrium state of a dilute gas is contained in the one-body distribution function $`f(𝐫,𝐯;t)`$ . Its temporal evolution is governed by the non-linear Boltzmann equation , which in standard notation reads

$`{\displaystyle \frac{}{t}}f+𝐯f+𝐠{\displaystyle \frac{}{𝐯}}f`$ $`=`$ $`{\displaystyle 𝑑𝐯_1𝑑\mathrm{\Omega }|𝐯𝐯_1|\sigma (|𝐯𝐯_1|,\theta )[f^{}f_1^{}ff_1]}`$ (23)
$``$ $`J[f,f].`$
The influence of the interaction potential appears through the dependence of the differential cross section $`\sigma `$ on the relative velocity $`|𝐯𝐯_1|`$ and the scattering angle $`\theta `$. In particular, for Maxwell molecules $`\sigma (|𝐯𝐯_1|,\theta )|𝐯𝐯_1|^1`$.

The knowledge of $`f`$ in a dilute gas is sufficient to determine its macroscopic state. For example, the hydrodynamic variables and their fluxes are just velocity moments of the distribution function. Thus, the local number density $`n`$ is defined by

$$n=𝑑𝐯f.$$
(24)
Analogously, the flow velocity $`𝐮`$ and the hydrostatic pressure $`p`$ are given by the following expressions:

$$𝐮=\frac{1}{n}𝑑𝐯𝐯f,$$
(25)
$$p=\frac{m}{3}𝑑𝐯V^2f,$$
(26)
where the peculiar velocity $`𝐕=𝐯𝐮`$ has been introduced as the velocity of a particle relative to the flow velocity. The momentum and energy fluxes are given by the pressure tensor $`𝖯`$ and the heat flux vector $`𝐪`$, respectively, as follows:

$$𝖯=m𝑑𝐯\mathrm{𝐕𝐕}f,$$
(27)
$$𝐪=\frac{m}{2}𝑑𝐯V^2𝐕f.$$
(28)
In the case of the stationary Poiseuille flow in a pipe, the Boltzmann equation (23) becomes

$$v_x\frac{f}{x}+v_y\frac{f}{y}+g\frac{f}{v_z}=J[f,f].$$
(29)
This equation is invariant under the transformations

$$(x,v_x)(x,v_x),(y,v_y)(y,v_y),(v_z,g)(v_z,g).$$
(30)
Strictly speaking, Eq. (29) must be supplemented by the appropriate boundary conditions at $`r=R`$ describing the interaction of the particles with the cylinder surface. However, in this work we are interested in the bulk region of the system ($`0<r<R\delta `$, where $`\delta `$ is the width of the boundary layer), which is not affected by the details of the boundary conditions. This implies that the radius $`R`$ of the cylinder must be sufficiently large, compared to the mean free path of a particle, to allow the existence of such a region (small Knudsen number). Here we will assume that this is the case and will look for the Hilbert-class or normal solution to Eq. (29), namely a solution where all the spatial dependence of the distribution function takes place through a functional dependence of $`f`$ on the hydrodynamic fields $`n`$, $`𝐮`$ and $`T`$.

Since the hydrodynamic variables of the gas and the associated fluxes are the first few moments of the function $`f`$, it is convenient to consider the hierarchy of moment equations stemming from the Boltzmann equation (29). A moment of an arbitrary degree $`k=k_1+k_2+k_3`$ is defined as

$$M_{k_1,k_2,k_3}(x,y)=𝑑𝐯V_x^{k_1}V_y^{k_2}V_z^{k_3}f(x,y,𝐯).$$
(31)
The invariance properties (30) imply the parity relations

$`M_{k_1,k_2,k_3}(x,y;g)`$ $`=`$ $`(1)^{k_1}M_{k_1,k_2,k_3}(x,y;g)`$ (32)
$`=`$ $`(1)^{k_2}M_{k_1,k_2,k_3}(x,y;g)`$
$`=`$ $`(1)^{k_3}M_{k_1,k_2,k_3}(x,y;g).`$
Although the problem at hand calls for the use of cylindrical coordinates, we will use Cartesian coordinates to construct the moment hierarchy. The reason is two-fold. First, the collision integrals are easier to express in Cartesian (or even spherical) coordinates than in cylindrical coordinates. Second, the use of Cartesian coordinates will serve us as a test of the calculations since, when obtaining the cylindrical components according to Eqs. (79) and (80), the relevant quantities must depend on $`x`$ and $`y`$ through the variable $`r=\sqrt{x^2+y^2}`$.

The Boltzmann hierarchy of moments is obtained by multiplying both sides of Eq. (29) by $`V_x^{k_1}V_y^{k_2}V_z^{k_3}`$ and then integrating over the velocity space. The resulting expression is

$`{\displaystyle \frac{}{x}}M_{k_1+1,k_2,k_3}+{\displaystyle \frac{}{y}}M_{k_1,k_2+1,k_3}`$ $`+`$ $`k_3{\displaystyle \frac{u_z}{x}}M_{k_1+1,k_2,k_31}+k_3{\displaystyle \frac{u_z}{y}}M_{k_1,k_2+1,k_31}`$ (33)
$``$ $`gk_3M_{k_1,k_2,k_31}=J_{k_1,k_2,k_3},`$
where

$$J_{k_1,k_2,k_3}=𝑑𝐯V_x^{k_1}V_y^{k_2}V_z^{k_3}J[f,f].$$
(34)
In the sequel, we will use the roman boldface $`𝐤`$ to denote the triad $`\{k_1,k_2,k_3\}`$ and the italic lightface $`k`$ to denote the sum $`k_1+k_2+k_3`$. Thus, $`M_𝐤M_{k_1,k_2,k_3}`$ is a moment of degree $`kk_1+k_2+k_3`$.

The collision integral $`J_𝐤`$ involves in general all the velocity moments of the distribution (including those of a degree higher than $`k`$) and its explicit expression in terms of those moments is unknown. An important exception is provided by the Maxwell interaction potential $`\phi (r)=K/r^4`$. In that case, the collision rate is independent of the velocity, i.e. $`|𝐯𝐯_1|\sigma (|𝐯𝐯_1|,\theta )=\sigma _0(\theta )`$, and thus $`J_𝐤`$ can be expressed as a bilinear combination of moments of degree equal to or smaller than $`k`$ :

$$J_𝐤=\underset{𝐤^{},𝐤^{\prime \prime }}{}^{}C_{𝐤^{},𝐤^{\prime \prime }}^𝐤M_𝐤^{}M_{𝐤^{\prime \prime }},$$
(35)
where the dagger denotes the constraint $`k^{}+k^{\prime \prime }=k`$. The coefficients $`C_{𝐤^{},𝐤^{\prime \prime }}^𝐤`$ are linear combinations of the eigenvalues

$$\lambda _r\mathrm{}=𝑑\mathrm{\Omega }\sigma _0(\theta )\left[1+\delta _{r0}\delta _\mathrm{}0\mathrm{cos}^{2r+\mathrm{}}\frac{\theta }{2}P_{\mathrm{}}\left(\mathrm{cos}\frac{\theta }{2}\right)\mathrm{sin}^{2r+\mathrm{}}\frac{\theta }{2}P_{\mathrm{}}\left(\mathrm{sin}\frac{\theta }{2}\right)\right]$$
(36)
of the linearized collision operator, where $`P_{\mathrm{}}(x)`$ are Legendre polynomials. The thermal conductivity and shear viscosity for Maxwell molecules (first obtained by Maxwell himself) are

$$\kappa (T)=\frac{5k_B}{2m}\frac{p}{n\lambda _{11}},$$
(37)
$$\eta (T)=\frac{p}{n\lambda _{02}},$$
(38)
where $`\lambda _{02}=\frac{3}{2}\lambda _{11}=0.436\times 3\pi \sqrt{2K/m}`$. The eigenvalues $`\lambda _r\mathrm{}`$ can also be expressed as linear combinations of the integrals

$$A_{2s}=𝑑\mathrm{\Omega }\sigma _0(\theta )\mathrm{sin}^{2s}\frac{\theta }{2}\mathrm{cos}^{2s}\frac{\theta }{2}.$$
(39)
In particular, $`\lambda _{02}=3A_2`$, $`\lambda _{04}=7A_2\frac{35}{4}A_4`$ and $`\lambda _{06}=\frac{27}{2}A_2\frac{357}{8}A_4+\frac{231}{8}A_6`$. The first few ratios $`A_{2s}/A_2`$ are $`A_4/A_2=0.15778\mathrm{}`$, $`A_6/A_2=0.031196\mathrm{}`$ and $`A_8/A_2=0.0066540\mathrm{}`$. The general expression (35) for $`J_𝐤`$ was elaborated by Truesdell and Muncaster , who also give $`J_𝐤`$ up to $`k=5`$. A simplification algorithm for this formula, along with explicit expressions of $`J_𝐤`$ for $`k=6`$$`8`$, has been recently proposed .

To simplify the subsequent analysis, it is suitable to use dimensionless quantities. First, let us define an effective collision frequency as

$$\nu =n\lambda _{02}.$$
(40)
Next, without loss of generality, we choose a reference frame stationary with the flow at $`r=0`$, i.e. $`u_0=0`$. Now we introduce the dimensionless variables

$`g^{}=v_0^1\nu _0^1g,u_z^{}=v_0^1u_z,p^{}=p_0^1p,T^{}=T_0^1T,`$
$`M_𝐤^{}=n_0^1v_0^kM_𝐤,J_𝐤^{}=n_0^1v_0^k\nu _0^1J_𝐤,𝐫^{}=\nu _0v_0^1𝐫.`$ (41)
In the above expressions, $`v_0=(k_BT_0/m)^{1/2}`$ is the thermal velocity and, as in Eqs. (13)–(18), the subscript $`0`$ means that the quantities are evaluated at $`r=0`$. The reduced (gravity) acceleration $`g^{}`$ has a clear physical meaning. The quantity $`h_0v_0^2/g`$ is the so-called scale height , i.e. the characteristic distance associated with the external (gravity) field. Thus, $`g^{}=(v_0/\nu _0)/h_0`$ represents the ratio between the mean free path and the distance through which a typical particle undergoes the action of $`𝐠`$. While the parameter $`g^{}`$ is a measure of the field strength on the scale of the mean free path, the Froude number

$$\text{Fr}=\frac{v_0^2}{gR}=\frac{h_0}{R}$$
(42)
measures the strength on the scale of the radius of the pipe. The relationship between both dimensionless parameters is

$$\text{Fr}=\frac{\text{Kn}}{g^{}},$$
(43)
where $`\text{Kn}(v_0/\nu _0)/R`$ is the Knudsen number.

The solution of the hierarchy of equations (33) is a very difficult task, even for Maxwell molecules, due to its nonlinear character and the fact that moments of degree $`k`$ are coupled to the first spatial derivatives of moments of degree $`k+1`$. On the other hand, we will restrict ourselves to cases where the parameter $`g^{}`$ is weak. This assumption is well justified in most of the practical situations, where values $`g^{}>10^2`$ would be abnormally large. For instance, in the case of air at room temperature under the action of terrestrial gravity, $`g^{}10^{12}`$. Consequently, we can solve the hierarchy (33) by means of a perturbation expansion in powers of $`g^{}`$, in a similar way as done in other works .

## 4 Expansion in powers of the external force

For the sake of clarity, the asterisks will be dropped in this Section, so all the quantities will be understood to be expressed in reduced units \[cf. Eqs. (41)\], unless stated otherwise.

Before performing the power series expansion of the moments, it is convenient to take into account the symmetry properties (32). For example, it is obvious that the profiles of the temperature and the hydrostatic pressure are even functions of $`g`$, while the flow velocity $`u_z`$ is an odd function. Therefore, we can write

$$p=1+p^{(2)}g^2+p^{(4)}g^4+\mathrm{},$$
(44)
$$T=1+T^{(2)}g^2+T^{(4)}g^4+\mathrm{},$$
(45)
$$u_z=u^{(1)}g+u^{(3)}g^3+\mathrm{}.$$
(46)
The expansion of the number density $`n=M_{000}=p/T`$ can be obtained from Eqs. (44) and (45). These hydrodynamic fields are (even) functions of $`x`$ and $`y`$ through the distance $`r=\sqrt{x^2+y^2}`$. As for the moments $`M_𝐤`$, they are even (odd) functions of $`g`$ if their index $`k_3`$ is even (odd). Thus,

$$M_𝐤=M_𝐤^{(0)}+M_𝐤^{(1)}g+M_𝐤^{(2)}g^2+M_𝐤^{(3)}g^3+\mathrm{},$$
(47)
where $`M_𝐤^{(s)}=0`$ if $`s+k_3=\text{odd}`$. The zeroth-order term $`M_𝐤^{(0)}`$ is the moment of the equilibrium distribution function normalized to $`p^{(0)}=1`$, $`T^{(0)}=1`$. Its expression is

$$M_𝐤^{(0)}=C_{k_1}C_{k_2}C_{k_3},$$
(48)
where $`C_k=(k1)!!=(k1)\times (k3)\times (k5)\times \mathrm{}\times 1`$ if $`k=\text{even}`$, being 0 if $`k=\text{odd}`$.

In the first stage of the calculations, the expansions (44)–(47) are inserted into Eq. (33). By equating the terms of degree $`s`$ in $`g`$ on both sides, we get

$`{\displaystyle \frac{}{x}}M_{k_1+1,k_2,k_3}^{(s)}`$ $`+`$ $`{\displaystyle \frac{}{y}}M_{k_1,k_2+1,k_3}^{(s)}+k_3{\displaystyle \underset{i=0}{\overset{[(s1)/2]}{}}}[{\displaystyle \frac{u^{(2i+1)}}{x}}M_{k_1+1,k_2,k_31}^{(s12i)}`$ (49)
$`+`$ $`{\displaystyle \frac{u^{(2i+1)}}{y}}M_{k_1,k_2+1,k_31}^{(s12i)}]k_3M^{(s1)}_{k_1,k_2,k_31}=J^{(s)}_𝐤,`$
where $`[\alpha ]`$ means the integer part of $`\alpha `$. According to Eq. (35), the right side of (49) is

$$J_𝐤^{(s)}=\underset{𝐤^{},𝐤^{\prime \prime }}{}^{}C_{𝐤^{},𝐤^{\prime \prime }}^𝐤\underset{i=0}{\overset{s}{}}M_𝐤^{}^{(si)}M_{𝐤^{\prime \prime }}^{(i)}.$$
(50)
Equation (50) is still very difficult to solve in general. First, it is nonlinear and couples moments of degree $`k`$ and lower to moments of degree $`k+1`$. Second, the coefficients $`M_𝐤^{(s)}`$ are unknown functions of the spatial variables $`x`$ and $`y`$. To overcome the latter obstacle, we assume that a self-consistent solution of Eq. (49) exists in which the coefficients in the series defined by (44)–(47) are polynomials in $`x`$ and $`y`$. This ansatz is justified by the results obtained in the case of the planar geometry from the BGK model  and the Boltzmann equation , as well as in the case of the cylindrical geometry from the BGK model . More specifically, we assume that the coefficients of the hydrodynamic profiles are of the form

$$p^{(s)}=\underset{i=1}{\overset{s1}{}}p_{2i}^{(s)}r^{2i},$$
(51)
$$T^{(s)}=\underset{i=1}{\overset{s}{}}T_{2i}^{(s)}r^{2i},$$
(52)
$$u^{(s)}=\underset{i=1}{\overset{s}{}}u_{2i}^{(s)}r^{2i}.$$
(53)
The task becomes more difficult when guessing the spatial dependence of the coefficients $`M_𝐤^{(s)}`$. It can be checked that the only representation that leads to consistent solutions of the hierarchy (49) is of the form

$$M_𝐤^{(s)}=\underset{i=0}{\overset{\sigma (s)}{}}\underset{j=0}{\overset{i}{}}\chi _𝐤^{(s,ij,j)}x^{ij}y^j,$$
(54)
where $`\sigma (s=1)=1`$ and $`\sigma (s)=2s`$ for $`s>1`$. From Eqs. (32) it follows that $`\chi _{k_1,k_2,k_3}^{(s,i,j)}=0`$ if $`s+k_3=\text{odd}`$ or $`i+k_1=\text{odd}`$ or $`j+k_2=\text{odd}`$. Insertion of Eq. (54) into Eq. (50) shows that $`J_𝐤^{(s)}`$ has the same structure as $`M_𝐤^{(s)}`$, namely

$$J_𝐤^{(s)}=\underset{i=0}{\overset{\sigma (s)}{}}\underset{j=0}{\overset{i}{}}J_𝐤^{(s,ij,j)}x^{ij}y^j.$$
(55)
The numerical coefficients $`p_{2i}^{(s)}`$, $`T_{2i}^{(s)}`$, $`u_{2i}^{(s)}`$ and $`\chi _𝐤^{(s,i,j)}`$ are determined by consistency. To do so, we insert Eqs. (51)–(54) into Eq. (49) and equate the coefficients of the terms $`x^iy^j`$ in both sides. This yields a hierarchy of linear equations for the unknowns. In the hierarchy, the coefficients of the moments of degree $`k+1`$, $`k`$ and $`k1`$ on the left-hand side are related to those of degree $`k`$, $`k1`$, $`k2`$, …on the right-hand side. In addition, the coefficients of the type $`\chi _{k_1+1,k_2,k_3}^{(s,i+1,j)}`$ and $`\chi _{k_1,k_2+1,k_3}^{(s,i,j+1)}`$ on the left-hand side determine those of the type $`\chi _{k_1,k_2,k_3}^{(s,i,j)}`$ on the right-hand side. In spite of the intricacy of the resulting hierarchy, it can be recursively solved by following steps similar to those made in the case of the planar geometry . To get explicit results to order $`s=3`$ it is necessary to consider the moment equations up to degree $`k=8`$.

The method is outlined in Appendix B and the main results are presented in the next Section.

## 5 Discussion

### 5.1 Hydrodynamic profiles

By proceeding along the recursive scheme described in Appendix B, we have obtained the explicit spatial dependence of the coefficients in the series expansions (44)–(46) through order $`g^4`$ and those of the series expansion (47) for moments of second ($`k=2`$) and third ($`k=3`$) degrees through order $`g^3`$. Here we display the results to third order. Also, to ease the physical interpretation of the results, they are expressed in cylindrical coordinates and in real units. The hydrodynamic profiles are

$$u_z(r)=u_0\frac{\rho _0r^2g}{4\eta _0}\left[1+\frac{1}{144}\frac{\rho _0^2g^2}{\kappa _0\eta _0T_0}r^4+\zeta _u\left(\frac{mg}{k_BT_0}\right)^2r^2+\zeta _u^{}\frac{\rho _0\eta _0^2g^2}{p_0^3}\right]+𝒪(g^5),$$
(56)
$$p(r)=p_0\left[1+\zeta _p\left(\frac{mg}{k_BT_0}\right)^2r^2\right]+𝒪(g^4),$$
(57)
$$T(r)=T_0\left[1\frac{1}{64}\frac{\rho _0^2g^2}{\kappa _0\eta _0T_0}r^4+\zeta _T\left(\frac{mg}{k_BT_0}\right)^2r^2\right]+𝒪(g^4),$$
(58)
where $`\zeta _u=\frac{4429}{12600}0.35`$, $`\zeta _u^{}=24.322\mathrm{}`$, $`\zeta _p=\frac{3}{10}`$ and $`\zeta _T=\frac{34}{175}0.19`$. We recall that the subscript $`0`$ in a quantity represents its value on the $`z`$-axis, i.e. at $`r=0`$. Elimination of $`r`$ among Eqs. (56)–(58) allows one to get the following nonequilibrium equation of state:

$$T=T_0\frac{\eta _0}{4\kappa _0}(u_zu_0)^2+\frac{\zeta _T}{\zeta _p}\frac{T_0}{p_0}(pp_0)+𝒪(g^4).$$
(59)
The non-zero elements of the pressure tensor are

$$P_{rr}(r)=p_0\left[1+\frac{1}{6}\zeta _p\left(\frac{mg}{k_BT_0}\right)^2r^2\zeta _P^{}\frac{\rho _0\eta _0^2g^2}{p_0^3}\right]+𝒪(g^4),$$
(60)
$$P_{\varphi \varphi }(r)=p_0\left[1+\frac{1}{2}\zeta _p\left(\frac{mg}{k_BT_0}\right)^2r^2\zeta _P^{}\frac{\rho _0\eta _0^2g^2}{p_0^3}\right]+𝒪(g^4),$$
(61)
$$P_{zz}(r)=p_0\left[1+\frac{7}{3}\zeta _p\left(\frac{mg}{k_BT_0}\right)^2r^2+2\zeta _P^{}\frac{\rho _0\eta _0^2g^2}{p_0^3}\right]+𝒪(g^4),$$
(62)
$$P_{rz}(r)=\frac{\rho _0g}{2}r\left[1+\frac{1}{192}\frac{\rho _0^2g^2}{\kappa _0\eta _0T_0}r^4+\frac{\zeta _p\zeta _T}{2}\left(\frac{mg}{k_BT_0}\right)^2r^2\right]+𝒪(g^5),$$
(63)
where $`\zeta _P^{}=1.7388\mathrm{}`$. Finally, the components of the heat flux are

$$q_r(r)=\frac{1}{16}\frac{\rho _0^2g^2}{\eta _0}r^3+𝒪(g^4),$$
(64)
$$q_z(r)=\frac{2}{5}\frac{\kappa _0mg}{k_B}\left[\zeta _0\zeta _q\frac{\rho _0^2g^2}{\kappa _0\eta _0T_0}r^4\zeta _q^{}\left(\frac{mg}{k_BT_0}\right)^2r^2\zeta _q^{\prime \prime }\frac{\rho _0\eta _0^2g^2}{p_0^3}\right]+𝒪(g^5),$$
(65)
where $`\zeta _0=1`$, $`\zeta _q=\frac{37}{28}0.29`$, $`\zeta _q^{}=3.2527\mathrm{}`$ and $`\zeta _q^{\prime \prime }=86.415\mathrm{}`$.

Equations (56)–(65) are the main results of this paper. They give the exact profiles in the bulk region of the hydrodynamic fields and their fluxes for a dilute gas of Maxwell molecules through third order in the external force. Comparison of Eqs. (6), (9) and (13)–(18) with Eqs. (56)–(65) show that the NS constitutive equations predict $`\zeta _u=\zeta _u^{}=\zeta _p=\zeta _T=\zeta _P^{}=\zeta _0=\zeta _q=\zeta _q^{}=\zeta _q^{\prime \prime }=0`$. These vanishing values of the coefficients summarize the main limitations of the NS description of the Poiseuille flow induced by an external force. Equation (65) shows that $`q_z0`$, even to first order in $`g`$ ($`\zeta _00`$), despite the absence of a thermal gradient along the longitudinal direction. This represents an obvious breakdown of Fourier’s law and is in fact a Burnett-order effect . To second order in $`g`$ the hydrostatic pressure is not uniform but increases radially ($`\zeta _p>0`$); analogously, there are normal stress differences, so that $`P_{rr}<P_{\varphi \varphi }<p<P_{zz}`$. These are non-Newtonian effects.

Also to second order in $`g`$, there exists a positive quadratic term ($`\zeta _T>0`$) in addition to the negative quartic term in the temperature profile. The former term is responsible for the fact that the temperature has a local minimum rather than a maximum at $`r=0`$. The maximum temperature is located at a distance $`r_{\text{max}}=\sqrt{32\zeta _T}\mathrm{}_02.49\mathrm{}_0`$ (independent of $`g`$), where $`\mathrm{}_0(\eta _0\kappa _0T_0)^{1/2}/p_0=\sqrt{\frac{2}{5}\text{Pr}}(v_0/\nu _0)`$ is the (local) mean free path at $`r=0`$; the relative difference between the maximum temperature and the temperature on the cylinder axis is $`(T_{\text{max}}T_0)/T_0=16\zeta _T^2(\mathrm{}_0/h_0)^20.604(\mathrm{}_0/h_0)^2`$, where we recall that $`h_0k_BT_0/mg`$. This result represents a dramatic violation of Fourier’s law: while the heat flows radially outwards, the temperature increases from $`r=0`$ to $`r=r_{\text{max}}`$. Therefore, within the region $`0rr_{\text{max}}`$, the heat flows from the colder to the hotter layers. This paradoxical effect, which is beyond the Burnett description , does not violate the conservation of energy (5) since the radial increase of $`rq_r`$ (which tends to cool the gas) is exactly compensated for by the viscous heating term $`rP_{rz}u_z/r`$. It is interesting to note that from Eqs. (58) and (64) one can obtain

$$\kappa \frac{T}{r}=q_r\frac{1}{9}r_{\text{max}}^2^2q_r+𝒪(g^4),$$
(66)
where $`^2X=r^1(rX/r)/r`$ is the Laplacian in cylindrical coordinates. Equation (66) is an extension of Fourier’s law (8) showing that the sign of the thermal gradient results from a competition between the radial component of the heat flux and its Laplacian . The latter dominates for $`r<r_{\text{max}}`$ and so $`T/r>0`$ in that region.

Third-order contributions appear in $`u_z`$ and $`P_{rz}`$ but they do not have a qualitatively important influence. They are responsible for a spatial variation of the flow velocity $`u_z`$ and the shear stress $`P_{rz}`$ more pronounced than expected from the NS equations. The same happens for the mass rate of flow defined by Eq. (19). Insertion of the hydrodynamic profiles (56)–(58) yields

$`\dot{M}(a)`$ $`=`$ $`{\displaystyle \frac{\pi \rho _0^2ga^4}{8\eta _0}}\left[1+{\displaystyle \frac{5}{384}}{\displaystyle \frac{\rho _0^2g^2}{\eta _0\kappa _0T_0}}a^4+{\displaystyle \frac{4\zeta _u+\zeta _p\zeta _T}{3}}\left({\displaystyle \frac{mg}{k_BT_0}}\right)^2a^2+\zeta _u^{}{\displaystyle \frac{\rho _0\eta _0^2g^2}{p_0^3}}\right]`$ (67)
$`+𝒪(g^5).`$
Comparison with Eq. (20) shows that the mass rate of flow grows with the radius $`a`$ more rapidly than predicted by the NS approximation, but otherwise the qualitative behavior is similar in both descriptions.

### 5.2 Generalized constitutive equations

Let us use Eqs. (56)–(65) to rewrite the fluxes under the form of generalized constitutive equations. We begin with the shear stress:

$`P_{rz}`$ $`=`$ $`\eta {\displaystyle \frac{u_z}{r}}[1{\displaystyle \frac{\zeta _u^{}}{15\zeta _T}}{\displaystyle \frac{\mathrm{}^2}{T}}^2T{\displaystyle \frac{4}{15}}({\displaystyle \frac{\zeta _u^{}}{15\zeta _T}}+8\zeta _u+6\zeta _T2\zeta _p)`$ (68)
$`\times {\displaystyle \frac{\mathrm{}^2}{k_BT/m}}\left({\displaystyle \frac{u_z}{r}}\right)^2]+𝒪(g^5).`$
Here $`\eta T`$ is the local shear viscosity at $`r`$, so $`\eta \eta _0TT_0=𝒪(g^2)`$. Analogously, $`\mathrm{}=\left(\eta \kappa T\right)^{1/2}/p`$ is the local mean free path and $`\mathrm{}\mathrm{}_0=𝒪(g^2)`$. In Eq. (68) the NS term is of order $`g`$, while the super-Burnett terms (i.e. of third order in the hydrodynamic gradients) are of order $`g^3`$. The normal stresses can be written as

$$P_{rr}=p\left[1+\frac{10\zeta _p}{3}\frac{\mathrm{}^2}{T}^2T\left(\frac{\zeta _P^{}}{15\zeta _p}+\frac{10\zeta _T}{3}\right)\frac{\mathrm{}^2}{p}^2p\right]+𝒪(g^4),$$
(69)
$$P_{\varphi \varphi }=p\left[1+2\zeta _p\frac{\mathrm{}^2}{T}^2T\left(\frac{\zeta _P^{}}{15\zeta _p}+2\zeta _T\right)\frac{\mathrm{}^2}{p}^2p\right]+𝒪(g^4),$$
(70)
$$P_{zz}=p\left[1\frac{16\zeta _p}{3}\frac{\mathrm{}^2}{T}^2T+\left(\frac{2\zeta _P^{}}{15\zeta _p}+\frac{16\zeta _T}{3}\right)\frac{\mathrm{}^2}{p}^2p\right]+𝒪(g^4).$$
(71)
Thus, in order to obtain the normal stresses to order $`g^2`$ we need to include at least the Burnett contributions in the Chapman–Enskog expansion.

The radial component of the heat flux can be expressed as

$$q_r=\kappa \frac{}{r}\left(T+4\zeta _T\mathrm{}^2^2T\right)+𝒪(g^4).$$
(72)
In this case, both the NS term and the super-Burnett term are of order $`g^2`$. As for the longitudinal heat flux, one has

$`q_z`$ $`=`$ $`{\displaystyle \frac{\eta ^2}{2\rho }}\left(^2u_z\right)\left(\theta _4{\displaystyle \frac{\zeta _u^{}+\zeta _q^{\prime \prime }}{5\zeta _T}}^2T\right)+{\displaystyle \frac{2}{5}}nk_B\mathrm{}^2{\displaystyle \frac{u_z}{r}}`$ (73)
$`\times \left(\theta _5{\displaystyle \frac{T}{r}}+{\displaystyle \frac{\theta _3}{3}}{\displaystyle \frac{T}{p}}{\displaystyle \frac{p}{r}}+\mu \mathrm{}^2{\displaystyle \frac{}{r}}^2T\right)+𝒪(g^5),`$
where $`\theta _3=3`$, $`\theta _4=3\zeta _0=3`$, $`\theta _5=\frac{1}{2}(1+64\zeta _q)=\frac{39}{4}`$ and

$$\mu =\frac{\zeta _u^{}+\zeta _q^{\prime \prime }}{15\zeta _T}+16\zeta _u+2\zeta _T(7+64\zeta _q)+4\zeta _q^{}8\zeta _p64.141\mathrm{}.$$
(74)
The Burnett term headed by $`\theta _4`$ is of of first order in $`g`$, while the Burnett terms with $`\theta _3`$ and $`\theta _5`$ are of third order. In addition, there are super-super-Burnett terms that are also of order $`g^3`$.

Equations (68)–(73) can be written in other equivalent forms by taking into account that the hydrodynamic gradients are not independent in our problem. For instance,

$$\frac{m}{k_BT}\left(\frac{u_z}{r}\right)^2=\frac{15}{4}\left(\frac{\zeta _T}{\zeta _p}\frac{1}{p}^2p\frac{1}{T}^2T\right)+𝒪(g^4),$$
(75)
$$^2p=4\zeta _p\frac{p}{T}\mathrm{}^2^4T+𝒪(g^4)=\frac{16}{15}\zeta _p\frac{mp}{k_BT}\mathrm{}^2\left(^2u_z\right)^2+𝒪(g^4).$$
(76)
In particular, Eq. (76) implies that the Burnett term $`^2p`$ and the super-super-Burnett terms $`^4T`$ and $`\left(^2u_z\right)^2`$ are of order $`g^2`$. Therefore, in order to obtain correctly the profiles (56)–(65) from the Chapman–Enskog method one would actually need to go beyond the apparent orders in the gradients of Eqs. (68)–(73).

The previous analysis shows that the Chapman–Enskog expansion and the expansion in powers of $`g`$ in the Poiseuille flow are quite different, even though the gradients are induced by the external force. Both expansions require a weak field, namely $`h_0\mathrm{}_0`$. On the other hand, while the Chapman–Enskog expansion requires that $`\mathrm{}_0^2T_0^1^2T1`$, i.e. $`r\mathrm{}_0`$, the expansion in powers of $`g`$ remains valid for distances smaller than or of the order of the mean free path. Since the non-monotonic behavior of the temperature profile occurs for $`0rr_{\text{max}}\mathrm{}_0`$, this effect is neglected by the NS description.

### 5.3 Comparison with the BGK solution

Once we have discussed the limits of validity of the NS theory (which otherwise provides a powerful tool to study the macroscopic state of many nonequilibrium flows) in our Poiseuille problem, it is worthwhile wondering whether the BGK model kinetic equation succeeds in capturing the main qualitative features of the Boltzmann solution. The solution of the BGK model for the cylindrical Poiseuille flow was obtained by Tij and Santos  to fourth order in $`g`$ for a general class of power-law repulsive potentials. The results agree exactly with the structure of Eqs. (56)–(65). On the other hand, most of the coefficients $`\zeta `$’s have different numerical values. The failure of the BGK model to reproduce the exact Boltzmann values for those coefficients is not surprising, given the simplicity of the model, which assumes that the practical effect of collisions is to make the distribution function relax to local equilibrium with a single characteristic rate $`\nu `$. In fact, it is well known  that the BGK model yields a value ($`\text{Pr}=1`$) for the Prandtl number $`\text{Pr}5k_B\eta /2m\kappa `$ different from the correct one ($`\text{Pr}=\frac{2}{3}`$). Table 1 compares the numerical values of the coefficients $`\zeta `$’s given by the Boltzmann equation for Maxwell molecules (BM) with those obtained from the BGK model. Except for $`\zeta _p`$ and $`\zeta _0`$, the coefficients are different in both theories, the ratio being bounded between $`0.38`$ and $`1.38`$. For instance, the location $`r_{\text{max}}`$ of the maximum temperature is in the BGK model 15% smaller than in the Boltzmann equation and the temperature difference $`T_{\text{max}}T_0`$ in the BGK model is almost half the Boltzmann value. A comparison of the coefficients appearing in the hydrodynamic fields to fourth order in $`g`$ \[not shown in Eqs. (57) and (58), but given in Appendix B\] shows that the discrepancies between the BGK and Boltzmann equations increase with the power of $`g`$. In particular, the series obtained from the BGK model seem to diverge more quickly than those obtained from the Boltzmann equation.

In order to illustrate graphically the differences between the NS, BGK and BM descriptions, we plot in Figs. 16 the hydrodynamic and flux profiles for the case $`g=0.1v_0^2/\mathrm{}_0`$ \[which corresponds to $`h_0/\mathrm{}_0=10`$\], when only terms through third order in $`g`$ are retained. Of course, higher order terms are not expected to be negligible for that particular high value of the external field , especially as one departs from the axis, as discussed in Subsection 5.5. However, the retained terms are sufficient to show the qualitative differences in the predictions of the three approaches, so that the main aim of Figs. 16 is to highlight those differences. It must be born in mind that the BGK curves in Figs. 16 correspond to Eqs. (56)–(65) with the BGK values for the coefficients $`\zeta `$’s (cf. Table 1) but with the correct value $`\text{Pr}=\frac{2}{3}`$ of the Prandtl number. As expected, the BGK model is in qualitative agreement with the BM results.

### 5.4 Comparison with the planar Poiseuille flow

Let us now add a few comments on the similarities between the cylindrical and the planar Poiseuille flows. The solution of the Boltzmann equation for Maxwell molecules to order $`g^2`$ in the planar geometry  yields profiles analogous to those of Eqs. (56)–(65). In particular, the temperature profile is

$$T(y)=T_0\left[1\frac{1}{12}\frac{\rho _0^2g^2}{\kappa _0\eta _0T_0}y^4+\zeta _T\left(\frac{mg}{k_BT_0}\right)^2y^2\right]+𝒪(g^4),$$
(77)
where now $`\zeta _T=1.0153`$. As a consequence, the temperature has a maximum with $`(T_{\text{max}}T_0)/T_0=3\zeta _T^2(\mathrm{}_0/h_0)^23.09(\mathrm{}_0/h_0)^2`$ at $`|y|=y_{\text{max}}=\sqrt{6\zeta _T}\mathrm{}_02.47\mathrm{}_0`$. Therefore, the maximum temperature occurs at a separation $`y_{\text{max}}`$ from the middle plane practically equal to the radial distance $`r_{\text{max}}`$ of the cylindrical case, but the relative temperature change $`(T_{\text{max}}T_0)/T_0`$ is about five times larger in the former case than in the latter. Interestingly, this effect is accurately captured by the BGK model , as explained by the fact that in the planar Poiseuille flow one has  $`\zeta _p^{\text{BM}}/\zeta _p^{\text{BGK}}=1`$, $`\zeta _T^{\text{BM}}/\zeta _T^{\text{BGK}}=1.0153/(19/25)1.34`$, $`\zeta _{P}^{}{}_{}{}^{\text{BM}}/\zeta _{P}^{}{}_{}{}^{\text{BGK}}=6.2602/(306/25)0.51`$, $`\zeta _0^{\text{BM}}/\zeta _0^{\text{BBGK}}=1`$, in close agreement with the corresponding ratios in the cylindrical Poiseuille flow (see Table 1).

### 5.5 Boundary conditions

When solving the Boltzmann hierarchy of moments \[cf. Eq. (33)\] by a perturbation expansion in powers of gravity we have taken the hydrodynamic quantities at the axis ($`r=0`$) as reference values. In so doing, we have avoided the need of imposing specific boundary conditions. The price to be paid is that it is difficult to determine how far from the axis the profiles (56)–(65) remain valid. Consider for instance the temperature profile (58). The first neglected term $`𝒪(g^4)`$ is a polynomial in $`r`$ of degree 8 \[cf. Eq. (109)\]. This term can actually be neglected versus the retained terms only if $`g_{}^{}{}_{}{}^{4}(r/\mathrm{}_0)^8g_{}^{}{}_{}{}^{2}(r/\mathrm{}_0)^4`$, i.e. if

$$g^{}(r/\mathrm{}_0)^21.$$
(78)
The same conclusion is obtained from the velocity profile (56) and the pressure profile (57). If $`r\mathrm{}_0`$, condition (78) is equivalent to $`g^{}1`$. On the other hand, a much stronger condition $`g^{}\text{Kn}^2`$ is needed to extend the profiles (56)–(65) to distances $`r`$ comparable with the pipe radius $`R`$. In order to find the hydrodynamic fields for $`rR`$ when condition (78) does not hold, one would need to solve the full Boltzmann equation with the appropriate (e.g. diffuse reflection) boundary conditions at $`r=R`$ corresponding to a given wall temperature $`T_w`$. This has been done by Aoki et al.  in the case of the BGK equation for the planar geometry. They derived the hydrodynamic equations and the matching conditions corresponding to the normal (or Hilbert) solution by an asymptotic expansion in powers of the Knudsen number Kn and for values of the Froude number of order $`\text{Fr}1/\text{Kn}`$. In view of Eq. (43), this implies that $`g^{}\text{Kn}^2`$. The successive hydrodynamic equations had to be solved numerically, but Aoki et al. verified that a Taylor series expansion around the center of the gap agreed with the perturbation solution found in Ref. .

An asymptotic method analogous to that worked out by Aoki et al.  would be much more difficult to carry out in the case of the Boltzmann equation for the cylindrical geometry. Nevertheless, it is natural to expect that the results derived in this paper would be recovered by performing a Taylor series expansion around $`r=0`$ of the solution for $`g^{}\text{Kn}^2`$. This in addition would provide the matching conditions for the profiles (56)–(58) in the region $`r/\mathrm{}_0g_{}^{}{}_{}{}^{1/2}\text{Kn}^1`$.

## 6 Conclusion

In this work, we have studied the (laminar) stationary Poiseuille flow in a cylindrical pipe induced by an external force $`m𝐠`$. We have considered a dilute gas of Maxwell molecules and have analyzed the hierarchy of moment equations associated with the nonlinear Boltzmann equation. A consistent solution has been found under the form of a perturbation expansion in powers of $`g`$ through third order. In principle, the method can be pushed to higher orders but, not only the algebraic intricacy of the solution scheme grows rapidly, but the results would not be of practical use given the suspected asymptotic character of the series. The discussion of the results has focused on the most relevant physical quantities, namely the hydrodynamic profiles (flow velocity, pressure and temperature) and the associated fluxes (pressure tensor and heat flux). The results obtained here show once again that the Navier–Stokes (NS) theory remains unable to envisage the correct hydrodynamic profiles, even at a qualitative level, for distances smaller than or of the order of the mean free path. The most important limitation of the NS description is that it predicts a monotonically decreasing temperature as one moves apart from the cylinder axis. In contrast, our solution to the Boltzmann equation shows that the temperature has a local minimum ($`T=T_0`$) at the axis ($`r=0`$) and reaches a maximum value $`T=T_{\text{max}}T_0\left[1+0.6(\mathrm{}_0/h_0)^2\right]`$ at a distance from the center $`r=r_{\text{max}}2.5\mathrm{}_0`$ of the order of the mean free path $`\mathrm{}_0`$. In addition, a longitudinal component of the heat flux exists ($`q_z0`$) in the absence of gradients along the longitudinal direction and normal stress differences ($`P_{rr}<P_{\varphi \varphi }<p<P_{zz}`$) are present. On the other hand, the BGK model provides results qualitatively similar to those found in this paper . Notwithstanding this, the profiles predicted by the BGK model require some corrections from a quantitative point of view, as expected. In particular, the BGK model underestimates the maximum value $`T_{\text{max}}`$ of the temperature, as well as its location $`r_{\text{max}}`$. From that point of view, one might say that the BGK model deviates from the NS predictions less than the Boltzmann equation, at least for Maxwell molecules. The same conclusion arises from a comparison between the Boltzmann and BGK results for the planar Poiseuille flow . Nevertheless, it is worth noting that the BGK model is accurate in describing the effect of the geometry of the Poiseuille flow on the profiles. For instance, it correctly predicts that in the planar case the maximum temperature is located at a separation $`y_{\text{max}}`$ from the middle plane practically equal to the radial distance $`r_{\text{max}}`$ of the cylindrical case, although the relative temperature change $`(T_{\text{max}}T_0)/T_0`$ is about five times larger in the former case than in the latter.

Although we have not considered interactions different from Maxwell’s, we expect the results found here to remain essentially valid, except that the numerical values of the coefficients may change. This expectation is supported by the results obtained with the BGK kinetic model . In that case one finds that the coefficients $`\zeta _p`$, $`\zeta _T`$, $`\zeta _P^{}`$ and $`\zeta _0`$ are universal, while the coefficients $`\zeta _u`$, $`\zeta _u^{}`$, $`\zeta _q`$, $`\zeta _q^{}`$ and $`\zeta _q^{\prime \prime }`$ change by 5%, 0.5%, 3%, 0.7% and 0.8%, respectively, when going from Maxwell molecules to hard spheres.

Let us conclude by asserting that the stationary Poiseuille flow driven by an external force, both in a slab and in a pipe, is a good and conceptually simple example in fluid dynamics showing the limitations of a purely continuum theory in contrast to a kinetic theory approach. Next, we want to emphasize the important role played by the BGK model kinetic equation to pave the way towards a more fundamental theory based on the Boltzmann equation.

A.S. gratefully acknowledges partial support from the Ministerio de Ciencia y Tecnología (Spain) through grant No. BFM2001-0718.

## Appendix A Cylindrical coordinates

The relationships between the cylindrical and Cartesian components of a vector $`𝐀`$ and a tensor $`𝖡`$ are

$$\left(\begin{array}{c}A_r\\ A_\varphi \\ A_z\end{array}\right)=𝖴\left(\begin{array}{c}A_x\\ A_y\\ A_z\end{array}\right),$$
(79)
$$\left(\begin{array}{ccc}B_{rr}& B_{r\varphi }& B_{rz}\\ B_{\varphi r}& B_{\varphi \varphi }& B_{\varphi z}\\ B_{zr}& B_{z\varphi }& B_{zz}\end{array}\right)=𝖴\left(\begin{array}{ccc}B_{xx}& B_{xy}& B_{xz}\\ B_{yx}& B_{yy}& B_{yz}\\ B_{zx}& B_{zy}& B_{zz}\end{array}\right)𝖴^{},$$
(80)
where

$$𝖴=\left(\begin{array}{ccc}x/r& y/r& 0\\ y/r& x/r& 0\\ 0& 0& 1\end{array}\right)$$
(81)
is a unitary matrix and $`𝖴^{}`$ is its transpose.

In general, the divergence of a vector $`𝐀`$ in cylindrical coordinates is

$$𝐀=\frac{1}{r}\frac{}{r}\left(rA_r\right)+\frac{1}{r}\frac{}{\varphi }A_\varphi +\frac{}{z}A_z.$$
(82)
The cylindrical components of the divergence of a tensor $`𝖡`$ are

$$\left(𝖡\right)_r=\frac{1}{r}\frac{}{r}\left(rB_{rr}\right)+\frac{1}{r}\frac{}{\varphi }B_{\varphi r}+\frac{}{z}B_{zr}\frac{1}{r}B_{\varphi \varphi },$$
(83)
$$\left(𝖡\right)_\varphi =\frac{1}{r}\frac{}{r}\left(rB_{r\varphi }\right)+\frac{1}{r}\frac{}{\varphi }B_{\varphi \varphi }+\frac{}{z}B_{z\varphi }+\frac{1}{r}B_{\varphi r},$$
(84)
$$\left(𝖡\right)_z=\frac{1}{r}\frac{}{r}\left(rB_{rz}\right)+\frac{1}{r}\frac{}{\varphi }B_{\varphi z}+\frac{}{z}B_{zz}.$$
(85)
In the Poiseuille problem the cylindrical components of the flow velocity $`𝐮`$, the heat flux $`𝐪`$ and the pressure tensor $`𝖯`$ depend on the radial variable $`r`$ only. Consequently,

$$𝗊=\frac{1}{r}\frac{}{r}\left(rq_r\right),$$
(86)
$$\left(𝖯\right)_r=\frac{1}{r}\frac{}{r}\left(rP_{rr}\right)\frac{1}{r}P_{\varphi \varphi },$$
(87)
$$\left(𝖯\right)_\varphi =0,$$
(88)
$$\left(𝖯\right)_z=\frac{1}{r}\frac{}{r}\left(rP_{rz}\right),$$
(89)
$$𝖯:𝐮=P_{rz}\frac{u_z}{r}.$$
(90)
In Eq. (88) we have taken into account that, by symmetry, $`P_{r\varphi }=P_{\varphi r}=0`$.

## Appendix B Solution of the hierarchy of moment equations

In this Appendix we outline the steps needed to solve the hierarchy (49) with the profiles (51)–(54). To do so, we proceed sequentially from a given order $`s`$ to the next order $`s+1`$. It is also necessary to use the consistency relation

$$p=\frac{1}{3}\left(M_{2,0,0}+M_{0,2,0}+M_{0,0,2}\right)$$
(91)
to any order in $`g`$. Let us start from the first-order equations.

### B.1 First order in $`g`$ $`(s=1)`$

Setting $`s=1`$, Eq. (49) becomes

$`{\displaystyle \frac{}{x}}M_{k_1+1,k_2,k_3}^{(1)}`$ $`+`$ $`{\displaystyle \frac{}{y}}M_{k_1,k_2+1,k_3}^{(1)}+k_3{\displaystyle \frac{u^{(1)}}{x}}M_{k_1+1,k_2,k_31}^{(0)}`$ (92)
$`+`$ $`k_3{\displaystyle \frac{u^{(1)}}{y}}M_{k_1,k_2+1,k_31}^{(0)}k_3M_{k_1,k_2,k_31}^{(0)}=J_𝐤^{(1)}.`$
Next, inserting (53) and (54) into (92), we obtain

$`\chi _{k_1+1,k_2,k_3}^{(1,1,0)}`$ $`+`$ $`\chi _{k_1,k_2+1,k_3}^{(1,0,1)}+2k_3M_{k_1+1,k_2,k_31}^{(0)}u_2^{(1)}x+2k_3M_{k_1,k_2+1,k_31}^{(0)}u_2^{(1)}y`$
$``$ $`k_3M_{k_1,k_2,k_31}^{(0)}=J_𝐤^{(1,0,0)}+J_𝐤^{(1,1,0)}x+J_𝐤^{(1,0,1)}y,`$
where $`J_𝐤^{(s,i,j)}`$ is defined by Eq. (55). Equating the coefficients of the same degree in $`x`$ and $`y`$, Eq. (LABEL:51) decouples into the following set of equations:

$$2k_3M_{k_1+1,k_2,k_31}^{(0)}u_2^{(1)}=J_{k_1,k_2,k_3}^{(1,1,0)},$$
(94)
$$2k_3M_{k_1,k_2+1,k_31}^{(0)}u_2^{(1)}=J_{k_1,k_2,k_3}^{(1,0,1)},$$
(95)
$$\chi _{k_1+1,k_2,k_3}^{(1,1,0)}+\chi _{k_1,k_2+1,k_3}^{(1,0,1)}k_3M_{k_1,k_2,k_31}^{(0)}=J_{k_1,k_2,k_3}^{(1,0,0)}.$$
(96)
Let us start with the moments of small degree ($`k=2`$). For example, by taking $`\{k_1,k_2,k_3\}=\{1,0,1\}`$ in (94) and $`\{k_1,k_2,k_3\}=\{0,1,1\}`$ in (95,) we obtain

$$2u_2^{(1)}=\chi _{1,0,1}^{(1,1,0)},$$
(97)
$$2u_2^{(1)}=\chi _{0,1,1}^{(1,0,1)}.$$
(98)
Now we insert the above equations into Eq. (96) with $`\{k_1,k_2,k_3\}=\{0,0,1\}`$. This allows us to obtain the coefficient of the flow velocity profile to first order in $`g`$:

$$u_2^{(1)}=\frac{1}{4}.$$
(99)
The process of solution continues in the same manner for the coefficients of moments of degree $`k>2`$. We use the relations (94) and (95) in a first stage to determine the coefficients $`\chi _𝐤^{(1,1,0)}`$ and $`\chi _𝐤^{(1,0,1)}`$ with $`k=4`$. In the second stage, the obtained results are inserted into (96) to evaluate the coefficients of degree $`k=3`$. The same steps can continue infinitely by following the routes $`k=6k=5`$; $`k=8k=7`$; …. It is not necessary to evaluate all the coefficients, but some of them are essential to the solution of the problem to the next order in $`g`$.

### B.2 Second order in $`g`$ $`(s=2)`$

As done in the first-order evaluation, let us take $`s=2`$ in Eq. (49) and insert Eqs. (51)–(54). We then obtain an equation more complicated than that of the first order. Of course, the symmetry property $`\chi _𝐤^{(s,i,j)}=0`$ if $`i+k_1=\text{odd}`$ or $`j+k_2=\text{odd}`$ must be applied. By equating the coefficients of the same degree in $`x`$ and $`y`$ we get the following set of algebraic equations:

1. $`i+j=4`$

$$0=J_𝐤^{(2,i,j)},$$
(100)
2. $`i+j=3`$

$$(i+1)\chi _{k_1+1,k_2,k_3}^{(2,i+1,j)}+(j+1)\chi _{k_1,k_2+1,k_3}^{(2,i,j+1)}=J_𝐤^{(2,i,j)},$$
(101)
3. $`i+j=2`$

$`(i+1)\chi _{k_1+1,k_2,k_3}^{(2,i+1,j)}`$ $`+`$ $`(j+1)\chi _{k_1,k_2+1,k_3}^{(2,i,j+1)}+2k_3u_2^{(1)}[\chi _{k_1+1,k_2,k_31}^{(1,i1,j)}`$ (102)
$`+`$ $`\chi _{k_1,k_2+1,k_31}^{(1,i,j1)}]=J^{(2,i,j)}_𝐤,`$
4. $`i+j=1`$

$`(i+1)\chi _{k_1+1,k_2,k_3}^{(2,i+1,j)}`$ $`+`$ $`(j+1)\chi _{k_1,k_2+1,k_3}^{(2,i,j+1)}+2k_3u_2^{(1)}[\chi _{k_1+1,k_2,k_31}^{(1,i1,j)}`$ (103)
$`+`$ $`\chi _{k_1,k_2+1,k_31}^{(1,i,j1)}]k_3\chi ^{(1,i,j)}_{k_1,k_2,k_31}=J^{(2,i,j)}_𝐤,`$
5. $`i+j=0`$

$$\chi _{k_1+1,k_2,k_3}^{(2,1,0)}+\chi _{k_1,k_2+1,k_3}^{(2,0,1)}k_3\chi _{k_1,k_2,k_31}^{(1,0,0)}=J_𝐤^{(2,0,0)}.$$
(104)

In the above equations the convention $`\chi _𝐤^{(2,i,j)}=0`$ if $`i<0`$ or $`j<0`$ is implicitly understood.

According to the system (101)–(104), the coefficients of degree $`k+1`$ with $`i+j=\alpha +1`$ determine those of degree $`k`$ with $`i+j=\alpha `$. For instance, Eq. (100) allows one to determine the coefficients with $`k=4`$ and $`i+j=4`$. The results obtained are inserted into (101) to determine the coefficients of degree $`k=3`$ and $`i+j=3`$. Then we go to Eq. (102) in order to determine those of degree $`k=2`$ and $`i+j=2`$. In general, the process of solution is done according to the following scheme:

$$\{\chi _k^{(2,i+j=4)}\}\{\chi _{k1}^{(2,i+j=3)}\}\{\chi _{k2}^{(2,i+j=2)}\}\{\chi _{k3}^{(2,i+j=1)}\}\{\chi _{k4}^{(2,i+j=0)}\},$$
(105)
where $`\{\chi _k^{(2,\alpha )}\}`$ denotes the set of coefficients of the same degree, namely $`\{\chi _k^{(2,\alpha )}\}\{\chi _𝐤^{(2,i,j)};k_1+k_2+k_3=k,i+j=\alpha \}`$. The expression of $`J_𝐤^{(2,i,j)}`$ is a linear combination of the coefficients $`\chi _𝐤^{(2,i,j)}`$, $`T_4^{(2)}`$, $`T_2^{(2)}`$ and $`p_2^{(2)}`$. Once a row $`k`$ of the chain (105) is solved, it is necessary to go to the following row $`k+2`$. This process continues until all the coefficients of the required profiles are completely determined. The results for the pressure and temperature are

$$p^{(2)}=\frac{3}{10}r^2,T^{(2)}=\frac{34}{175}r^2\frac{1}{240}r^4.$$
(106)
### B.3 Third and fourth orders in $`g`$ ($`s=3,4`$)

The process to third order in $`g`$ is practically the same as the previous one. Since the pressure and the temperature are even functions of $`g`$, only the coefficients of $`u_z`$ appear at this level. The result is

$$u^{(3)}=u_2^{(3)}r^2\frac{\mathrm{4\hspace{0.17em}429}}{\mathrm{50\hspace{0.17em}400}}r^4\frac{1}{\mathrm{2\hspace{0.17em}160}}r^6,$$
(107)
where the coefficient $`u_2^{(3)}=6.0806\mathrm{}`$ is related to the ratio $`A_4/A_2`$.

The calculation algorithm is essentially the same for the higher orders in $`g`$. However, the equations become more and more cumbersome and their solution requires a considerable computational effort. Here we give the final results for the pressure and the temperature to fourth order:

$$p^{(4)}=p_2^{(4)}r^2+p_4^{(4)}r^4+\frac{7}{\mathrm{7\hspace{0.17em}200}}r^6,$$
(108)
$$T^{(4)}=T_2^{(4)}r^2+T_4^{(4)}r^4\frac{\mathrm{30\hspace{0.17em}931}}{\mathrm{13\hspace{0.17em}608\hspace{0.17em}000}}r^6\frac{23}{\mathrm{691\hspace{0.17em}200}}r^8,$$
(109)
where $`p_2^{(4)}=24.160\mathrm{}`$, $`p_4^{(4)}=0.062666\mathrm{}`$, $`T_2^{(4)}=59.511\mathrm{}`$ and $`T_4^{(4)}=0.52856\mathrm{}`$. The high values of $`|p_2^{(4)}|`$ and $`|T_2^{(4)}|`$ suggest that the series (44)–(47) are only asymptotic, in agreement with the situation in the case of the BGK model .

Along with the hydrodynamic fields, we have obtained all the moments of degrees $`k=3`$ and $`k=4`$ to third order in $`g`$. This allows us to get the pressure tensor and the heat flux. Their expressions in cylindrical coordinates are given in Sec. 5.


text:
# Variability of sub-mJy radio sources

## 1 Introduction

Synoptic surveys for Gamma-ray bursts (GRBs), and subsequent ground-based observational follow-up at radio through optical wavelengths, has highlighted the importance of transient celestial phenomena (Masetti 2001). The new parameter space of the transient cosmos has been emphasized in the design of future telescopes, such as the optical Large Synoptic Survey Telescope (Tyson & Angel 2001), and the radio Square Kilometer Array (van Haarlem 1999).

While it is well documented that flat spectrum radio sources can be variable (Aller et al. 1985), the areal density of such sources has not been well quantified through multi-epoch, wide field blind surveys. At high flux density levels ($`>10`$ mJy at 1.4 GHz), one can make a rough estimate of the areal density of variable radio sources by simply assuming that all flat spectrum sources are variable. For instance, the areal density of all sources $`10`$mJy is $`2\times 10^3`$ arcmin-2, and the fraction of flat spectrum sources is about 10$`\%`$, implying an areal density of variable radio sources of $`2\times 10^4`$ arcmin-2 (Gruppioni et al. 1999; White et al. 1997; Windhorst et al. 1985; Hopkins et al. 2000, 2002). This number is consistent with the (null) results of Frail et al. (1994) in their search for highly variable mJy-level sources associated with GRBs.

Source populations at these high flux densities are dominated by AGN. Below about 1 mJy the slope of the source counts flattens, and star forming galaxies are thought to dominate the faint source population (Windhorst et al. 1985; Georgakakis et al. 1999; Hopkins et al. 2000, 2002). Hence, when considering the areal density of variable sub-mJy radio sources, one cannot simply extrapolate the results from high flux density source samples to low flux densities.

Knowledge of the areal density of variable sub-mJy radio sources is critical for setting the back-ground, or ‘confusion’, level for studies of faint variable source populations, such as GRBs (Frail et al. 1997). A recent comparison of the NVSS and FIRST surveys by Levinson et al. (2002) sets a conservative upper limit of $`3\times 10^6`$ arcmin-2 to the areal density of ‘orphan’ GRB radio afterglows (i.e. GRBs for which the $`\gamma `$-ray emission is not beamed toward us) with S$`{}_{1.4}{}^{}6`$ mJy. They also argue that the areal density of radio supernovae will be considerably smaller. While the area of the sky covered by Levinson et al. (2002) was much larger that the study presented herein, their flux density limit was higher than any GRB radio afterglow yet recorded. In this paper we present a smaller area study, but we consider variable sources at flux density levels ($`0.1`$mJy) applicable to typical GRB radio afterglows.

In general, variability of radio sources at the sub-mJy level is an essentially unexplored part of parameter space – a part of parameter space which may fundamentally drive the design of future radio telescopes, such as the SKA (Carilli et al. 2002). Herein we present the first study to delve into this part of parameter space, by exploring systematically the variability of the sub-mJy radio source population at 1.4 GHz. We examine variability on timescales of 17 months and 19 days. Note that the LSST will probe similar variability timescales in the optical, with sampling on weekly to yearly timescales.

## 2 Observations

Observations were made using the using the VLA at 1.4 GHz in the B configuration (maximum baseline = 10km). The region observed is in the Lockman Hole centered at: (J2000) 10h 52m 56.00s, 57d 29 06.0′′. Table 1 summarizes the observations. Column 1 gives the observing date, column 2 gives the observed hour angle range, and column 3 gives the rms in the final image. These observations are part of a larger multiwavelength program to study the evolution of dusty star forming galaxies (Bertoldi et al. in prep).

Standard wide field imaging techniques were employed in order to generate an unaberrated image of the full primary beam of the VLA (FWHM = 32). The absolute flux density scale was set using 3C286. We then generated a CLEAN component model of the field using self-calibrated data taken on Sept. 5, 2002. The data from all days were then self-calibrated in amplitude and phase without gain renormalization using this model. This process should ensure that all the data are on the same flux scale. Images made before and after this process showed that the absolute flux scale changed by at most 1$`\%`$. We also checked to see if variable sources could be removed (or added) due to using a single model to self-calibrate data from all the different days. Components at the 0.1 to 0.2 mJy level were added to the self-calibration model at random positions in the field, and the self-calibration process was repeated. In no case was a new source generated. This gives us confidence that the self-calibration process is robust to small perturbations in the model, i.e. that the problem is over-constrained and that the input self-calibration model is dominated by non-varying sources.

Images for each day were generated using the wide field imaging capabilities in the AIPS task IMAGR (Perley 1999). To remove problems with ‘beam squint’ (slightly different pointing centers for right and left circular polarizations) the right and left polarizations were imaged separately. The images were then summed, weighted by the rms on each image. The final image using data from all the observing days is shown in Figure 1. The rms noise on this image is 7$`\mu `$Jy beam-1 and the restoring CLEAN beam is circular with FWHM = 4.5′′.

## 3 Analysis

We searched for source variations over 17 months by comparing images made in April 2001 with those made in August/September 2002, and also on timescales of 19 days by comparing images made on August 17 and September 5, 2002. We searched for variable sources out to the 10$`\%`$ point of the primary beam, and we limited the analysis to sources with variations $`>2\%`$. Images from all epochs were convolved to 6′′ resolution to mitigate differences that might occur due to the non-linear process of deconvolution. The rms on the images used for variability analysis was 12.5 $`\mu `$Jy for the 17 month comparison, and 17$`\mu `$Jy for the 19 day comparison.

The images from different epochs were both subtracted and averaged. A fractional variability image was then generated by dividing these two images, blanking at the 5$`\sigma `$ level (in absolute value). For the 17 months variability analysis $`\sigma `$ was 17$`\mu `$Jy on the differenced image, and 9$`\mu `$Jy on the averaged image. The corresponding numbers for the 19 day analysis were 26$`\mu `$Jy and 15$`\mu `$Jy, respectively.

The difference image for the 17 month analysis is shown in Figure 2. Some artifacts are seen around the brighter extended sources arising from differences in deconvolution and residual calibration errors. Beyond these artifacts, the difference images are remarkably free of sources. This result gives us confidence that the imaging process (self-calibration and deconvolution) does not generate spurious sources at the $`5\sigma `$ level, and tells us right away that the radio sky is not highly variable at the 100$`\mu `$Jy level.

Variable sources were identified in the divided image. Five variable sources were found at the $`5\sigma `$ level in the 17 month comparison, while four sources were found in the 19 day analysis. We then returned to the original images to find the flux densities of the sources at each epoch. Table 2 lists source positions (columns 1 and 2), flux densities at the two epochs (columns 3 and 4), and the distance from the phase center (column 5) for the 17 month analysis, and Table 3 lists the corresponding values for the 19 day analysis. Note that the source J1051+5734 was seen to vary on both timescales. In fact, this source varied between September 5 and September 9 from 1.53 mJy to 0.71 mJy. The September 2002 value listed in Table 2 is the weighted average of these two measurements.

We next consider the sensitivity of our observations to variable sources at some absolute level, $`\mathrm{\Delta }`$S. The analysis is complicated by the roughly Gaussian roll-off of the primary beam of the VLA, with FWHM $`32^{}`$. Identification of variable sources was done using the non-primary beam corrected maps in order to have uniform noise across the field. Of course, the final flux densities and noise levels for the variable sources were corrected for the primary beam attenuation.

Given a $`\mathrm{\Delta }`$S, one can calculate the maximum primary beam correction, f, for which a change in flux density $`\mathrm{\Delta }`$S could be detected at the 5$`\sigma `$ level, where $`\sigma `$ is the noise at the field center on the difference image: f = 5$`\sigma `$/$`\mathrm{\Delta }`$S. The value of f then sets the distance from the pointing center, R, to which such variation could have been detected, given the primary beam shape of the VLA.

Values for $`\mathrm{\Delta }`$S, f, and R are listed in columns 1, 2, and 3, respectively in Table 4. Column 4 lists the number of variable sources over 17 months that meet these criteria. Column 5 lists the number of sources within the specified radius with flux density, S $``$ $`\mathrm{\Delta }`$S/2. This flux density sets the limit for 100$`\%`$ variability, e.g. a source could be 100$`\mu `$Jy on day 1, and 0 $`\mu `$Jy on day 2, leading to a value on the difference image of 100 $`\mu `$Jy and a value on the average image of 50$`\mu `$Jy. Column 6 lists the areal density of sources with S$``$$`\mathrm{\Delta }`$S/2 on our Lockman hole image, while column 7 lists the corresponding values from the recent study of Fomalont et al. (2002) comprised of a number of different fields.

We have investigated the source sizes using the B array observations at 4.5′′ resolution. Two of the sources are partially resolved. For J1051+5708 Gaussian fitting to the profile yields a peak surface brightness, I$`{}_{\nu }{}^{}=8.92\pm 0.05`$ mJy beam-1, a total flux density, S$`{}_{\nu }{}^{}=10.00\pm 0.09`$ mJy, and a (deconvolved) size of $`1.9^{\prime \prime }\times 1.1^{\prime \prime }`$, with major axis position angle (PA) = $`37^o`$. The corresponding numbers for J1055+5718 are: I$`{}_{\nu }{}^{}=8.82\pm 0.06`$ mJy beam-1, S$`{}_{\nu }{}^{}=12.35\pm 0.13`$ mJy, and $`3.4^{\prime \prime }\times 2.2^{\prime \prime }`$ at PA $`=81^o`$. The rest of the sources are unresolved, with upper limits between 1′′ and 2′′ depending on signal-to-noise.

The nature of this analysis is such that we are not sensitive to very rapid variations, e.g. timescales of minutes or less. For instance, a 50 mJy flare of 1 min duration would average down to about 100$`\mu `$Jy over 7 hours. While our nominal sensitivity would be adequate to detect such an event, the existence of such a transient source in the visibility data would lead to errors in the self-calibration and imaging process which would manifest themselves clearly on the images. Of course, it is possible that such a bright, short timescale event was mis-identified as interference in the data editing process, and removed.

## 4 Discussion

For the analysis on 17 month timescales we could have detected sources with $`\mathrm{\Delta }`$S $`100\mu `$Jy out to a radius of 7.8, but none were detected. This sets an upper limit to the areal density of such sources of $`5\times 10^3`$ arcmin-2. Another interesting point is that within this radius there are 46 sources between 50 and 100 $`\mu `$Jy on the averaged image. Hence, we set an upper limit of about 2$`\%`$ to the fraction of sources in this flux density regime that are highly ($`50\%`$) variable.

These results are grossly consistent with the idea that below 1 mJy at 1.4 GHz the radio source population is dominated by star forming galaxies, as compared to AGN which dominate at high flux densities (Hopkins 2000, 2002). The exact distribution of AGN vs. starbursts vs. other source types as a function of radio flux density is not fully determined at this time, and is an area of active current research (Hopkins et al. 2002; Georgakakis et al. 1999; Richards 2001; Fomalont et al. 2002). As a rough guide we consider the models of Hopkins et al. (2000). They suggests that at $``$ 10 mJy the source population is 90$`\%`$ steep spectrum radio sources, 10$`\%`$ flat spectrum radio sources, and $`1\%`$ star forming galaxies. Again, at high flux density levels the flat spectrum sources correspond to the variable radio source population. At 100 $`\mu `$Jy the models suggests that the proportions change to roughly 80$`\%`$ star forming galaxies, 15$`\%`$ steep spectrum AGN, and 5$`\%`$ flat spectrum AGN.

These models have not considered truly transient source populations, such as GRB radio afterglows, which have timescales of 10’s of days. The statistics of radio afterglows of GRBs are such that in an area of 7.8 radius one would expect to see $`0.005`$ sources above 0.1 mJy at 1.4 GHz at any given time, assuming GRBs are highly beamed (see below). Hence, the fact that we did not see such a source is not surprising. More importantly, the results presented herein allow us to set a limit on the variable source confusion level for GRB radio afterglow searches, again at a flux density level relevant to the observed population. For example, radio searches are needed for localizing an important subclass of afterglows known as ‘dark GRBs’ (e.g., Djorgovski et al. 2001), for which optical emission from the GRB is not detected. It has been suggested that the absence of optical emission from these sources is the result of either dust obscuration, or the Gunn-Peterson effect, i.e. Ly$`\alpha `$ absorption by the neutral intergalactic medium. This latter effect would place the sources at $`z>6.3`$ (Fan et al. 2003). An upper limit to the areal density of $`5\times 10^3`$ arcmin-2 for variable sources at the 100$`\mu `$Jy level implies a lower limit of 200 arcmin2 to the area that can be safely searched at 1.4 GHz before one such source is expected to be detected by chance. For comparison, the typical GRB error circle is $``$30 arcmin2, but larger error circles are not uncommon.

This result can also be used to derive a rough limit on the mean beaming angle for GRBs (Perna & Loeb 1998). From a sample of 25 radio afterglows (Frail et al. in prep) we estimate that 10-25% will be visible above the 100 $`\mu `$Jy level at 1.4 GHz, with an average lifetime of one month. The GRB event rate is approximately 600 per year (Fishman & Meegan 1995) so we expect to find only $`8.4\times 10^8`$ arcmin-2. However, if GRBs are highly beamed, as recent studies seem to suggest (Frail et al. 2001), then our upper limit to areal density of $`5\times 10^3`$ arcmin-2 implies a beaming factor $`f_b<6\times 10^4`$, or a mean jet opening angle $`\theta _j>1^{}`$. This value is not very constraining compared to existing limits, but to our knowledge this is the first time a survey for variability has been done at the appropriate flux level for radio afterglows. More stringent limits will require sensitive, larger area surveys with existing or planned instruments (Totani & Panaitescu 2002).

A final point we consider is cosmic variance. It is possible that the Lockman hole region we have sampled was just statistically-poor in transient sources. A rough indication of the effect of cosmic variance comes from the overall source counts (columns 6,7 in Table 4). The source counts we derive agree to within 20$`\%`$ with those found by Fomalont et al. (2002) in other areas of the sky. In general, it has been found that for 30 fields-of-view the maximum field-to-field scatter in the sub-mJy source counts is about a factor two (Fomalont et al. 2002). We consider this an upper limit to the effect cosmic variance has on the areal density of variable source presented herein.

The National Radio Astronomy Observatory (NRAO) is operated by Associated Universities, Inc. under a cooperative agreement with the National Science Foundation. We thank F. Owen, E. Fomalont, R. Becker, and J. Condon for useful comments concerning this work, and the referee for careful and insightful criticism.

af = 5$`\sigma /\mathrm{\Delta }`$S, where $`\sigma =17\mu `$Jy = rms on difference maps at the field center.

bR = distance from pointing center out to which $`\mathrm{\Delta }`$S can be detected to $`5\sigma `$.

cNvar = number of variable sources $`\mathrm{\Delta }`$S within given radius over 17 months.

dNsrcs = number of sources with flux density, $`\mathrm{S}\mathrm{\Delta }\mathrm{S}/2`$ within R. This S sets the 100$`\%`$ variability detection limit.

eNsrcs/area = number of sources with $`\mathrm{S}\mathrm{\Delta }\mathrm{S}/2`$ per arcmin2 from our Lockman Hole field.

fNsrcs/area = number of sources with $`\mathrm{S}\mathrm{\Delta }\mathrm{S}/2`$ per arcmin2 from the study of Fomalont et al. (2002): N($`\mathrm{S}_{\mathrm{mJy}})=0.026\times \mathrm{S}^{1.1}\mathrm{arcmin}^2`$.

Figure Captions

FIG. 1.— Radio image of the Lockman Hole at 1.4 GHz with a resolution of 4.5′′. The rms noise on the image is 7$`\mu `$Jy. The grayscale range is from -0.15 mJy to 0.25 mJy.

FIG. 2.— The difference image between observations in April 2001 and Aug/Sept 2002 at 6′′ resolution. The grayscale range is from -0.2 mJy to 0.2 mJy.


text:
# References

Supersymmetric grand unified theories (GUTs) provide an especially attractive framework for physics beyond the standard model (and MSSM), and it is therefore natural to ask if there exists in this framework a compelling, perhaps even an intimate connection with inflation. In ref.  one possible approach to this question was presented. In its simplest realization, inflation is associated with the breaking at scale $`M`$ of a grand unified gauge group $`G`$ to $`H`$. Indeed, inflation is ‘driven’ by quantum corrections which arise from the breaking of supersymmetry by the vacuum energy density in the early universe. The density fluctuations, it turns out, are proportional to $`(M/M_{\mathrm{Planck}})^2`$, where $`M_{\mathrm{Planck}}1.2\times 10^{19}`$ GeV denotes the Planck mass. From the variety of $`\delta T/T`$ measurements, especially by the Wilkinson Microwave Anisotropy Probe (WMAP) , the symmetry breaking scale $`M`$ is of order $`10^{16}`$ GeV, essentially identical to the scale of supersymmetric grand unification.

Because of the logarithmic radiative corrections that drive inflation, the spectrum of scalar density fluctuations turns out to be essentially flat. For the simplest models, the scalar spectral index was found to be $`n_s=0.98(\pm 0.01)`$ , in excellent agreement with a variety of observations  including the recent WMAP data. The variation $`dn_s/d\mathrm{ln}k`$ of the spectral index is found to be small ($`10^3`$).

In some recent papers  it was shown how the above scheme can be extended to five dimensional supersymmetric models. There are good reasons for discussing such models. Consider, for instance, the case of $`G=SO(10)`$ (or $`SU(5)`$) in four dimensions. The presence of dimension five baryon number violating operators mediated through Higgsino exchange implies in the ‘minimal’ scheme a proton life time $`\tau _{pK^+\overline{\nu }}10^{30\pm 2}`$ yrs. This may be in conflict with the recent lower bounds ($`\tau _p>1.9\times 10^{33}`$ yrs) for $`pK^+\overline{\nu }`$ determined by the Superkamiokande experiment . There are other serious issues such as the notorious doublet-triplet (DT) splitting problem, which have led people to investigate five (and higher) dimensional theories compactified on suitable orbifolds that provide a relatively painless way of implementing the DT splitting. Furthermore, dimension five proton decay can be easily eliminated which is an attractive feature of the five dimensional framework.

In this paper we present a realistic model of inflation based on five dimensional supersymmetric $`SO(10)`$ compactified on an orbifold $`S^1/(Z_2\times Z_2^{})`$. There are two fixed points (branes) where the gauge symmetries are $`SO(10)`$ and $`SU(4)_c\times SU(2)_L\times SU(2)_R`$, respectively . Through the spontaneous breaking of $`SO(10)`$ to $`SU(5)`$, the effective low energy symmetry corresponds to the MSSM gauge group. The inflationary scenario will be associated with the symmetry breaking $`SO(10)SU(5)`$. In its simplest realization the scalar spectral index $`n_s=0.98\pm 0.01`$. The variation of $`n_s`$ with respect to the wave number $`k`$ is small but in principle measurable ($`dn_s/d\mathrm{ln}k10^3`$), which will be tested by ongoing and future observations. The WMAP data combined with other observations appears to prefer, but does not require, a far more significant running, $`dn_s/d\mathrm{ln}k=0.031_{0.018}^{+0.016}`$. Our model predicts a negligible tensor-to-scalar ratio $`r10^8`$. After inflation is over, the inflaton (which belongs to $`\mathrm{𝟏𝟔}`$, $`\overline{\mathrm{𝟏𝟔}}`$ of $`SO(10)`$) decays into right-handed neutrinos whose out of equilibrium leads to the observed baryon asymmetry via leptogenesis . The symmetry breaking scale $`M`$ is estimated from inflation to be of order $`10^{16}`$ GeV, which leads to righ-handed neutrino masses of the correct magnitude, of order $`10^{14}`$ GeV or less, that can yield a mass spectrum for the light neutrinos suitable for neutrino oscillations .

In this paper we also provide a new resolution of the well known monopole problem. In the class of models discussed, unless care is exercised, superheavy monopoles can be produced at the end of inflation leading to cosmological disaster. In ref. , the problem was circumvented by including sizable higher order (non-renormalizable) terms in the superpotential, such that the GUT symmetry is broken along an inflationary trajectory. In our new scenario we show how the problem is solved in a five dimensional framework without invoking non-renormalizable terms.

The four dimensional inflationary model is best illustrated by considering the following superpotential which allows the breaking of a gauge symmetry $`G`$ down to $`H`$, keeping supersymmetry intact :

$`W_{\mathrm{infl}}=\kappa S(\varphi \overline{\varphi }M^2).`$ (1)
Here $`\varphi `$ and $`\overline{\varphi }`$ represent superfields whose scalar components acquire non-zero vacuum expectation values (VEVs), which break $`G`$ to $`H`$. The singlet superfield $`S`$ provides the scalar field that drives inflation. Note that by invoking a suitable $`R`$ symmetry $`U(1)_R`$, the form of $`W`$ is unique at the renormalizable level. For example, $`W`$ and $`S`$ can be assigned an $`R`$-charge of unity, while the $`R`$-charges of $`\varphi `$, $`\overline{\varphi }`$ are zero. It is gratifying to realize that $`R`$ symmetries naturally occur in (higher dimensional) supersymmetric theories and can be appropriately exploited.

From $`W`$, it is straightforward to show that the supersymmetric minimum corresponds to non-zero (and equal in magnitude) VEVs for $`\varphi `$ and $`\overline{\varphi }`$, while $`S=0`$. (After supersymmetry breaking $`\stackrel{`}{a}`$ la $`N=1`$ supergravity, $`S`$ acquires a VEV of order $`m_{3/2}`$ (gravitino mass)).

An inflationary scenario is realized in the early universe with both $`\varphi `$, $`\overline{\varphi }`$ and $`S`$ displaced from their present day minima. Thus, for $`S`$ values in excess of the symmetry breaking scale $`M`$, $`\varphi `$, $`\overline{\varphi }`$ VEVs vanish, the gauge symmetry is restored, and a potential energy density proportional to $`M^4`$ dominates the universe. With supersymmetry thus broken, there are radiative corrections from the $`\varphi `$-$`\overline{\varphi }`$ supermultiplets that provide logarithmic corrections to the potential which drives inflation. In one loop approximation ,

$`V\kappa ^2M^4\left[1+{\displaystyle \frac{\kappa ^2𝒩}{32\pi ^2}}\left(2\mathrm{l}\mathrm{n}{\displaystyle \frac{\kappa ^2|S|^2}{\mathrm{\Lambda }^2}}+(z+1)^2\mathrm{ln}(1+z^1)+(z1)^2\mathrm{ln}(1z^1)\right)\right],`$ (2)
where $`z=x^2=|S|^2/M^2`$, $`𝒩`$ is the dimensionality of the representations to which $`\varphi `$, $`\overline{\varphi }`$ belong, and $`\mathrm{\Lambda }`$ denotes a renormalization mass scale. The logarithmic loop corrections in Eq. (2) enable the inflaton field to slowly roll down to the supersymmetric vacuum state. From Eq. (2) the microwave CMB anisotropy on the Hubble scale $`l`$ is found to be

$`\left({\displaystyle \frac{\delta T}{T}}\right)_l{\displaystyle \frac{8\pi }{\sqrt{𝒩}}}\left({\displaystyle \frac{N_l}{45}}\right)^{1/2}\left({\displaystyle \frac{M}{M_{\mathrm{Planck}}}}\right)^2x_l^1y_l^1f(x_l^2)^1.`$ (3)
Here, $`y_lx_l(17/12x_l^2+\mathrm{})`$, $`f(x_l^2)^11/x_l^2`$, for $`S_l`$ sufficiently larger than $`M`$, and $`N_l5060`$ denotes the e-foldings needed to resolve the horizon and flatness problems.

Comparison of the expression for $`\delta T/T`$ in Eq. (3) with the WMAP result shows that the gauge symmetry breaking scale $`M`$ is around $`10^{16}`$ GeV , which is tantalizingly close to the GUT scale inferred from the evolution of the MSSM gauge couplings. Thus, it is natural to embed this kind of inflationary scenario within a GUT framework. However, in this case we must make sure that cosmological problems associated with topological defects such as monopoles do not arise. In addition, while constructing a realistic inflationary model based on supersymmetric GUT, we would also like to resolve the notorious DT splitting problem. These, as we will see, are most easily carried out in a five dimensional framework.

The scalar spectral index $`n_s`$ is approximately given by

$`n_s1{\displaystyle \frac{1}{N}}`$ (4)
where $`N`$ denotes the number of e-foldings experienced by the scale under consideration. For the horizon scale, in particular, $`N_h5060`$, so that

$`n_s^{(h)}=0.98\pm 0.01.`$ (5)
It should be noted that the inclusion of supergravity corrections can, in some cases, lead to a spectral index larger than unity \[for a recent discusion and additional references, see ref. .\]. The galactic scale corresponds to $`N_g4050`$, and given Eq. (4), we conclude that the variation of $`n_s`$ with $`k`$ is quite tiny, $`dn_s/d\mathrm{ln}k10^3\mathrm{few}\times 10^4`$. It would be interesting to test this prediction against the ongoing and future observations.

As shown in refs. , a combination of the the gravitino constraint on the reheat temperature ($`T_R10^{10}`$ GeV ) as well as leptogenesis requires that the dimensionless superpotential coupling $`\kappa 10^3`$. Thus, the vacuum energy density during inflation is of order $`10^6M_{\mathrm{GUT}}^4`$, so that the tensor-to-scalar ratio $`r10^8`$, which could be hard to detect in any forseeable experiment.

After inflation is over the universe converts to a radiation dominated epoch through the superpotential couplings $`\gamma _{ij}\overline{\varphi }\overline{\varphi }\mathrm{𝟏𝟔}_i\mathrm{𝟏𝟔}_j/M_P`$, where $`\mathrm{𝟏𝟔}_i`$ ($`i=1,2,3`$) denote the three chiral families of $`SO(10)`$ (with $`R`$-charge $`=1/2`$), $`\gamma _{ij}`$ is a dimensionless coupling, and $`M_P`$ ($`2.4\times 10^{18}`$ GeV) is the reduced Planck mass. That is, the inflaton decay produces right-handed neutrinos whose out of equilibrium decay produces the observed baryon asymmetry via leptogenesis along the lines previously discussed in .

Next we present a realistic five dimensional $`SO(10)`$ model in which the inflationary scenario described by the superpotential $`W`$ in Eq. (1) can be realized. We assume compactification on an orbifold $`S^1/(Z_2\times Z_2^{})`$, such that on the two fixed points (branes) we have the gauge symmetries $`SO(10)`$ and $`SU(4)_c\times SU(2)_L\times SU(2)_R`$ respectively. To realize the MSSM gauge group at low energies, we introduce two pairs of the Higgs hypermultiplets $`\mathrm{𝟏𝟔}_H`$ and $`\overline{\mathrm{𝟏𝟔}}_H`$ in the bulk with $`Z_2\times Z_2^{}`$ parities,

$`\mathrm{𝟏𝟔}_H=(\overline{\mathrm{𝟒}},\mathrm{𝟏},\mathrm{𝟐})_H^{++}+(\mathrm{𝟒},\mathrm{𝟐},\mathrm{𝟏})_H^+,`$ (6)
$`\mathrm{𝟏𝟔}_{}^{𝐜}{}_{H}{}^{}=(\overline{\mathrm{𝟒}},\mathrm{𝟏},\mathrm{𝟐})_H^c+(\mathrm{𝟒},\mathrm{𝟐},\mathrm{𝟏})_H^{c+},`$ (7)
$`\overline{\mathrm{𝟏𝟔}}_H=(\mathrm{𝟒},\mathrm{𝟏},\mathrm{𝟐})_H^{++}+(\overline{\mathrm{𝟒}},\mathrm{𝟐},\mathrm{𝟏})_H^+,`$ (8)
$`\overline{\mathrm{𝟏𝟔}}_{}^{𝐜}{}_{H}{}^{}=(\mathrm{𝟒},\mathrm{𝟏},\mathrm{𝟐})_H^c+(\overline{\mathrm{𝟒}},\mathrm{𝟐},\mathrm{𝟏})_H^{c+}.`$ (9)
The relevant superpotentials on the two branes, B1 ($`SO(10)`$ brane) and B2 ($`SU(4)_c\times SU(2)_L\times SU(2)_R`$ brane) are:

$`W_{B1}=\kappa S\left(\mathrm{𝟏𝟔}_H\overline{\mathrm{𝟏𝟔}}_HM_1^2\right),`$ (10)
$`W_{B2}=\kappa S\left(c_1H^c\overline{H}^c+c_2\mathrm{𝟏𝟏}^{}M_2^2\right)+c_3\mathrm{\Sigma }H^c\overline{H}^c,`$ (11)
where $`H^c(\overline{\mathrm{𝟒}},\mathrm{𝟏},\mathrm{𝟐})_H^{++}`$, $`\overline{H}^c(\mathrm{𝟒},\mathrm{𝟏},\mathrm{𝟐})_H^{++}`$, and $`c_1`$, $`c_2`$, $`c_3`$ are dimensionless couplings. In $`W_{B2}`$, we exhibit only the chiral multiplets with $`(++)`$ parities of $`\mathrm{𝟏𝟔}_H`$, $`\overline{\mathrm{𝟏𝟔}}_H`$ which contain massless modes, since the heavy KK modes would be decoupled. Since the inflaton $`S`$ is a bulk superfield, it participates in both superpotentials. In Eq. (11), a pair of singlet superfields $`\mathrm{𝟏}`$, $`\mathrm{𝟏}^{}`$ and a superfield $`\mathrm{\Sigma }`$ in the adjoint representation $`(\mathrm{𝟏𝟓},\mathrm{𝟏},\mathrm{𝟏})`$ with suitable $`U(1)_R`$ charges are introduced on B2.

During inflation, $`S`$ and $`\mathrm{\Sigma }`$ develop VEVs ($`S>M_1,M_2`$), while $`\mathrm{𝟏𝟔}_H=\overline{\mathrm{𝟏𝟔}}_H=H^c=\overline{H}^c=\mathrm{𝟏}=\mathrm{𝟏}^{}=0`$. As shown in refs. , positive vacuum energies localized on the branes could trigger exponential expansion of the three space, in the presence of a brane-localized Einstein-Hilbert term. Due to a non-zero VEV of $`\mathrm{\Sigma }`$ during inflation, the $`SU(4)_c`$ factor in $`SU(4)_c\times SU(2)_L\times SU(2)_R`$ is spontaneously broken to $`SU(3)_c\times U(1)_{BL}`$, and the accompanying monopoles are inflated away.

In this brane model, $`\mathrm{𝟏𝟔}_H`$, $`\overline{\mathrm{𝟏𝟔}}_H`$ on B1, and $`\mathrm{𝟏}`$, $`\mathrm{𝟏}^{}`$ on B2 play the role of $`\varphi `$, $`\overline{\varphi }`$ in Eq. (1). With the (localized) VEVs of the scalar components of $`\mathrm{𝟏𝟔}_H`$, $`\overline{\mathrm{𝟏𝟔}}_H`$ along the $`SU(5)`$ singlet direction (i.e. $`\nu _H^c`$, $`\overline{\nu }_H^c`$) at B1 after inflation, the $`SO(10)`$ gauge symmetry breaks to $`SU(5)`$. On the other hand, at B2 only the singlets $`\mathrm{𝟏}`$, $`\mathrm{𝟏}^{}`$ rather than $`H^c`$, $`\overline{H}^c`$ develop VEVs at the minimum of the potential. Since $`\mathrm{\Sigma }`$ becomes heavy by VEVs of $`\mathrm{𝟏𝟔}_H`$, $`\overline{\mathrm{𝟏𝟔}}_H`$ on B1, the VEV $`\mathrm{\Sigma }`$ vanishes after inflation, and so the symmetry $`SU(4)_c\times SU(2)_L\times SU(2)_R`$ on B2 is restored. Consequently, the effective low energy theory after inflation is the desired MSSM ($`=\{SU(5)\}\{SU(4)_c\times SU(2)_L\times SU(2)_R\}`$). We note that the symmetry breaking process $`SU(3)_c\times U(1)_{BL}\times SU(2)_L\times SU(2)_RSU(3)_c\times SU(2)_L\times U(1)_Y`$ does not create any unwanted topological defects such as monopoles, and so we have formulated a realistic 5D model in which the monopole problem is solved without introducing non-renormalizable terms.

While the Goldsotne fields $`(\overline{\mathrm{𝟑}},\mathrm{𝟏})_{2/3}^{++}`$, $`(\mathrm{𝟏},\mathrm{𝟏})_1^{++}`$ (also $`(\mathrm{𝟑},\mathrm{𝟏})_{2/3}^{++}`$, $`(\mathrm{𝟏},\mathrm{𝟏})_1^{++}`$) of $`H^c`$ ($`\overline{H}^c`$) are absorbed by the appropriate gauge bosons, the superhiggs mechanism leaves intact the massless supermultiplets $`(\overline{\mathrm{𝟑}},\mathrm{𝟏})_{1/3}`$, $`(\mathrm{𝟑},\mathrm{𝟏})_{1/3}`$, which can acquire masses of order $`m_{3/2}`$ from their couplings to $`S`$ after supersymmetry breaking. To eliminate this pair from the low energy theory, we can introduce on B1 a 10-plet with couplings $`\mathrm{𝟏𝟔}_H\mathrm{𝟏𝟔}_H\mathrm{𝟏𝟎}`$ and $`\overline{\mathrm{𝟏𝟔}}_H\overline{\mathrm{𝟏𝟔}}_H\mathrm{𝟏𝟎}`$ (thus, $`\mathrm{𝟏𝟎}`$ has an $`R`$-charge of unity), and/or a $`(\mathrm{𝟔},\mathrm{𝟏},\mathrm{𝟏})`$ field ($`D`$) on B2 with couplings $`H^cH^cD`$ and $`\overline{H}^c\overline{H}^cD`$. Then, the pair acquires superheavy masses proportional to $`\nu _H^c`$ or $`\overline{\nu }_H^c`$, and the low energy spectrum is precisely the MSSM one.

Note that we introduced the Higgs 16-plets in the bulk rather than on the $`SO(10)`$ brane B1 in order to avoid unwanted states associated with the pseudo-Goldstone symmetry of the superpotential. Recall that the orbifold compactification breaks $`SO(10)`$ down to $`SU(4)_c\times SU(2)_L\times SU(2)_R`$.

To resolve the DT splitting problem, the Higgs 10-plet ( or $`(\mathrm{𝟏},\mathrm{𝟐},\mathrm{𝟐})`$) should be introduced in the bulk (on B2). By suitable $`Z_2\times Z_2^{}`$ parity assignments, the MSSM Higgs doublets are kept light, while the color triplets become superheavy.

In summary, we have taken the approach that a satisfactory inflationary scenario should:

(i) resolve the flatness and horizon problem;

(ii) resolve cosmological problems associated with topological defects;

(iii) give rise to the observed $`\delta T/T`$ fluctuations;

(iv) provide a satisfactory explanation of the origin of the observed baryon asymmetry;

(v) be well grounded in particle physics.

While four dimensional $`SO(10)`$ models of inflation are hard to construct, especially if a resolution of DT splitting problem is also desired, things are much easier if we consider five dimensional $`SO(10)`$. In the model we have discussed, the gauge symmetry during inflation is broken to $`SU(3)_c\times U(1)_{BL}\times SU(2)_L\times SU(2)_R`$ so that monopoles are inflated away. When inflation ends, the unbroken gauge symmetry turns out to be $`SU(3)_c\times SU(2)_L\times U(1)_Y`$. The scalar spectral index $`n_s`$ is very close (or even equal ) to unity, $`dn_s/d\mathrm{ln}k10^3`$, and the tensor-to-scalar ratio $`r10^8`$. The observed baryon asymmetry naturally follows from leptogenesis.

Acknowledgments

We thank Zurab Tavartkiladze and Yasunori Nomura for useful discussions. The work is partially supported by DOE under contract number DE-FG02-91ER40626.