format: text

text:
# Abstract

## Abstract

In this paper we show how, under certain restrictions, the hydrodynamic equations for the freely evolving granular fluid fit within the framework of the time dependent Landau–Ginzburg (LG) models for critical and unstable fluids (e.g. spinodal decomposition). The granular fluid, which is usually modeled as a fluid of inelastic hard spheres (IHS), exhibits two instabilities: the spontaneous formation of vortices and of high density clusters. We suppress the clustering instability by imposing constraints on the system sizes, in order to illustrate how LG-equations can be derived for the order parameter, being the rate of deformation or shear rate tensor, which controls the formation of vortex patterns. From the shape of the energy functional we obtain the stationary patterns in the flow field. Quantitative predictions of this theory for the stationary states agree well with molecular dynamics simulations of a fluid of inelastic hard disks.

## Key words

Granular fluid; instabilities; pattern formation; hydrodynamic equations; time dependent Landau-Ginzburg theory.

## 1 Introduction

Granular matter  consists of small or large macroscopic particles. When out of equilibrium, its dynamics is controlled by dissipative interactions, and distinguished in quasi-static flows or granular solids on the one hand, and rapid flows or granular fluids  on the other hand.

Typical realizations of granular solids are sand piles, avalanches, Saturn’s rings, grain silos, stress distributions. Here particles remain essentially in contact, and the dynamics is controlled by gravity, friction and surface roughness. In this paper we concentrate on granular fluids. Typical examples are driven granular flows, such as Poisseuille flow , vibrated beds , or rapid flows with some form of continuous energy input . Here the dynamics is controlled by inelastic binary collisions, separated by ballistic motion of the particles. The forces are of short range and repulsive, and the system is frequently modelled as a collection of smooth inelastic hard spheres (IHS)  of diameter $`\sigma `$ and mass $`m`$. Momentum is conserved during collisions, which makes the system a fluid, but energy is not conserved. In a collision, on average, a fraction $`ϵ`$ of the relative kinetic energy of the colliding pair is lost, where $`ϵ`$ is referred to as the degree of inelasticity. In the literature  $`ϵ=1\alpha ^2`$, usually expressed in terms of the coefficient of restitution $`\alpha `$. Its detailed definition does not concerns us here.

Here we focus on the idealized limiting case of a freely evolving rapid flow without energy input and with nearly elastic collisions, and therefore slowly cooling. This system shows  an interesting instability. When prepared in a spatially homogeneous equilibrium state, the system does not stay there, but slowly develops patterns, both in the flow field (vortices), and in the density field (clusters), the so called clustering instability. For the two-dimensional case the analogies with spinodal decomposition have already been pointed out in the literature .

The search for the proper macroscopic description of unstable granular fluids has been pursued by many authors \[1-13,18,20,22-32\]. Recently, two new points of view have been presented, namely by Ben-Naim et al.  and by Soto et al. . The first one needs to be discussed in more detail because macroscopic equations for granular fluids like the Burgers equation appear in that paper, as well as in the present one. These authors conjecture that the (vector) Burgers equation describes the flow velocity $`𝐮(𝐫,t)`$ of a granular fluid, or at least of a dilute granular gas –the authors are not very explicit on this point– on large space and time scales for arbitrary values of the inelasticity ($`0<ϵ1`$) in $`d`$dimensions, i.e.

$$_t𝐮=𝐮𝐮+\nu ^2𝐮,$$
(1)
where $`\nu `$ is the kinematic viscosity, and the solutions of interest have to satisfy $`\times 𝐮=\mathrm{𝟎}`$. An interesting aside is that the rotation-free solution of the vector Burgers equation can be expressed as the gradient of a scalar field, $`𝐮=h(𝐫,t)`$, where $`h(𝐫,t)`$ satisfies the equation,

$$_th=\nu ^2h+\frac{1}{2}|h|^2.$$
(2)
This equation is the famous KPZ-equation, named after Khadar, Parisi and Zhang . It describes the growth dynamics of solid surfaces, where $`h`$ is the height function.

It is well known that the Burgers equation in the inviscid limit ($`\nu 0`$) describes the ”sticky dust” or adhesion model of perfectly inelastic point particles with $`ϵ=1`$ . The conjectured validity of the Burgers equation for the inelastic hard sphere fluid implies a universality hypothesis, i.e. the large $`(𝐫,t)`$behavior of a granular gas is in the universality class of sticky dust, independent of its inelasticity $`ϵ`$. Several implications of this conjecture in dimensions $`d2`$ have been criticized in Ref . For the one-dimensional case the authors present rather convincing evidence based on molecular dynamics (MD) simulations using $`N=10^6`$ particles, and on scaling arguments. In addition, Boldyrev  has studied the one-dimensional randomly driven Burgers equation with a pressure term $`(p)/\rho `$ with $`pn^a`$ included. He has shown that properties, like the structure factor or pair correlation function, are not affected by the pressure term if $`a<2`$. On account of this argument, the $`p`$-term in the one-dimensional Navier-Stokes equation can be neglected on the largest $`(r,t)`$ scales for a dilute gas in the inviscid limit, making the Burgers equation the appropriate macroscopic equation. Note that in the freely cooling granular fluid the viscosity decreases as $`\sqrt{T}`$ as $`T0`$. In higher dimensions the systems used in MD simulations ($`N=5\times 10^4`$ in two dimensions  and $`N=10^6`$ in three dimensions ) are too small to draw any conclusions on large scale behavior. Here analytic or scaling arguments to support the conjecture are lacking.

Although equations of the form (1) will frequently appear in the present paper, we emphasize that Ref.  refers only to the largest possible scales (where the thermodynamic limit has been taken). In this paper we explicitly restrict ourselves to small systems in order to suppress the clustering instability. This has been done to simplify the problem. Moreover, Ref.  only refers to dilute granular gases, whereas this paper covers also liquid densities.

A second new development, which is somewhat similar to ours, is given by Soto et al. . These authors also study granular fluids, contained in systems that are sufficiently small, such that the clustering instability is suppressed. In these small systems the growth of vortices is very slow. Consequently, the remaining hydrodynamic modes are enslaved by the slowly growing unstable vorticity modes and the amplitude of these vorticity modes remains very small. Under this condition, they obtained the amplitude equations for the slowest vorticity modes, which can be derived from a potential function.

The question of interest in the present paper is: can the models for granular fluids be fitted into the generic classification of Landau–Ginzburg–type models, as given by Hohenberg and Halperin , to describe critical dynamics and hydrodynamic instabilities? The goal of this article is to illustrate how, under certain restrictions, the standard nonlinear hydrodynamic equations for the IHS fluid  can be cast into a Landau–Ginzburg–type equation of motion for the order parameter, which can be derived from an energy functional and, more specifically, to point out which terms in the original hydrodynamic equations are responsible for the quartic terms in the Landau-Ginzburg energy functional.

The plan of the paper is as follows. In Sec.2, we start with the hydrodynamic equations. The decay of the total energy at short times and the results of a linear stability analysis are briefly reviewed. In Sec.3, we introduce an assumption of incompressible flows under certain restrictions on system size or time regime. Then, under these assumptions, the hydrodynamic equations are reduced to a closed equation for a scaled flow field. It is shown in Sec.4 that this equation for a scaled flow field can be cast into the form of a time–dependent Landau–Ginzburg equation for an appropriate order parameter. The shape of the energy functional is discussed and possible stationary solutions are presented. Finally, in Sec.5, we make a quantitative comparison of the theoretical predictions at large times with molecular dynamics simulations of inelastic hard disks. We end with some conclusions in Sec.6.

## 2 Dynamic Equations and Instabilities

The starting point are the hydrodynamic equations. They are not only needed here to recapitulate our present theoretical understanding of this system, but also to formulate the new extensions to be discussed in this paper.

The macroscopic time evolution of the IHS fluid on large spatial and temporal scales  can be described by the nonlinear hydrodynamic equations for the local density $`n(𝐫,t)`$, the local flow field $`𝐮(𝐫,t)`$ and the local temperature $`T(𝐫,t)`$, supplemented with a sink term $`\mathrm{\Gamma }`$ accounting for the energy loss through inelastic collisions, i.e.

$`D_tn`$ $`=`$ $`n𝐮,`$
$`D_t𝐮`$ $`=`$ $`{\displaystyle \frac{1}{mn}}p+2\nu \text{D},`$
$`D_tT`$ $`=`$ $`{\displaystyle \frac{2p}{dn}}𝐮+b_T^2T+2b_{}\text{D}:\text{D}\mathrm{\Gamma }.`$ (3)
In this paper the inelasticity is always assumed to be small. For later convenience the macroscopic equations are given for a $`d`$-dimensional systems. Here $`D_t=_t+𝐮`$ is the time derivative in a comoving frame. The local energy density of the IHS fluid is $`e=\frac{1}{2}mnu^2+\frac{d}{2}nT`$, and $`p`$ is the pressure. The shear rate $`D_{\alpha \beta }`$ is the symmetrized dyadic, $`\{_\alpha u_\beta \}`$, which is also made traceless. The coefficient $`b_T=2\kappa /dn`$ is proportional to the heat conductivity $`\kappa `$, and $`b_{}=2m\nu /d`$ to the shear viscosity $`\nu `$. For simplicity of presentation the bulk viscosity has been set equal to zero, and the transport coefficients $`\nu `$ and $`\kappa `$, which depend on the local density and temperature, are taken at some fixed reference values, $`\overline{n}(t)`$ and $`\overline{T}(t)`$, to be specified later. The four terms in the energy balance equation account for work done by the pressure, for heat conduction, for nonlinear viscous heating and collisional dissipation. Gradients in the flow velocity considerably slow down the collisional cooling process through viscous heating.

On the basis of kinetic theory one can derive that the rate of collisional energy loss, $`\mathrm{\Gamma }=2\gamma _0\omega T`$, is proportional to the collision frequency $`\omega `$ multiplied by the fraction of energy $`ϵT`$ lost per collision , where $`\gamma _0=ϵ/2d=(1\alpha ^2)/2d`$. In general, the collision frequency $`\omega (T)`$ is proportional to the root mean square velocity $`v_0=\sqrt{2T/m}`$, and its explicit form for hard sphere fluids can be found in Refs. . When the system is prepared initially in a homogeneous equilibrium state, it evolves at short times in a spatially homogeneous cooling state (HCS) with a time dependent temperature. Combination of Eqs. (3) and the expression for $`\mathrm{\Gamma }`$ given above, yields $`_tT\gamma _0T^{3/2}`$. This leads to Haff’s homogeneous cooling law  for the mean energy per particle in the IHS fluid,

$$E(t)=\frac{d}{2}T(t)=\frac{E_0}{(1+\gamma _0\omega _0t)^2}=E_0\mathrm{exp}(2\gamma _0\tau ),$$
(4)
which is needed for later comparison. Here $`E_0=(d/2)T_0`$ is the initial energy, and $`t_0=1/\omega _0`$ with $`\omega _0=\omega (T_0)`$ is the mean free time in the initial state. In general the number of collisions per particle $`\tau (t)`$ in a time $`t`$ is defined through $`d\tau =\omega (T(t))dt`$. In the HCS integration of this relation yields $`\mathrm{exp}(\gamma _0\tau )=(1+\gamma _0\omega _0t)`$.

However, this state is unstable against spatial fluctuations in density $`n(𝐫,t)`$, temperature $`T(𝐫,t)`$ and flow velocity $`𝐮(𝐫,t)`$. The present theoretical understanding of these instabilities  is essentially based on a linear stability analysis of the hydrodynamic fluctuations in the density, $`\delta n=n\overline{n}`$, temperature, $`\delta T=T\overline{T}`$, and flow velocity $`𝐮`$. This is done by using the rescaled Fourier modes $`\delta n_𝐤`$, $`\delta \stackrel{~}{T}_𝐤=\delta T_𝐤/\overline{T}`$ and $`\delta \stackrel{~}{𝐮}_𝐤𝐮_𝐤/\sqrt{\overline{T}}`$, where an overline denotes a spatial average, $`\overline{a}=(1/V)_V𝑑𝐫a(𝐫)`$, and $`\overline{T}`$ is the global granular temperature. Fourier transforms are defined as $`f_𝐤=_V𝑑𝐫e^{i𝐤𝐫}f(𝐫)`$. The rescaled eigenmodes are described by $`\delta a_𝐤(\tau )=\mathrm{exp}(z_\lambda \tau )\delta a_𝐤(0)`$. The exponential growth rates of unstable ($`z_\lambda (k)>0`$) and stable ($`z_\lambda (k)<0`$) modes are shown in Fig. 1 as a function of the wave number $`k`$.

This figure shows that the transverse flow field $`𝐮_𝐤`$ or shear mode ($`\lambda =`$) with a wave number $`k<k_{}^{}`$ is unstable, and develops vortices. On the other hand density fluctuations couple weakly, in order $`𝒪(k)`$, to the heat modes ($`\lambda =H`$), and $`z_H(k)`$ in Fig.1 shows that these fluctuations are unstable in the range $`k<k_H^{}`$, and linearly stable in the range $`k>k_H^{}`$ , i.e. remain at thermal noise level. The stability thresholds $`k_{}^{}`$ and $`k_H^{}`$ are defined as the the root of $`z_\lambda (k)=0`$ for $`\lambda =\{,H\}`$, and are marked as black dots in the figure. The figure also shows that the growth rate $`z_{}(k)`$ for the vorticity mode is much larger than the growth rate for the heat mode $`z_H(k)`$, which couples to the density fluctuations. This explains why vortices appear long before the density clusters start to appear.

An intuitive explanation for the appearance and growth of large scale vortices is that a binary collision destroys a fraction of the kinetic energy of relative motion of the colliding pair. The cumulative effect of many successive collisions is that they make the particles move locally in a more parallel and coherent fashion. This creates locally patches of vorticity. These patches grow in size by selective suppression (stronger damping) of short wavelength modes .

An intuitive explanation of the appearance of high density clusters goes as follows . Suppose one prepares the system initially in a spatially homogeneous equilibrium state, and there occurs locally a spontaneous negative pressure fluctuation (‘depression’). The resulting particle flow from the surroundings tries to compensate for the local depression. This creates an excess local density, which increases the collision frequency, and in turn decreases the temperature. This creates again a depression, and the process keeps repeating itself, thus creating cold dense clusters, surrounded by a hot dilute gas. Moreover, these arguments also suggest that the pressure fluctuations and the corresponding gradients remain substantially smaller than those in density and temperature. In fact it would be very interesting if one could test this suggestion by measuring the pressure locally, i.e. in cells of the typical size of a density cluster, through molecular dynamics simulations, or Direct Monte Carlo Simulations (DSMC) of the Boltzmann or Enskog equation, or by the Lattice Boltzmann method.

The shape of the dispersion relations for the growth rates also explains the suppression of instabilities through a reduction of the system size . Let $`k_0(N)=2\pi /L`$ be the smallest wave number, allowed in a system of linear extent $`L`$ at fixed density and with periodic boundary conditions. In systems with $`k_H^{}<k_0(N)<k_{}^{}`$, there are growing vorticity modes, but all density fluctuations are stable according to a linear stability analysis. In systems with $`k_0(N)<k_H^{}`$ the fluctuations in the density and in the flow field are unstable . However, systems with $`k_0(N)>k_{}^{}`$ do not show any instability. The smallest allowed wave numbers $`k_0(N)`$ for different numbers of particles, $`N=50,\mathrm{\hspace{0.17em}100},\mathrm{\hspace{0.17em}200},\mathrm{},\mathrm{\hspace{0.17em}1000}`$ at fixed packing fraction $`\varphi `$, are shown as stars in Fig.1. This information on small systems will be used in later sections.

We finally remark that several authors have also studied nonlinear terms in the macroscopic equations for granular fluids, such as the viscous heating term , and the nonlinear convective term . The inclusion of the combined effects of viscous heating and collisional cooling is essentially the new mechanism driving the dynamics of dissipative granular fluids in the time regime, directly following the homogeneous cooling state.

In an unpublished report  one of us has developed a systematic expansion method for solving the nonlinear hydrodynamic equations for granular fluids. The method is based on the separation of time scales of the relatively slow process of vorticity diffusion and the relatively fast process of heat conduction for all $`k`$-dependent excitations, allowed in the system ($`k>k_0(N)`$). This picture is consistent with the relative sizes of the growth rates $`z_{}(k)\gamma _0`$ and $`z_H(k)\gamma _0`$, as shown in Fig. 1.

The separation of the time scales permits a systematic approximation scheme, based on the Chapman-Enskog method or multi-time scale expansion method . In zeroth order this method yields a closed equation of motion for the rescaled flow field $`\stackrel{~}{𝐮}`$, as well as an equation for the global granular temperature $`\overline{T}(t)`$. In first approximation (see Ref. ) this method yields equations for the density and temperature fluctuations, $`\delta n`$ and $`\delta T`$, and a higher order correction to the Navier-Stokes equation.

In the present paper we will not explain the rather technical calculations based on the multi-time scale method, but we will try to elucidate the essential physics by deriving the zeroth order results in a more intuitive fashion. The systematic derivation will be published elsewhere.

## 3 Incompressible flows

As discussed in the previous section, instabilities and pattern formation occur in two different local fields, $`𝐮`$ and $`n`$, and on two different time scales, namely first patches of vorticity appear, and only much later density clusters appear. As explained above, the appearance of density clusters can even be further delayed, or all together suppressed by decreasing the system size . These observations suggest to analyze first the simplest case where fluctuations in density and temperature remain small. This would happen in the time regime following the unstable homogeneous cooling state, and possibly in the full time regime for sufficiently small systems, as suggested by the linear stability analysis. Of course the stability of solutions on the largest time scales is determined by the full nonlinear equations. The conditions for small $`n`$ and $`T`$fluctuations might be realized in incompressible flows, where $`𝐮=0`$. Then the density remains constant in the comoving frame, and the temperature balance equation greatly simplifies. As is well known from standard fluid dynamics and from the theory of turbulence , flows of elastic fluids are quite incompressible. This implies,

$$𝐮=0\mathrm{or}u_𝐤\widehat{𝐤}𝐮_𝐤=0,$$
(5)
where $`u_𝐤`$ is the longitudinal Fourier mode. Moreover, MD simulations and the theory of hydrodynamic fluctuations in granular flows  show that the incompressibility assumption, $`u_𝐤=0`$, remains valid down to very small wave numbers, satisfying $`k\xi _{}\stackrel{>}{}1`$ , and ultimately breaks down at the largest wavelengths, where $`\xi _{}1/\gamma _0`$ is the largest intrinsic dynamic correlation length in IHS fluid. It satisfies the inequality, $`\xi _{}\xi _{}(\nu /\omega \gamma _0)^{1/2}`$ for nearly elastic systems. Both correlation lengths are indicated in Fig.1 and defined more explicitly in .

Therefore, as a zeroth approximation to our nonlinear theory, we make the incompressibility assumption, $`𝐮=0`$, following Refs. , and the equations of motion become,

$`D_tn`$ $`=`$ $`0,`$
$`D_t𝐮`$ $`=`$ $`{\displaystyle \frac{1}{mn}}p+\nu ^2𝐮,`$
$`D_tT`$ $`=`$ $`b_T^2T+b_{}\left[𝐮+(𝐮)^{}\right]:𝐮2\gamma _0\omega T.`$ (6)
Consequently $`n(𝐫,t)`$ is constant in a comoving frame, and the set of nonlinear equations (6) can not describe the growth of inhomogeneities in the density field, but supposedly describe the time evolution of the system as long as the density fluctuations are small.

As a consequence of the incompressibility assumption the local density stays constant in time, and neither depressions, nor hot and cold regions can develop. Therefore, the pressure gradient in (6) remains negligibly small as well, and the Navier-Stokes equation becomes,

$$_t𝐮=𝐮𝐮+\nu ^2𝐮,$$
(7)
where the flow velocity $`𝐮=𝐮_{}`$ is purely rotational with $`𝐮=0`$. We also note that the divergence of Eq.(7), in combination with $`𝐮=0`$, reduces to $`(_\alpha u_\beta )(_\beta u_\alpha )=0`$ at all times. As already discussed in the introduction, the Burgers equation for a rotation-free flow field has been conjectured  as the large scale macroscopic equation for a d-dimensional dilute granular gas of arbitrary inelasticity $`ϵ`$.

Next we consider the temperature balance equation, which involves two processes: the diffusion process of heat conduction, where Fourier modes $`T_𝐤`$ decay with a rate $`b_Tk^2`$, and the global process of collisional cooling and nonlinear viscous heating, describing the decay of $`T_𝐤`$ for $`k0`$, or equivalently of $`\overline{T}(t)(1/V)_V𝑑𝐫T(𝐫,t)`$, referred to as global temperature.

To split off the behavior of $`\overline{T}(t)`$ from Eq.(6) we take the spatial average of the $`T`$–equation, yielding for the global temperature,

$$_t\overline{T}=b_{}\overline{|𝐮|^2}2\gamma _0\omega \overline{T},$$
(8)
where overlines denote spatial averages. We have introduced the notation $`|\text{A}|^2=_{\alpha \beta }|A_{\alpha \beta }|^2`$, for a second rank tensor A. The transport coefficient $`b_{}`$ and collision frequency $`\omega `$ are functions of the spatially average density $`\overline{n}=N/V`$ and temperature $`\overline{T}(t)`$ (see ’reference values’ below Eq. (3)). This is allowed as long as the local fluctuations $`\delta n=n\overline{n}`$ and $`\delta T=T\overline{T}`$ are not getting too large through nonlinear effects.

The new evolution equations (7) and (8) contain the time dependent coefficients $`b_{}`$, $`\omega `$ and $`\nu `$, which are proportional to $`\overline{v}_0(t)(2\overline{T}(t)/m)^{1/2}`$. Therefore, it is convenient to introduce the scaled field $`\stackrel{~}{𝐮}=𝐮/\overline{v}_0`$, and the scaled time $`\tau `$, defined as $`d\tau =\omega (\overline{T})dt`$. The final macroscopic evolution equations then become,

$`_\tau \mathrm{ln}\overline{T}`$ $`=`$ $`{\displaystyle \frac{4}{d}}𝒟_{}\overline{|\stackrel{~}{𝐮}|^2}2\gamma _0,`$
$`_\tau \stackrel{~}{𝐮}+l_0\stackrel{~}{𝐮}\stackrel{~}{𝐮}`$ $`=`$ $`𝒟_{}^2\stackrel{~}{𝐮}{\displaystyle \frac{1}{2}}(_\tau \mathrm{ln}\overline{T})\stackrel{~}{𝐮}`$ (9)
$`=`$ $`\gamma _0\stackrel{~}{𝐮}+𝒟_{}^2\stackrel{~}{𝐮}{\displaystyle \frac{2}{d}}𝒟_{}\overline{|\stackrel{~}{𝐮}|^2}\stackrel{~}{𝐮},`$
where $`l_0=\overline{v}_0/\omega `$ is the (time-independent) mean free path. The last equation is a closed equation for $`\stackrel{~}{𝐮}`$, and the physically consistent solutions need to obey the relation $`𝐮=0`$. The global temperature is slaved by the flow field. The rescaled vorticity diffusion coefficient, defined as $`𝒟_{}=\nu /\omega `$, is independent of time. The first term on the right hand side of the equation for $`\stackrel{~}{𝐮}`$ accounts for the instability, the second for the vorticity diffusion and the last one for the saturation effects, caused by the nonlinear viscous heating. It slows down the growth of unstable $`𝐤`$-modes, and may ultimately lead to a steady state for $`\stackrel{~}{𝐮}`$. For large times the above equations (7) and (8) are no longer valid, and nonlinear effects will induce density inhomogeneities, even in small systems with $`k_H^{}<k_0(N)<k_{}^{}`$.

## 4 Spontaneous symmetry breaking

In this section we will drop the nonlinear convective term $`\stackrel{~}{𝐮}\stackrel{~}{𝐮}`$, but at the end of our analysis we admit out of all possible solutions only those that satisfy Eq. (9) with the convective term included.

The final equation for the rescaled $`𝐮`$field has the form of the Landau-Ginzburg equation of motion for a non-conserved order parameter. This can be made more explicit by introducing the order parameter $`\text{S}=\stackrel{~}{𝐮}`$, and applying $``$ to the $`𝐮`$-equation in (9) with the result,

$`_\tau \text{S}`$ $`=`$ $`\left(\gamma _0+𝒟_{}^2\right)\text{S}{\displaystyle \frac{2}{d}}𝒟_{}\overline{|\text{S}|^2}\text{S}`$ (10)
$`=`$ $`V\delta [\text{S}]/\delta \text{S}^{},`$
where the last line contains the functional derivative of the energy functional $`[\text{S}]`$, defined as

$$[\text{S}]=\frac{1}{2}\gamma _0\overline{|\text{S}|^2}+\frac{1}{2}𝒟_{}\overline{|\text{S}|^2}+\frac{1}{2d}𝒟_{}\left(\overline{|\text{S}|^2}\right)^2.$$
(11)
In our further considerations it is more convenient to use Fourier modes. Moreover $`𝐮_𝐤=𝐮_𝐤`$, because of the assumption of incompressibility, and $`\text{S}_𝐤𝐤\stackrel{~}{𝐮}_𝐤`$. Then we obtain the equation of motion,

$`_\tau \text{S}_𝐤`$ $`=`$ $`V^2\delta [\text{S}]/\delta \text{S}_𝐤^{}`$ (12)
$`=`$ $`\left\{\gamma _0𝒟_{}k^2{\displaystyle \frac{2𝒟_{}}{dV^2}}{\displaystyle \underset{𝐪}{}}|\text{S}_𝐪|^2\right\}\text{S}_𝐤`$
with an energy functional,

$$[\text{S}]=\frac{1}{2V^2}\underset{𝐤}{}(\gamma _0+𝒟_{}k^2)|\text{S}_𝐤|^2+\frac{𝒟_{}}{2d}\left(\frac{1}{V^2}\underset{𝐤}{}|\text{S}_𝐤|^2\right)^2,$$
(13)
where the wave number $`𝐤=\mathrm{𝟎}`$ does not contribute.

The energy in (11) and (13) resembles a Landau free energy form for a tensorial order parameter $`\text{S}=\stackrel{~}{𝐮}`$ with tr $`\text{S}=\stackrel{~}{𝐮}=0`$, which is in fact the shear rate or rate of deformation tensor. It has a quartic term $`S^4`$, and a quadratic term $`S^2`$ with a coefficient that vanishes at a critical threshold $`k_{}^{}=(\gamma _0/𝒟_{})^{1/2}`$. It differs from the standard Landau free energy in that the quartic term contains summations over two totally independent wave numbers.

These results are very interesting. If the energy functional has a minimum, then there is a fixed point solution, $`\text{S}_𝐤(\mathrm{})`$, that is approached for large times. They are found by setting the right hand side of (12) equal to zero, i.e.,

$$\left\{\gamma _0𝒟_{}k^2\frac{2𝒟_{}}{dV^2}\underset{𝐪}{}|\text{S}_𝐪|^2\right\}\text{S}_𝐤=0.$$
(14)
We will show that depending on the parameter, $`k_0=2\pi /L`$, being above or below the threshold value $`k_{}^{}\sqrt{\gamma _0/D_{}}`$, the fixed point value of the order parameter, $`\{\text{S}_𝐤(\mathrm{})\}`$, is vanishing or non-vanishing. A stable fixed point with some non-vanishing Fourier components indicates that the system approaches a non-equilibrium steady state with a stationary pattern in the flow field, and spontaneous symmetry breaking has occurred.

Consider the right hand side of (12), and observe that the expression between curly brackets is necessarily negative for $`k>k_{}^{}`$, and the Fourier mode $`\text{S}_𝐤`$ decays to zero. If the smallest possible wave number satisfies, $`k_0>k_{}^{}`$, all $`\text{S}_𝐤`$ decay to zero, and there is no spontaneous symmetry breaking. The system remains spatially homogeneous. However, if $`k_0<k_{}^{}`$, then there exists the possibility that the expression inside brackets in Eq.(12) vanishes for a non-vanishing value of $`\text{S}_𝐤(\mathrm{})`$, i.e. there is an extremum determined by the condition,

$$\frac{1}{V^2}\underset{𝐪}{}|\text{S}_𝐪|^2=\frac{d}{2}(\gamma _0𝒟_{}k^2)/𝒟_{}.$$
(15)
If the extremum is a saddle point then there are directions of exponentially growing solutions. However, if $`kk_{}^{}`$, then the right hand side of (15) may become zero and even positive. There is the possibility of stationary and of exponentially growing solutions. The fixed point $`\{\text{S}_𝐤(\mathrm{})=0\text{ for any }𝐤\text{ }\}`$ is a saddle point. One can also show  that all fixed point solutions with non-vanishing $`\text{S}_𝐤`$ for any $`|𝐤|k_0`$ are saddle points, and that the only non-vanishing solution $`\{\text{S}_𝐤0\text{ if }|𝐤|=k_0,\text{S}_𝐤=0\text{ otherwise}\}`$ is a stable fixed point with an infinitely degenerate minimum, given symbolically by the Mexican hat shape as illustrated in Fig. 3. The condition (15) for $`𝐮_𝐤`$ with $`|𝐤|=k_0=2\pi /L`$ in $`d`$-dimensions is then

$$\frac{1}{V^2}\underset{\alpha =1}{\overset{d}{}}\left|\stackrel{~}{𝐮}_{𝐤_{0\alpha }}\right|^2=\frac{d}{4}(\gamma _0𝒟_{}k_0^2)/𝒟_{}k_0^2,$$
(16)
where $`𝐤_{0\alpha }=k_0\widehat{𝐞}_\alpha `$ and $`\widehat{𝐞}_\alpha `$ is a unit vector in the direction $`\alpha `$ ($`\alpha =\{1,2,\mathrm{}d\}`$). As the solutions of these equations have to satisfy Eq. (5), the Fourier components can be written as

$`\stackrel{~}{𝐮}_{𝐤_{0\alpha }}={\displaystyle \frac{1}{2}}V{\displaystyle \underset{\beta (\alpha )}{}}\widehat{𝐞}_\beta A_{\alpha \beta }e^{i\theta _{\alpha \beta }},`$ (17)
where the phases $`\theta _{\alpha \beta }`$ and amplitudes $`A_{\alpha \beta }`$ are $`2d(d1)`$ real numbers. On account of (16) the amplitudes satisfy the relations

$`A_0^2\underset{1\alpha \beta d}{{\displaystyle }}A_{\alpha \beta }^2=d(\gamma _0𝒟_{}k_0^2)/𝒟_{}k_0^2.`$ (18)
The stable stationary solutions in real space are then

$`\stackrel{~}{𝐮}_0(𝐫)=\underset{\alpha \beta }{{\displaystyle }}A_{\alpha \beta }\widehat{𝐞}_\beta \mathrm{cos}\left(k_0r_\alpha +\theta _{\alpha \beta }\right).`$ (19)
Out of this set of solutions we select those that satisfy the full nonlinear Eq. (9) with the convective term included, i.e. we determine the $`d(d1)`$ amplitudes $`A_{\alpha \beta }`$ such that the relation, $`\stackrel{~}{𝐮}_0\stackrel{~}{𝐮}_0=\mathrm{𝟎}`$, is satisfied. This yields $`d`$ conditions. For the two-dimensional case only two fixed point solutions remain, instead of infinitely many degenerate minima, i.e.

$`\stackrel{~}{𝐮}_0(x)`$ $`=`$ $`\widehat{𝐞}_yA_0\mathrm{cos}(k_0x+\theta _x)`$
$`\stackrel{~}{𝐮}_0(y)`$ $`=`$ $`\widehat{𝐞}_xA_0\mathrm{cos}(k_0y+\theta _y).`$ (20)
The symmetry of the steady state is spontaneously broken, and spontaneous fluctuations in the initial state determine which of these two minima will be reached.

In summary, the equation, describing the growth dynamics of vortices in granular fluids, without the convective term is described by a time dependent Landau-Ginzburg equation for a non-conserved order parameter, $`\text{S}=\stackrel{~}{𝐮}`$, derived from an energy functional with a continuous set of degenerate minima, having the shape of a Mexican hat. The last observation may suggest that unstable granular fluids have some resemblance to Model H  which has a similar energy surface in the neighborhood of its stable fix points. However, this is not the case. Addition of the nonlinear convective term to Eq.(10) selects out of this infinite set of minima only subsets of admissible solutions. In two dimensions only two distinct minima survive. Therefore, the unstable two-dimensional IHS fluid has a greater resemblance to spinodal decomposition for a non-conserved order parameter, which belongs to the universality class of Model A .

It should be noted that the complete solution of Eq. (9) with the convective term included may have a larger set of physically acceptable solutions. However, we have not been able to find any.

## 5  MD-Simulations

In this section we present results obtained from the zeroth approximation and compare the theoretical predictions with the results of computer simulations of small systems in two dimensions. At the end of this section, some results of the first order approximation in Ref.  are presented without derivation and compared with computer simulations.

A snapshot of a typical configuration for a small system with $`k_H^{}<k_0<k_{}^{}`$ at large times $`\tau \tau _{cr}`$ is shown in Fig. 3. A shear flow with variation of the $`u_x`$component in the $`y`$direction is observed. Individual $`v_x`$components of the particle velocities are plotted in Fig. 4 versus their $`y`$coordinate. A fit to a sinusoidal curve (solid line) shows that the solution described in (20) is realized.

We assume that the number of collisions $`\tau `$ in the simulations of Figs. 3 and 4 is sufficiently large, so that the components $`\stackrel{~}{𝐮}_𝐤`$ are very close to their fixed point values. Then, the relation, $`𝒟_{}\overline{|\stackrel{~}{𝐮}|^2}=\frac{1}{2}d(\gamma _0𝒟_{}k_0^2)`$, follows from (16) and (9), and $`_\tau \mathrm{ln}\overline{T}2𝒟_{}k_0^2`$. The global temperature for large times, i.e. $`\tau \tau _{cr}=L^2/𝒟_{}`$, becomes,

$$\overline{T}T_e\mathrm{exp}(2𝒟_{}k_0^2\tau ).$$
(21)
Equation (21) gives the temperature as a function of $`\tau `$. The integration constant cannot be determined from the theory, but a fit of (21) to the simulation data in Fig. 5 for $`\tau >200`$ yields $`T_e3.87\times 10^4`$. We also note that the mode coupling theory of Ref.  gives the same decay rate for the total energy, when this theory is applied to a small system. In that case the Fourier sum $`(1/V)_𝐪`$ can not be replaced by $`𝑑𝐪/(2\pi )^d`$. On the contrary, the sum is essentially given by its first term only.

The relation between $`\tau `$ and $`t`$ is defined through $`d\tau =\omega (\overline{T})dt`$, where the collision frequency $`\omega `$ is proportional to the square root of $`\overline{T}`$, i.e.

$$\frac{d\tau }{dt}=\omega (\overline{T})=\omega _0\sqrt{\frac{\overline{T}}{T_0}}\omega _0\sqrt{\frac{T_e}{T_0}}\mathrm{exp}\left[𝒟_{}k_0^2\tau \right],$$
(22)
where $`\omega _0=\omega (T_0)`$ is the collision frequency at the initial time. As $`\omega `$ can be measured directly as a function of $`t`$ in MD simulations, it yields an independent determination of $`T_e`$ with the result $`T_e3.68\times 10^4`$. Integration of (22) yields,

$$\mathrm{exp}(𝒟_{}k_0^2\tau )\omega _0\sqrt{\frac{T_e}{T_0}}𝒟_{}k_0^2(tt_e),$$
(23)
valid for $`(tt_e)`$ large. After eliminating $`\tau `$ from Eq. (21) and (23) we obtain the global temperature for asymptotically large time $`t`$ as,

$$\overline{T}\frac{T_0}{(\omega _0𝒟_{}k_0^2)^2}\frac{1}{(tt_e)^2},$$
(24)
where visual inspection of Fig. 6 shows that $`t_e10^4`$. The temperature shows algebraic decay with the same exponent as in Haff’s law (4), but the prefactor in Haff’s law does neither depend on the system size, nor on the dimensionality. In (24) it is proportional to $`L^4N^{4/d}`$, whereas the prefactor in Haff’s law is independent of the system size (see Fig.6).

Once we have $`\overline{T}`$, we can calculate the energy per particle, defined in  as,

$$E=\frac{1}{n}\left[\overline{\frac{m}{2}n𝐮^2+\frac{d}{2}nT}\right]=\overline{T}\left[\overline{\stackrel{~}{𝐮}^2}+\frac{d}{2}\right],$$
(25)
where local density variations are supposedly small and the relation $`𝐮=\overline{v}_0\stackrel{~}{𝐮}`$ has been used. Combination of this result with (16) yields,

$$E\frac{d\gamma _0}{2𝒟_{}k_0^2}\overline{T}.$$
(26)
This shows that at long times, the energy and temperature are proportional. Therefore, energy decays exponentially when expressed in terms of $`\tau `$ (see (21)), with the same decay rate as the global temperature. This property is observed by MD simulations, as shown in Fig. 5. From (26) and (24) the decay of the energy per particle in terms of $`t`$ in the long time limit is given by

$$E(t)\frac{d\gamma _0T_0}{2\omega _0^2(𝒟_{}k_0^2)^3}\frac{1}{(tt_e)^2}B\frac{N^{6/d}}{(tt_e)^2}.$$
(27)
Hence, if we vary the number of particles $`N`$, at fixed packing fraction $`\varphi `$, the amplitude of $`t^2`$-decay of the energy is proportional to $`N^3`$ in two-dimensional systems, or to $`N^{6/d}`$ in $`d`$-dimensional systems. Comparison of Eq. (27) with simulations is shown in Fig. 6. It would also be interesting to compare the coefficient $`B`$ in (27) with the simulation results of Ref. for the IHS fluid which have been performed in five and six dimensions.

Finally, a systematic expansion to be presented in , shows that the stationary solution of local temperature and local density are given by

$`\delta \stackrel{~}{T}{\displaystyle \frac{T(y)}{\overline{T}}}1`$ $``$ $`{\displaystyle \frac{b_TV_0^2}{2d\omega }}\mathrm{cos}[2(k_0y+\theta _y)]`$
$`\delta \stackrel{~}{n}{\displaystyle \frac{n(y)}{\overline{n}}}1`$ $``$ $`{\displaystyle \frac{b_TV_0^2\overline{T}\left(\frac{p}{T}\right)}{2d\omega \overline{n}\left(\frac{p}{n}\right)}}\mathrm{cos}[2(k_0y+\theta _y)],`$ (28)
where $`\overline{n}=N/V`$, and $`\theta _y`$ is the same phase factor as given in (20). It means that the density and temperature inhomogeneities show a period which is half the period of the shear flow profile. The relation between the spatial periods has already been observed in Ref. in MD simulations of a fluid of two-dimensional hard disks, and in Ref.  in Direct Monte Carlo simulations of the Boltzmann equation for an IHS gas. Figure 7 shows the observed density and temperature inhomogeneities. A fit to a sinusoidal curve supports the temperature and density profiles given by (28).

As shown in Eq.(28), the temperature and density inhomogeneities are in opposite phase, implying that dense regions are cold, and the dilute regions hot. The amplitudes are such that the overall pressure is constant, as we have assumed in the course of the paper.

## 6 Conclusions

Under the restrictions imposed in our derivations, as discussed in Sec.3, we have shown that the unstable dynamics and formation of vortex patterns in the flow field of a freely evolving fluid of inelastic hard spheres (IHS) can be cast in the form of a time-dependent Landau-Ginzburg-type model for a non-conserved order parameter. In the two-dimensional case (but not for $`d3`$) the growth of the vortex pattern in the IHS fluid is qualitatively similar to spinodal decomposition for a non-conserved scalar order parameter, referred to as model A in the Hohenberg–Halperin classification . The analogy between unstable IHS fluid in two-dimensions and spinodal decomposition has already been pointed out in Ref. . A difference between our model and model A is that the energy functional (11) contains a quartic term with summations over two independent wave numbers. This implies a non-local interaction of the order parameter $`𝖲=\stackrel{~}{𝐮}`$. The non-local interaction is caused by the scaled field $`\stackrel{~}{𝐮}=𝐮/\overline{v}_0𝐮/\sqrt{\overline{T}}`$. Because the global temperature $`\overline{T}`$ is determined by $`𝖲`$ in all space (see Eq.(8)), the evolution of a local order parameter $`𝖲`$ is affected by $`𝖲`$ at any other points in space through $`\overline{T}`$.

We have shown that nonlinear viscous heating gives rise to a quartic term in the energy functional, which is responsible for saturation of unstable vorticity modes. Unstable vorticity modes initially grow as predicted by the linear stability analysis, but eventually saturate because of nonlinear viscous heating. The influence of nonlinear viscous heating on the formation of temperature inhomogeneities and clustering has been pointed out by Goldhirsch and Zanetti , and investigated in more detail by Brey et al. , by comparing the results of a hydrodynamic model with viscous heating, using direct Monte Carlo simulation of the Boltzmann equation.

The theoretical predictions of our zeroth order theory for the flow field and the decay of the energy are in quantitative agreement with MD simulations of small systems as shown in Sec.5. They support the intuitive arguments used in deriving the zeroth order description, presented in Sec.3. Moreover, the results presented here can be obtained as the lowest order approximation in a multi-time scale expansion . This theory makes no assumption of incompressibility. There are some properties of those small systems that can not be described by the zeroth order theory. Density and temperature inhomogeneities are only found in the next order of approximation.

It is worth mentioning that results of our zeroth and first order approximations are consistent with the results of a nonlinear analysis by Soto et al.  of systems which are close to the stability thresholds $`k_{}^{}`$. If the smallest wave number of the system $`k_0`$ is slightly smaller than $`k_{}^{}`$, vorticity modes with the wave number $`k_0`$ grow so slowly that the rest of hydrodynamic modes are enslaved by these vorticity modes. Besides, the amplitudes of these vorticity modes at long times are expected to saturate because of nonlinear effects and remain small. On the basis of this consideration, they obtained amplitude equations for the vorticity modes with $`k_0`$. The inhomogeneous density and temperature are given as functions of the vorticity modes with $`k_0`$.

Finally, we emphasize that the periodic boundary conditions, with which the simulations were carried out, obviously play a crucial role in the formation of a shear flow profile and of the density and temperature profiles, presented in (20) and (28) respectively. Indeed, once the typical size of the vortex patterns becomes comparable to the length $`L`$ of the system for $`\tau >\tau _{cr}=L^2/𝒟_{}`$, the artificial periodic boundary conditions start to affect the evolution of the system. Consequently, the long time regime $`\tau \tau _{cr}`$ described by the stationary solution is of less physical interest than the regime of unstable growth $`\tau \tau _{cr}`$, where the typical size of the vortex patterns remains small compared to the length $`L`$ of the system. A thermodynamically large system is always in the unstable growth regime. Unfortunately, the only analytic ’large’ time results for the unstable growth regime of granular fluids in thermodynamically large systems, have been obtained from fluctuating (linear) hydrodynamic equations or from mode coupling theory , but not from truly nonlinear theories. An exception is a dilute gas of inelastic point particles in one dimension , where strong evidence supports the conjecture that the large space-time behavior is described by the adhesion model and the Burgers equation.

## 7 Acknowledgements

M.E. acknowledges stimulating discussion with R. Desai and R. Kapral. J.W. and R.B. acknowledge support of the foundation ”Fundamenteel Onderzoek der Materie (FOM)”, which is financially supported by the Dutch National Science Foundation (NWO). J.W. also acknowledges support of a Huygens scholarship. R.B. wants to thank the Institute for Theoretical Physics of Universiteit Utrecht for its hospitality. R.B. is supported by grant DGES-PB97-0076 (Spain).


text:
# Zone-boundary excitations in coupled Haldane spin chain compounds PbNi2V2O8 and SrNi2V2O8.

## I Introduction

The quasi-one-dimensional (quasi-1D) $`S=1`$ antiferromagnet (AF) PbNi2V2O8 was recently shown to be the first example of a Haldane-gap system to undergo long-range magnetic ordering upon doping with non-magnetic impurities. This unique behavior was understood through inelastic neutron scattering measurements, that suggested that three dimensional (3D) inter-chain interactions in PbNi2V2O8 are almost (but not quite) strong enough to destroy the Haldane singlet ground state even in the undoped material. In fact, in the very similar isostructural compound SrNi2V2O8, inter-chain coupling exceeds the critical value, and, unlike the Pb-based material, the undoped Sr-Ni vanadate orders magnetically at $`T=7`$ K.

The introduction of non-magnetic vacancies into gapped quantum spin chains produces free $`S=1/2`$ degrees of freedom on either side of each impurity site. In the presence of arbitrarily weak inter-chain interactions, these free spins will order in three dimensions at sufficiently low temperature. However, next-nearest neighbor (nnn) interactions within the chains can prevent this from happening. Naively, in the presence of AF nnn coupling, one can expect the liberated spins around each impurity to bind into non-magnetic singlets (dimers), preserving the disordered ground state of the system.

There are at least two indications to that nnn interactions may be active in PbNi2V2O8 and SrNi2V2O8. First, the spin chains in both materials are not linear, but are formed by Ni2+ ions arranged in 4-step spirals propagating along the $`c`$ axis of the orthorhombic lattice (for a detailed discussion of chain geometry and inter-chain interactions see Ref. ). In this “twisted” structure nnn superexchange pathways involving equatorial oxygens of the NiO6 tetrahedra can be quite effective. Second, there remains some discrepancy between the in-chain exchange constant $`J`$ deduced from the Haldane gaps measured with neutron scattering, assuming only nearest-neighbor (nn) interactions ($`J=9.0`$ meV), and from high-temperature ($`TJ`$) susceptibility data ($`J=8.2`$ meV). The former technique probes low-energy excitations around $`q_{}=\pi `$ (AF zone-center), while the latter is sensitive to excitations in the entire Brillouin zone. The mismatch may be interpreted in terms of an additional coupling with a characteristic interaction length of twice the spin-spin separation in the chains, that is known to have different effects on zone-center and zone-boundary excitation energies.

The most direct way to resolve this issue would be to measure the dispersion of magnetic excitations in the entire Brillouin zone in a single crystal sample. Unfortunately, no single-crystal or aligned powder samples of either PbNi2V2O8 or SrNi2V2O8, suitable for inelastic neutron scattering studies, are available to date. As will be discussed below, some information on nnn coupling can be obtained using powder samples, if the measurements are extended all the way to the zone-boundary energy. Previous neutron scattering studies employed cold neutrons to investigate low-energy excitations in the two compounds. In the present paper we report thermal-neutron studies specifically aimed at measuring high-energy excitations. We then engage in a semi-quantitative discussion of the role nnn interactions play in vacancy-induced ordering of weakly-coupled Haldane spin chains, and apply the results to the particular case of PbNi2V2O8.

## II Experimental procedures and results

Both vanadates crystallize in tetragonal structure, space group $`I41cd`$, with lattice constants $`a=12.249(3)`$ Å, $`c=8.354(2)`$ Å for PbNi2V2O8, and $`a=12.1617`$ Å, $`c=8.1617`$ Å for SrNi2V2O8, respectively. All measurements were performed on large ($`5`$ g) powder samples of PbNi2V2O8 and SrNi2V2O8, similar to those used in previous cold-neutron studies. The experiments were carried out at the HB3 3-axis spectrometer installed at the High Flux Isotope Reactor at Oak Ridge National Laboratory. Neutrons with a fixed final energy of 13.5 meV were used with a 5 cm thick Pyrolitic Graphite (PG) filter positioned after the sample to eliminate higher-order beam contamination. PG(002) reflections were used for monochromator and analyzer, in combination with $`48^{}40^{}40^{}120^{}`$ collimators. Sample environment was in all cases a standard He-4-flow cryostat. The data were taken at $`T=1.5`$ K.

Typical data sets collected in constant-$`Q`$ scans for PbNi2V2O8 and SrNi2V2O8 are shown in Fig. 1 and Fig. 2, respectively. Figure 3 shows the bulk of the data collected for the Sr system (only 5 scans, at $`|Q|/c^{}=`$ 2.5, 3, 3.5, 4 and 5, were collected for the Pb-based material). The measured intensities were scaled to correct for the previously determined energy-dependent $`\lambda /2`$ contribution to the incident beam, that affects the monitor rate. Within the incident neutron energy range used, the monitor efficiency changes by roughly 30%. In addition, a flat background (17 counts/min), measured at 30 meV energy transfer (above the magnetic excitation band), was subtracted from the data.

For both materials all scans have the characteristic M-shape. In a powder sample one measures the spherical average of the dynamic structure factor $`S(𝑸,\omega )`$. For large momentum transfers, a large portion of $`𝑸`$ space is sampled simultaneously, and the powder cross section is in many ways similar to the density of states function $`n(\omega )`$. In particular, one expects to see a peak in the powder cross section wherever the density of states is large, i.e., at the gap energy and at the zone boundary. As discussed previously, the peak around 5 meV energy transfer corresponds to excitations close to the Haldane gap at a momentum transfer $`q_{}=\pi +2\pi n`$ ($`n`$-integer) along the chain axis. The weaker higher-energy peaks, represent excitation near the top of the band, i.e., in the vicinity of $`q_{}=\pi /2+\pi n`$ ($`n`$-integer).

## III Data analysis

### A Model cross section

The convoluted nature of inelastic powder data make it notoriously difficult to interpret (for an example see Ref. , particularly Section IV). The only way to extract information on the dynamic structure factor from such measurements is to assume a “reasonable” model and refine the parameter values to best-fit the data. The main assumption of the model used for PbNi2V2O8 and SrNi2V2O8 to date is that these compounds are described as weakly-interacting $`S=1`$ AF chains with magnetic anisotropy. This picture is based on susceptibility, high field and specific heat measurements, and is supported by the observed unique shape of constant-$`E`$ neutron scans. In the present work we shall adopt this model as well.

The form of the neutron scattering cross section for this model was derived in Ref. , and is written in the single mode approximation (SMA). Inter-chain interactions were accounted for within the Random Phase Approximation (RPA). This cross section function performed rather well in describing low-energy neutron data. Previous studies of other Haldane-gap compounds, such as NENP and CsNiCl3 , have demonstrated that the SMA is, in fact, applicable in most of the Brillouin zone. The previously derived SMA cross section (Eqs. 5–13 of Ref. ) could thus be safely used in interpreting our high-energy results. However, the expression for the dispersion relation (Eqs. 14 and 15 of Ref. ) had to be slightly modified. These modifications are discussed in the following paragraphs.

#### 1 Isolated chains with nearest-neighbor interactions

In Ref.  the dispersion of excitations in non-interacting Haldane spin chains was assumed to be sinusoidal and symmetric around the classical AF zone-boundary $`q_{}=\pi /2`$. In fact, since translational symmetry is not broken in the Haldane ground state, the excitation energy extrapolates to a larger value at $`q_{}0`$ than at $`q_{}\pi `$. This effect can be taken into account empirically by including an extra term in the dispersion relation. The dynamic susceptibility of a single Haldane spin chain is then written as:

$`\chi ^{\beta \beta }(q_{},\omega )=\chi _\pi ^{\beta \beta }{\displaystyle \frac{1\mathrm{cos}q_{}}{2}}{\displaystyle \frac{\mathrm{\Delta }_{\beta \beta }^2}{(\mathrm{}\omega _q_{}^{\beta \beta })^2(\mathrm{}\omega +iϵ)^2}},`$ (1)
$`(\mathrm{}\omega _q_{}^{\beta \beta })^2=\mathrm{\Delta }_{\beta \beta }^2+v^2\mathrm{sin}^2q_{}+\alpha ^2\mathrm{cos}^2{\displaystyle \frac{q_{}}{2}}.`$ (2)
In these formulas, which replace Eq. 22 in Ref. , $`\mathrm{\Delta }_{\beta \beta }`$ is the gap energy for the Haldane gap mode of a particular polarization $`\beta \beta `$, $`\chi _\pi ^{\beta \beta }`$ is the static staggered susceptibility, and $`v`$ is the spin wave velocity. The spin chains PbNi2V2O8 and SrNi2V2O8 contain four equivalent spins per crystallographic period along the $`c`$ axis. Throughout this work the actual wave vector transfer will be denoted as $`𝑸=h𝒂^{}+k𝒃^{}+l𝒄^{}`$. The symbol $`q_{}𝑸𝒄/4`$ in Eq. 2 will stand for the reduced wave vector transfer along the chains.

The energetics of a Haldane-gap AF is by now very well established through numerical work by several authors. For an $`S=1`$ Heisenberg chain with single-ion anisotropy of type $`D(S^z)^2`$ the approximate values of parameters are as follows:

$`\mathrm{\Delta }_{xx}`$ $`=`$ $`\mathrm{\Delta }_{yy}\mathrm{\Delta }{\displaystyle \frac{2}{3}}D,`$ (3)
$`\mathrm{\Delta }_{zz}`$ $``$ $`\mathrm{\Delta }+{\displaystyle \frac{4}{3}}D,`$ (4)
$`\mathrm{\Delta }`$ $``$ $`0.41J,`$ (5)
$`v`$ $``$ $`2.49J,`$ (6)
$`\chi _\pi ^{\beta \beta }`$ $``$ $`1.26{\displaystyle \frac{v}{\mathrm{\Delta }_{\beta \beta }^2}}.`$ (7)
The parameter $`\alpha `$ in Eq. 2 characterizes the asymmetry of the dispersion relation. Existing experimental data for NENP and CsNiCl3  gives $`\alpha =1.45J0.58v`$ and $`\alpha =1.1J0.44v`$, respectively. The latter results inspires more confidence, as it was obtained with a much higher experimental resolution and gives a dispersion relation very similar to that obtained by direct diagonalization in finite-length chains. Throughout the rest of the paper we shall assume $`\alpha `$ fixed at $`\alpha 0.44v`$.

#### 2 Next-nearest-neighbor interactions

Next-nearest-neighbor in-chain interactions, that we believe can be relevant for PbNi2V2O8 and SrNi2V2O8, can be accounted for within the RPA. Isotropic nnn coupling of magnitude $`J^{}`$ gives an additional term in the r.h.s. of Eq. 2 that can be written as:

$$\chi _\pi ^{\beta \beta }\mathrm{\Delta }_{\beta \beta }^2J^{}\frac{1\mathrm{cos}q_{}}{2}\mathrm{cos}2q_{}.$$
(8)
The primary effect of nnn interactions is to modify the gap energies:

$`\mathrm{\Delta }_{\beta \beta }^2=\mathrm{\Delta }_{\beta \beta ,0}^2\left(1+\chi _\pi ^{\beta \beta }J^{}\right),`$ (9)
where $`\mathrm{\Delta }_{\beta \beta ,0}`$ are the gaps in the absence of nnn coupling. The relative effect of nn interactions on the zone-boundary energy and spin wave velocity is considerably smaller, and can be ignored. To a good approximation, in analyzing the experimental data one can use the dispersion relation 2, but treat the gap energies and $`v`$ as independent parameters, rather than assuming the relation 6 valid only in the case $`J^{}=0`$. The ratio $`\xi =v/\mathrm{\Delta }`$ can then be determined and used to obtain an estimate for $`J^{}`$. From Eqs. 37,9 for $`J^{}J`$ one gets:

$$\xi \xi _0(19.33\frac{J^{}}{J})$$
(10)
Here $`\xi _0=v/\mathrm{\Delta }_06.07`$ is equal to the equal-time correlation length in the absence of nnn coupling.

\[

\]

#### 3 Inter-chain interactions

For the coupling geometry found in PbNi2V2O8 and SrNi2V2O8, switching on inter-chain interactions within the RPA gives an additional term

$$\chi _\pi ^{\beta \beta }\mathrm{\Delta }_{\beta \beta }^2J_1\mathrm{cos}(\pi l/2)\left[\mathrm{cos}(\pi h)+\mathrm{cos}(\pi k)\right]\frac{(1\mathrm{cos}\pi l/2)}{2}$$
(11)
in the r.h.s. of Eq. 2. In this formula, which corresponds to the last terms in Eqs. 14 and 15 of Ref. , $`J_1`$ is the magnitude of nearest-neighbor inter-chain interactions. For our present purposes, the possible small anisotropy of $`J_1`$ will be ignored.

### B Parameters for PbNi2V2O8 and SrNi2V2O8

Our understanding of the physics of PbNi2V2O8 and SrNi2V2O8 is based on the knowledge of the following characteristics: the intrinsic gap energies $`\mathrm{\Delta }_{zz}`$ and $`\mathrm{\Delta }_{xx}=\mathrm{\Delta }_{yy}`$, the magnitude of inter-chain interactions $`J_1`$ and the spin wave velocity $`v`$. Can these four parameters be unambiguously determined from the available experimental data?

First, consider those quantities that characterize low-energy excitations near $`q_{}=\pi `$. For PbNi2V2O8 the actual minima of the 3D dispersion, $`E_{\mathrm{min},\beta \beta }`$ were independently measured in high-field magnetization studies of aligned powder samples. These energies are directly related to $`\mathrm{\Delta }_{\beta \beta }`$ and $`J_1`$. The third experimental value need to unambiguously “untangle” these three quantities was probed in previous high-resolution cold neutron measurements. Indeed, the well-defined lower-energy 5 meV peak in constant-$`q`$ powder scans is located approximately at the top of the dispersion band perpendicular to the chains at $`q_{}=\pi `$, and its position is directly related to $`\mathrm{\Delta }_{\beta \beta }`$ and $`J_1`$. From the combined high-field and neutron measurements $`\mathrm{\Delta }_{\beta \beta }`$ and $`J_1`$ can thus be reliably determined within the established model. For PbNi2V2O8 these parameters are: $`\mathrm{\Delta }_{zz}=3.1(3)`$ meV, $`\mathrm{\Delta }_{xx}=\mathrm{\Delta }_{yy}=4.0(3)`$ meV and $`J_1=0.17(2)`$ meV.

For SrNi2V2O8, which is ordered at low temperature, the model of weakly-coupled Haldane spin chains is, strictly speaking, not applicable. Nevertheless, as the ordered moment is expected to be very small, the model is still a good approximation, except in the direct proximity of the magnetic Bragg peaks. From the analysis of the low-energy powder spectrum for SrNi2V2O8 the following parameters were previously determined: $`\mathrm{\Delta }_{zz}=2.8(4)`$ meV, $`\mathrm{\Delta }_{xx}=\mathrm{\Delta }_{yy}=3.9(3)`$ meV and $`J_1=0.17(3)`$ meV.

The exact value of the remaining “high-energy” parameter $`v`$ has little impact on the low-energy powder spectrum, and was not treated as an independent variable in Ref. . In contrast, $`v`$ determines the high-temperature magnetic susceptibility of the system, that, provided $`k_\mathrm{B}T\mathrm{\Delta }`$, are insensitive to the details of the low-energy spectrum. This type of experiments yields $`v20.4`$ meV for PbNi2V2O8. The present study provides the most direct measure of the zone-boundary energy, and hence the spin wave velocity $`v`$, defined by the location of the higher-energy peak and cutoff in the constant-$`q`$ powder scans. In analyzing these data we shall use the cross section defined in Eqs. 5–13 of Ref.  in combination with the dispersion relation given by Eqs. 2 and 11 above, using only $`v`$ as an adjustable parameter of the model.

### C Fits to experimental data.

In the data analysis, the powder-average of the parameterized cross section was calculated numerically using a Monte Carlo routine. The result was further numerically convoluted with the $`QE`$ resolution function of the 3-axis spectrometer, calculated in the Cooper-Nathans approximation. Performing a true global fit to the data, using a single set of parameters to describe all measured scans, was not possible because of the scattering-angle dependence of neutron transmission in the sample. The effect is due to a large sample size and appreciable incoherent cross sections of Ni and V nuclei. Transmission in the sample was studied by measuring the intensity of incoherent elastic scattering, and found to change by $`30`$%, depending on scattering geometry. Within each constant-$`Q`$ scan the correction is small and could be ignored. However, in analyzing the entire data sets, collected in a broad range of momentum transfers (scattering angles), a separate intensity prefactor had to be used for each const-$`Q`$ scan. Thus, for each compound, the data were analyzed in a semi-global fit: a single adjustable parameter $`v`$ was used to simultaneously describe all measured scans, each with its own adjustable prefactor. The parameters were refined to best-fit the data using a least-squares algorithm.

Excellent fits were obtained with $`v=20.6(7)`$ meV for PbNi2V2O8 and $`v=21.2(4)`$ meV for SrNi2V2O8. Simulations based on the refined parameter values are shown in solid lines in Figs. 13. The fits correspond to $`\chi ^2=4.7`$ for PbNi2V2O8 and $`\chi ^2=3.8`$ for SrNi2V2O8. Given that the overall shape of the measured scans is reproduced by our model quite well, the large values of $`\chi ^2`$ are most likely due to the fact that experimental statistical error bars are less significant than unavoidable systematic errors, such as multi-phonon scattering in the sample. While such effects can lead to an inaccurate background subtraction, they are not likely to influence the measured zone-boundary energy, as the latter is accurately pin-pointed by the well-defined higher-energy cutoff in the energy scans. The relevant parameters used to fit the data, both fixed and adjustable, are summarized in Table LABEL:table for both systems. The refined values for $`v`$ are to be compared to previous estimates based on the measured gap energies and Eq.5, assuming only nearest-neighbor interactions, $`v=22.5(1.0)`$ meV and $`v=21.5(0.9)`$ meV, for PbNi2V2O8 and SrNi2V2O8, respectively. For comparison, the dashed lines in Figs. 1 are simulations based on these latter value for PbNi2V2O8.

It was important to verify that any uncertainty that may exist for the dispersion asymmetry parameter $`\alpha `$ does not have a significant impact on our determination of $`v`$. This was done by analyzing the PbNi2V2O8 data using a larger value $`\alpha =0.58v`$, as previously experimentally determined for NENP, rather than $`\alpha =0.44v`$. Essentially the same quality of fits was obtained, and the refined spin wave velocity $`v=19.9(8)`$ meV is within the error bar of that obtained with $`\alpha =0.44v`$, as described above. This demonstrates that the data analysis procedure employed is fairly insensitive to dispersion asymmetry.

We see that any deviations from the model involving only nearest-neighbor interactions are small in either material. While for SrNi2V2O8, within the error bars, the discrepancy is undetectable, for PbNi2V2O8 it appears to be more substantial. It is gratifying to see that for the latter system, the bandwidth measured in this work is in excellent agreement with estimates based on the high-temperature susceptibility (see above). As explained in the previous section, the reduction of zone-boundary energy is a sign of antiferromagnetic nnn in-chain interactions. Using Eq. 10 for PbNi2V2O8 we obtain $`J^{}0.4(2)`$ meV. Detecting such small values of $`J^{}`$ is only possible due to the large staggered susceptibility of a Haldane spin chain (the gap is relatively small compared to $`v`$), and hence the large prefactor in Eq. 10.

## IV Discussion

To summarize, the previously proposed nearest-neighbor model for PbNi2V2O8 and SrNi2V2O8 should be an adequate description of these compounds in most cases, in the entire Brillouin zone. For PbNi2V2O8 the data suggest, barely outside the error bars, a small antiferromagnetic nnn in-chain coupling of a few tenths of a meV. But how strong should this coupling be to have an impact on long-range ordering in a doped material? We shall now discuss the general significance of nnn interactions to spin-vacancy-induced order in weakly-coupled integer-spin quantum chains, and then consider the particular case of PbNi2V2O8.

### A Interactions between end-chain spins

The problem of spin-vacancy-induced ordering in coupled integer-spin quantum chains with only nearest-neighbor in-chain interactions has been studied theoretically by Shender and Kivelson. The introduction of a non-magnetic impurity in place of a $`S=1`$ ion produces two free $`S=1/2`$ spins located on the two ends of the severed chain. In each of the two fragments the end-chain spin is delocalized on a length scale defined by the intrinsic equal-time correlation length $`\xi 6`$ of a Haldane spin chain. In principle, existing 3D interactions in the system should lead to long-range magnetic ordering of these liberated end-chain spins at a non-zero temperature.

The most interesting result of Ref.  is that if inter-chain interactions are below the threshold value at which even the undoped system orders in three dimensions, there must be a critical spin-vacancy concentration $`x_c`$ below which the system remains disordered at $`T=0`$. The mechanism of this behavior is illustrated in Fig. 4, which show spins in the vicinity of one non-magnetic impurity. For simplicity, the delocalized nature of end-chain spins is not shown. By virtue of inter-chain coupling $`J_1`$ the end-chain spin on one of the fragments of the severed chain induces a “pocket” of staggered magnetization on a neighboring chain. The size of the pocket is determined by the correlation length $`\xi `$. The induced staggered moment interacts (again via $`J_1`$) with the end-chain spin on the other fragment. This adds up to an effective 2nd-order coupling between end-chain spins of the order of $`|J_{\mathrm{eff}}^{}|z\chi _\pi J_1^2`$, where $`z`$ is the inter-chain coordination number. In Ref.  it is argued that this coupling is antiferromagnetic and will bind the two end-chain spins into a non-magnetic singlet. As a result, the system will remain in a singlet ground state despite the introduction of a vacancy. Another way of looking at this phenomenon is that the effective nnn interaction $`J_{\mathrm{eff}}^{}`$ “patches” the broken Haldane chain at the impurity site with an antiferromagnetic bond, albeit a weak one.

We believe that this argument is, in general, incorrect. Indeed, from Fig. 4 it is clear that, no matter what the sign of $`J_1`$ is, $`J_{\mathrm{eff}}^{}`$ is actually ferromagnetic ($`J_{\mathrm{eff}}^{}<0`$) and, if anything, binds the two end-chain spins into a magnetic triplet, delocalized over the length scale $`\xi `$. Thus, in the absence of nnn interactions, long-range ordering should occur for arbitrary small impurity concentration $`x`$. The ordering temperature will, of course, be exponentially small with $`x`$.

This conclusion has to be modified in the presence of actual in-chain nnn AF interactions. For $`J^{}|J_{eff}^{}|`$ the arguments of Ref.  can be repeated verbatim to show that the system will remain disordered at $`T=0`$ for sufficiently small $`x`$. On the other hand, the doped system will order at a non-zero temperature for $`J^{}|J_{eff}^{}|`$. At any given small concentration there will be a quantum phase transition from a disordered to ordered state at a critical strength of $`J^{}`$ given by

$$J_c^{}z\chi _\pi J_1^2.$$
(12)
These estimates are semi-quantitative at best, yet they provide a simple physical picture and emphasize the relevant energy scales.

### B Application to PbNi2V2O8

In PbNi2V2O8, that remains in a gapped state at $`T0`$, long-range ordering can be induced by substituting some of the $`S=1`$ Ni2+ ions with non-magnetic Mg2+ ions. The transition was observed for all Mg-concentrations $`x`$ ranging from 0.1 down to 0.01. Though the Neel temperature $`T_\mathrm{N}`$ was found to decrease rapidly with decreasing $`x`$, there is, to date, no clear evidence of a critical concentration $`x_c`$ below which the system remains disordered. The fact that in PbNi2V2O8 doping-induced ordering is related to end-chain spins has been elegantly demonstrated by comparing the effects of Mg2+ ($`S=0`$) and Cu2+ ($`S=1/2`$) substitution on the Ni-sites. The latter is exactly half as effective, in terms of the effect on magnetic susceptibility, since the Cu-spin and the two end-chain spins become bound in a doublet state. In contrast with $`S=0`$ impurity substitution, the introduction of a $`S=1/2`$ impurity produces only one, not two, free spins in the system.

To estimate the critical strength of nnn AF interactions that would interfere with long-range ordering in weakly-doped PbNi2V2O8, we have to estimate $`J_{\mathrm{eff}}^{}`$. Using the experimental values of $`J`$, and assuming $`z=2`$, from Eqs. 12 and 57 we obtain $`J_c=|J_{\mathrm{eff}}^{}|=0.18(6)`$ meV. Given the experimental error bars and the semi-quantitative nature of Eq. 12, we can only conclude that in PbNi2V2O8 $`J^{}`$ is of the same order of magnitude as the critical value $`J_c^{}`$. The question of whether there is a critical doping level below which PbNi2V2O8 remains disordered, stays open. Hopefully, extending the $`T_\mathrm{N}(x)`$ measurements to lower concentrations, more accurate single-crystal neutron scattering measurements of the gap energies, and the development of a quantitative theoretical model will resolve this interesting issue in the future.

## V conclusions

i) This work provides a direct and accurate measurement of the zone-boundary energy in PbNi2V2O8 and SrNi2V2O8. The previously proposed model, based on weakly-interacting $`S=1`$ quantum spin chains with nearest-neighbor interactions, is shown to be an adequate description of these compounds in the entire dynamic range of magnetic excitations. ii) Unlike thought previously, weakly-coupled Haldane spin chains should order in three dimensions at arbitrary small concentrations of non-magnetic impurities. iii) Next-nearest neighbor AF interactions can inhibit long-range ordering at small impurity concentrations, but only if they exceed a certain critical value $`J_c^{}`$. iv) In the particular case of PbNi2V2O8 nnn interactions, if at all present, are no larger that several tenths of a meV, and may be of the same order of magnitude as $`J_c^{}`$.

###### Acknowledgements.

We would like to thank Y. Uchiyama and I. Tsukada for the collaboration in the early stage of this study, A. Tsvelik for fruitful discussions, and B. Taylor, R. Roth and K. Mohanty for technical support. Work at Brookhaven National Laboratory was carried out under Contract No. DE-AC02-98CH10886, Division of Material Science, U.S. Department of Energy. The work at the University of Tokyo was partially supported by Grant-in-Aid for COE Research “SCP Project” from the Ministry of Education, Science, Sports and Culture. Oak Ridge National Laboratory is managed for the U.S. D.O.E. by Lockheed Martin Energy Research Corporation under contract DE-AC05-96OR22464.


text:
# HIPASS Detection of an Intergalactic Gas Cloud in the NGC 2442 Group

## 1 Introduction

The nearby southern galaxy NGC 2442 ($`V_{}=1449`$ $`\mathrm{km}\mathrm{s}^1`$) provides a striking example of well-developed, but asymmetric spiral arms (Figure HIPASS Detection of an Intergalactic Gas Cloud in the NGC 2442 Group; see also Panel 207 of Sandage & Bedke (1994)). The northern arm bends back on itself quite tightly, and is sharply bisected by one continuous dust lane, while the southern arm is shorter, more open, and criss-crossed by numerous dust features. There have been several attempts to ascribe this asymmetry to a past interaction with a neighboring galaxy, but no clear culprit has yet emerged. Sandage & Bedke (1994) suggested the nearby E0 galaxy NGC 2434 ($`V_{}=1390`$ $`\mathrm{km}\mathrm{s}^1`$, $`\mathrm{separation}=16\stackrel{}{\mathrm{.}}8`$) or perhaps the more distant spiral NGC 2397 ($`V_{}=1363`$ $`\mathrm{km}\mathrm{s}^1`$, $`\mathrm{separation}=85\mathrm{}`$); Elmegreen et al. (1991) also blamed NGC 2434, but Mihos & Bothun (1997) argued on the basis of its proximity, and slightly disturbed morphology, that the small SB0/a galaxy AM 0738–692 is a more likely candidate. In an H I aperture synthesis study of NGC 2442 with the Australia Telescope Compact Array (ATCA), Houghton (1998) failed to detect any H I associated with AM 0738-692, but did find evidence of weakly-disturbed H I in the irregular galaxy ESO 059-G006 ($`V_{}=1346`$ $`\mathrm{km}\mathrm{s}^1`$) located some $`17\stackrel{}{\mathrm{.}}1`$ south of NGC 2442.

A deep optical image (Figure HIPASS Detection of an Intergalactic Gas Cloud in the NGC 2442 Group) showing NGC 2442 and most of the galaxies mentioned above has been produced by adding three contrast-enhanced IIIaJ plates from the UK Schmidt Telescope. This image also highlights a diffuse extension of the northern arm, which extends almost halfway towards ESO 059-G006 before becoming indistinguishable from the diffuse Galactic reflection nebulosity which pervades the field.

Distance estimates in the literature for NGC 2442 range from 14 Mpc (Sérsic & Donzelli, 1993), to 17.1 Mpc (Tully, 1988). We adopt $`D=15.5`$ Mpc ($`1\mathrm{}=4.5`$ kpc) in what follows. Garcia (1993) assigns only three other galaxies to the group (LGG 147) including NGC 2442: NGC 2397; NGC 2434; and PGC 20690 ($`=\mathrm{AM0720}723`$), a Magellanic irregular some $`3\stackrel{}{\mathrm{.}}5`$ south of NGC 2442. The same four galaxies were similarly grouped together by Tully (1988) and designated 53–17. However, in light of the above discussion, and several recent redshift determinations, it is clear that the NGC 2442 group is somewhat more populous than once thought. Table 1 lists all the galaxies within $`90\mathrm{}`$ (400 kpc projected distance) of NGC 2442 whose redshifts are known. The galaxies ESO 059-G010 and ESO 059-G007 are probably also close neighbors, but no redshifts are yet available.

In this paper, we report the discovery of significant amounts of diffuse H I to the northwest of NGC 2442, and discuss how this may help shed new light on the nature of its unusual appearance.

## 2 Observations and Results

The H I Parkes All-Sky Survey (HIPASS) is a blind H I survey of the entire local southern sky ($`\delta <2\mathrm{°}`$, $`1200<V_{}<12700`$ $`\mathrm{km}\mathrm{s}^1`$), carried out between February 1997 and March 2000. The survey uses a 13-beam receiver mounted on the Parkes 64 m telescope to scan the sky in $`8\mathrm{°}`$ strips of declination, revisiting each point on five separate occasions, for an effective integration time of 460 s beam-1. Full details of the multi-beam receiver system, survey strategy, and on-line data analysis can be found in Staveley-Smith et al. (1996) and Barnes et al. (2001). The HIPASS datacubes-11-11-11Spectra from the HIPASS data, gridded at $`8\mathrm{}`$ intervals, are now publicly accessible from http://www.atnf.csiro.au/research/multibeam/release/. used in this analysis have a final beamwidth of $`15\stackrel{}{\mathrm{.}}5`$, velocity resolution 18 $`\mathrm{km}\mathrm{s}^1`$, r.m.s. noise level of $`13`$ mJy beam-1, and a spatial pixel size of $`4\mathrm{}`$.

One of the first major data products to emerge from HIPASS is the HIPASS Bright Galaxy Catalog (Koribalski et al., 2001). In the course of measuring the centroid, total flux, and extent of H I emission from the 1000 brightest galaxies in the HIPASS data, we noticed an unusual extension to the integrated intensity (zeroth moment of the H I datacube) image of NGC 2442, with no obvious optical counterpart (Figure HIPASS Detection of an Intergalactic Gas Cloud in the NGC 2442 Group). The peak H I column density in this region is just 0.1 M pc-2, or $`1.3\times 10^{19}`$ cm-2. In accordance with the IAU standard, we assign this new object the provisional designation HIPASS J0731–69.

From an analysis of the zero moment map, as well as the total H I profile (Figure HIPASS Detection of an Intergalactic Gas Cloud in the NGC 2442 Group), we find this new source to have a flux integral of $`18\pm 3`$ Jy $`\mathrm{km}\mathrm{s}^1`$ (see Barnes et al. (2001) for a discussion of the uncertainties in total fluxes from the HIPASS data). Although the H I profile is somewhat asymmetric (and most unlike the “double-horned” or Gaussian profile expected of a spiral or a dwarf galaxy) we derive a heliocentric velocity for the peak of the line emission of $`1481\pm 18`$ $`\mathrm{km}\mathrm{s}^1`$, and profile width at 20 % of the peak intensity $`W_{20}=270\pm 60`$ $`\mathrm{km}\mathrm{s}^1`$. For comparison, the equivalent HIPASS parameters for NGC 2442 itself, from a robust moment analysis and from the H I profile (also slightly asymmetric – Fig. HIPASS Detection of an Intergalactic Gas Cloud in the NGC 2442 Group), are a flux integral of $`64\pm 4`$ Jy $`\mathrm{km}\mathrm{s}^1`$, heliocentric velocity of $`1456\pm 40`$ $`\mathrm{km}\mathrm{s}^1`$, and $`W_{20}=550\pm 12`$ $`\mathrm{km}\mathrm{s}^1`$.

The similarity in H I profile shapes raises the concern that HIPASS J0731–69 may just be a “ghost” of NGC 2442. We examined the five separate scans which pass through this region (each of which is displaced from the previous scan by $`7\mathrm{}`$; Barnes et al. 2001), and found a firm detection in each one of them, confirming that HIPASS J0731–69 is not an artifact of either the scanning, or the gridding process.

Of the 1000 or so brightest HIPASS sources examined by Koribalski et al. (2001), NGC 2442 was the only object found to have such a notable extension, with no optical counterpart. Most other extended sources turn out to be associated with known loose groups of galaxies. At least two sizable H I clouds lying outside the Zone of Avoidance, but apparently devoid of any stars, have also been found by visual inspection of the HIPASS data (Kilborn et al., 2000; Koribalski, 2001), but with velocities in the range $`400500`$ $`\mathrm{km}\mathrm{s}^1`$, they are probably part of, or just outside of the Local Group.

In an attempt to further characterize the nature of HIPASS J0731–69, we used the ATCA in its 750A configuration on 2000 August 25 UT to map the field of this emission. A mosaic of four pointings was observed over 10 hours, extending $`1\mathrm{°}`$ west and north of NGC 2442, and including NGC 2442 itself. These observations quite clearly show the H I in NGC 2442, but fail to detect any sign of HIPASS J0731–69, since the $`3\sigma `$ column density achieved with the limited number of baselines and amount of time available (i.e., $`10^{20}`$ cm-2) is still an order of magnitude greater than the peak column density found in the HIPASS data. We infer from this that the H I seen by HIPASS is indeed somewhat diffuse, and not clumpy (on the arcminute scales to which this ATCA configuration is most sensitive) as one might perhaps expect if it was associated with any particular optically-visible galaxy.

There are no cataloged galaxies within a $`15\mathrm{}`$ radius of the central position of HIPASS J0731–69 in the NASA Extragalactic Database (NED). We have examined the Digitised Sky Survey for potential optical counterparts to the newly-discovered gas. Although a total of seven independent objects were found within $`15\mathrm{}`$ of the central position, none are more than $`1\mathrm{}`$ in diameter, making it rather unlikely that any are new members of the NGC 2442 group (Garcia, 1993). The most likely candidate object, an isolated spiral galaxy at $`\alpha =07^\mathrm{h}31^\mathrm{m}10\stackrel{\mathrm{s}}{\mathrm{.}}4`$, $`\delta =68\mathrm{°}55\mathrm{}59\mathrm{}`$ (J2000), has been found to have a redshift inconsistent with that of the NGC 2442 group (I. Perez, private communication).

A channel map sequence (Figure HIPASS Detection of an Intergalactic Gas Cloud in the NGC 2442 Group) shows that even with the large HIPASS beamwidth and coarse velocity resolution, HIPASS J0731–69 is resolved into multiple components. Two components, centered on $`(\alpha ,\delta )=(07^\mathrm{h}34^\mathrm{m},69\mathrm{°}04\mathrm{})`$ and $`(07^\mathrm{h}31^\mathrm{m},69\mathrm{°}20\mathrm{})`$ (J2000), begin to appear at $`V_{}1325`$ $`\mathrm{km}\mathrm{s}^1`$ and then grow separately, before merging at a velocity near 1440 $`\mathrm{km}\mathrm{s}^1`$. A Gaussian fit to this merged component yields a centroid of $`\alpha =07^\mathrm{h}31^\mathrm{m}39\stackrel{\mathrm{s}}{\mathrm{.}}4`$, $`\delta =69\mathrm{°}01\mathrm{}36\mathrm{}`$ (J2000). At velocities near 1340 $`\mathrm{km}\mathrm{s}^1`$and 1380 $`\mathrm{km}\mathrm{s}^1`$, there are the first signs of a direct connection of the emission from the northeast and the southwest components respectively, to that from NGC 2442, and by 1418 $`\mathrm{km}\mathrm{s}^1`$, both components are apparently joined with NGC 2442 as well as with each other. The merged component of HIPASS J0731–69 again appears well connected with NGC 2442 for most velocities in the range 1500–1600 $`\mathrm{km}\mathrm{s}^1`$. However, it must be borne in mind that the separation of the two main components, as well as their distance from NGC 2442, are only $`12\times `$ the HIPASS beamwidth, so any apparent physical connections between each of them must be regarded with caution until higher resolution observations can be made.

## 3 Discussion

Since the bulk of this new gas lies at a projected separation of nearly $`40\mathrm{}`$ from NGC 2442 (and only slightly further from NGC 2397), it is probably best characterized as an intergalactic cloud associated with the NGC 2442 group, though both its velocity range, and Fig. HIPASS Detection of an Intergalactic Gas Cloud in the NGC 2442 Group, suggest the closest connection is to NGC 2442. One way to set limits on the actual separation of HIPASS J0731–69 from NGC 2442 would be to search for H$`\alpha `$ emission from the cloud’s surface, due to ultraviolet photons escaping from the galaxy. H$`\alpha `$ emission levels of 2–4 milliRayleighs are expected for a cloud within 200 kpc of the galaxy (Bland-Hawthorn et al., 1998), and emission measures as low as 1 milliRayleigh are now detectable with 4 m-class telescopes (J. Bland-Hawthorn, private communication).

Similar examples in nearby galaxy groups of intergalactic gas clouds without prominent optical counterparts include: the “Virgo Cloud” H I 1225+01 (Giovanelli & Haynes, 1989); the “Leo Ring” (Schneider et al., 1983); the H I cloud in the NGC 3256 Group (English, 1994); and the H I cloud superimposed on the blue compact dwarf FCC 35 in the Fornax Cluster (Putman et al., 1998). Table 2 compares the global properties of these four systems with HIPASS J0731–69. In each of the first two cases, faint dwarf galaxies apparently associated with these gas clouds have since been identified (Chengalur, Giovanelli, & Haynes, 1995; Schneider, 1989). In the absence of either a contrast-enhanced photo or a deep CCD image of the area covered by HIPASS J0731–69, we cannot rule out the possibility of a “Malin 1”-like object (Bothun et al., 1987), with an extended low surface brightness disk and very compact bulge component. However, the failure to detect gas clumping with the ATCA makes it rather less likely that any optical (galactic) counterpart exists.

The complex velocity field of HIPASS J0731–69 illustrated by Fig. HIPASS Detection of an Intergalactic Gas Cloud in the NGC 2442 Group, and the rather asymmetric H I profile of Fig. HIPASS Detection of an Intergalactic Gas Cloud in the NGC 2442 Group, are inconsistent with a uniformly rotating disk structure. Rather, it is more likely that HIPASS J0731–69 is an extended, stream-like feature of some kind. The H I velocity field for NGC 2442 of Houghton (1998) shows line-of-sight velocities ranging from just below 1100 $`\mathrm{km}\mathrm{s}^1`$ in the northeastern quadrant, to over 1700 $`\mathrm{km}\mathrm{s}^1`$ in the southwest. However, gas in the northern arm with velocities in the approximate range of $`14501600`$ $`\mathrm{km}\mathrm{s}^1`$ appears to be being drawn out and away from the rest of the disk in an easterly direction, toward the southwest component of HIPASS J0731–69. In addition, sections of the inner northern arm are displaced in velocity from the background disk by $`50`$ $`\mathrm{km}\mathrm{s}^1`$. Both these distortions are consistent with gas being pulled out of the plane, and their respective velocities and alignments lead us to hypothesize that they may, with better sensitivity and resolution, prove to be contiguous with the various components of HIPASS J0731–69. Similar anomalous velocities and double line profiles in the northern arm were seen optically by Bajaja, Agüero, & Paolantonio (1999) and by Mihos & Bothun (1997).

HIPASS J0731–69 spans $``$180 kpc, and contains $`10^9`$ M of atomic hydrogen, making it comparable in scale to both the Virgo Cloud and the Leo Ring. The Leo Ring contains as much gas as either of the nearby spiral galaxies M95 or M96, while the Magellanic system (LMC, SMC, Stream, Bridge, and Leading Arm) contains $`10^9`$ M of H I (Putman et al., 2001), or one-third as much gas as the Milky Way. The “forked” velocity structure of HIPASS J0731–69 is also not unlike that of the Magellanic Stream/Bridge/Leading Arm system (Putman et al., 1998), or the sections of the Leo Ring studied by Schneider (1985), which he was able to model quite well by a Keplerian orbit. The spatial resolution of the HIPASS data is insufficient to allow us to test any orbit models, but future mapping at higher resolution may enable the internal dynamics of this system to be understood.

Was the gas in HIPASS J0731–69 once part of NGC 2442? The H I mass-to-blue luminosity ratio for NGC 2442 alone ($`\mathrm{log}\mathrm{M}_{\mathrm{HI}}/\mathrm{L}_B=0.57`$, using our H I mass and $`B_\mathrm{T}=11.1`$ (Sérsic & Donzelli, 1993)), is identical to the median value for Sbc galaxies as a whole (Roberts & Haynes, 1994). However, compared with other Sbc galaxies, NGC 2442 is slightly sub-luminous optically, and its H I content is also just below the median. Including HIPASS J0731–69 would then make the total atomic gas content of NGC 2442 about average for an Sbc galaxy (this of course assumes that the disturbance to NGC 2442 has not radically changed its apparent Hubble type). The scatter in each parameter exhibited by Sbc galaxies ($`0.5`$ dex) means that global parameters are not particularly helpful in settling the question of the origin of HIPASS J0731–69.

We now briefly review how our discovery of HIPASS J0731–69 impacts upon the two main theories for the origin of the asymmetry in NGC 2442.

### 3.1 Tidal encounter?

As mentioned in the Introduction, there are several candidate galaxies for a recent past interaction with NGC 2442. What kind of close encounter might have caused NGC 2442 to lose up to 25% of its own H I content (or strip an equivalent amount of gas from some other galaxy), while only mildly affecting its global kinematics, infrared luminosity, and star formation rate (Mihos & Bothun, 1997)?

Large tidal tails and bridges on scales of 100 kpc or more are not uncommon in compact groups (e.g., the Leo triplet; Haynes, Giovanelli, & Roberts 1979), and in galactic mergers (see e.g., the merger “sequence” imaged in H I by Hibbard & van Gorkom (1996)). In such cases however, there are usually two or more galaxies of comparable mass involved, with separations of this same order. In addition, there is growing evidence for the formation of “tidal dwarf galaxies” from the gas in these tails (Duc & Mirabel, 1999). The resolution of the HIPASS data is not really adequate for us to determine whether HIPASS J0731–69 is an isolated cloud (or group of clouds), or if there is a bridge of material linking it directly to NGC 2442. We note that Houghton (1998), in her ATCA study of NGC 2442, observed a mild “stretching” of the H I disk towards the west, but did not find any signs of an H I tail within the $`34\mathrm{}`$ primary beam, making the former more likely.

Mihos & Bothun (1997) provide a model for a close passage of AM 0738–692 $``$200 Myr ago, which appears to be reasonably successful at explaining many of the morphological peculiarities of NGC 2442. There is no specific prediction of a large gas cloud like HIPASS J0731–69 in this scenario (but then neither was there any particular reason to expect such a feature). According to this model, the northern arm in particular is not driven by a spiral density wave, but is more like a tidal tail. Another good example of this phenomenon is the interacting pair NGC 6872/IC 4970 (Mihos, Bothun, & Richstone, 1993) which, just as in NGC 2442, exhibits large velocity gradients across the width of the tidal arm as material originally drawn away from the disk begins to fall back in, is compressed and shocked.

This mechanism accounts quite naturally for the prominent dust lane, enhanced star formation, and unusual line profiles observed in the northern arm of NGC 2442. However, in the Mihos & Bothun (1997) model, it is hard to believe that such an innocuous object as AM 0738–692 could once have contained as much gas as is in HIPASS J0731–69, or that it could have removed such a large quantity of gas to such a large distance. If NGC 2442 has been the victim of a recent tidal encounter, then the existence of HIPASS J0731–69 makes it far more likely that the culprit is a more massive galaxy to the northwest, such as the elliptical galaxy NGC 2434 or perhaps even the spiral + Magellanic Irregular pairing of NGC 2397 and NGC 2397A.

### 3.2 Ram pressure stripping?

The NGC 2442 system shares many of the characteristics of certain asymmetric galaxies, which are suspected of being ram-pressure stripped of some of their H I as they pass through the intracluster medium. Among the best examples of this phenomenon are NGC 2276 (Gruendl et al., 1993), NGC 4273 (Davis et al., 1995), NGC 7421 (Ryder et al., 1997), and the Virgo cluster galaxies NGC 4388 (Veilleux et al., 1999), NGC 4654 (Phookun & Mundy, 1995), and NGC 4522 (Vollmer et al., 2000). The most distinctive features of these objects are: (a) strongly asymmetric H I profiles; (b) a truncated optical and/or gas distribution on one side of the galaxy, often having a “bow shock”-like morphology; (c) a gaseous tail or “wake” in the opposite direction, possibly out of the disk plane; and (d) mild disturbance to an otherwise regular rotational velocity pattern. In some cases, the presence of a diffuse, hot gas component is confirmed by X-ray satellite observations, but the intracluster medium is not always so obvious.

The deep optical image of NGC 2442 (Fig. HIPASS Detection of an Intergalactic Gas Cloud in the NGC 2442 Group) shows quite a sharp edge to the northern edge of the stellar disk, when compared with the diffuse extensions to the southeast and southwest. The global asymmetry and sharp cutoff to the north side is even more pronounced in H$`\alpha `$ (Ryder & Dopita, 1993; Mihos & Bothun, 1997). The northern arm may represent the bow shock, where the disk is ploughing into the densest part of the inter-group medium, resulting in the observed velocity jumps. Disk rotation would stretch out the shock front, with the leading edge at the eastern end of the northern arm, and material flowing downstream along this arm to the west and south (this assumes the southern side of the disk of NGC 2442 is in the foreground). HIPASS J0731–69 would then represent a stream of gas trailing out behind (or completely detached from) NGC 2442, perhaps from a time when the galaxy experienced a particularly strong ram pressure. Both NGC 2442, and one of its companions to the south (ESO 059-G006), show distortions in their H I velocity fields (Houghton, 1998) which could just as easily be due to their passage through some inter-group medium, as to an interaction. A ROSAT HRI image barely detects NGC 2442 itself however, let alone any hot, diffuse gas.

It is interesting to contrast the global H I profile of NGC 2442 (Fig. HIPASS Detection of an Intergalactic Gas Cloud in the NGC 2442 Group) with the global 12CO ($`10`$) profile of Bajaja et al. (1995); the asymmetry in the atomic gas distribution is reversed in the molecular gas distribution, with the bulk of the CO emission apparently coming from the northeast sector of the disk, at lower velocities. Bajaja et al. (1995) did not survey the entire disk of NGC 2442, but were unable to detect 12CO ($`10`$) much beyond the point where the northern arm bends back on itself. Since the molecular gas is predominantly located in the inner disk, where the total mass surface density is highest (and ISM stripping least effective), any asymmetry in the CO distribution would tend to imply the action of tidal forces, rather than ran-pressure (Combes et al., 1988; Boselli et al., 1994). Since the ratio of molecular-to-atomic gas in NGC 2442 ($`M_{\mathrm{H2}}/M_{\mathrm{HI}}0.7`$) is close to the average for Sbc galaxies (Young & Scoville, 1991), this too weakens slightly the case for ram pressure stripping.

Clearly, HIPASS J0731–69 is an important new clue in the puzzle of NGC 2442’s disturbed appearance, as well as being a significant object in its own right on account of its apparent rarity in the HIPASS survey. Followup observations with more compact configurations of the ATCA are planned, as well as higher spectral resolution observations with the Parkes multi-beam system, and will be crucial in telling us more about the origin of HIPASS J0731–69 and its relationship to NGC 2442.

The Parkes telescope and Australia Telescope Compact Array are both part of the Australia Telescope, which is funded by the Commonwealth of Australia for operation as a national facility managed by CSIRO. We are grateful to the staff at the ATNF Parkes and Narrabri observatories, and to members of the Zone of Avoidance survey team, for assistance with HIPASS and follow-up observations. We especially wish to thank I. Perez, N. Killeen, T. Getts, and J. Bland-Hawthorn. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The Digitized Sky Surveys were produced at the Space Telescope Science Institute under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. The UK Schmidt Telescope was operated by the Royal Observatory Edinburgh, with funding from the UK Science and Engineering Research Council (later the UK Particle Physics and Astronomy Research Council), until 1988 June, and thereafter by the Anglo-Australian Observatory.