text:
# Transversal inhomogeneities in dilute vibrofluidized granular fluids

## I Introduction

Granular materials are assemblies of macroscopic particles dissipating their energy through inelastic collisions . They exhibit a very rich phenomenology that is only partially understood. One of the most peculiar behaviors of granular systems, which has attracted a lot of attention in the last years, is their tendency to spontaneously develop strong spatial inhomogeneities. In many different situations, the density shows a sharp profile that is not induced by the boundary conditions. This phenomenon is often referred to as a clustering effect , since high density regions coexist in the system with regions where the density is very low.

In vibrated granular systems, clustering effects show up in many cases as a spontaneous symmetry breaking in the direction parallel to the vibrating wall. Consider a gas enclosed in a box that is being supplied energy through a vibrating wall located at $`x=0`$. There are no other external forces acting on the system. The box is divided into two equal compartments by a wall normal to the $`x`$ axis starting at a certain height. At sufficiently low average density, the hydrodynamic fields are symmetric to both sides of the partition, but above a critical average density, which depends on the value of the restitution coefficient, an asymmetry in the number of particles at each side of the container occurs . This asymmetry has been shown to be associated with a bifurcation of the solution of the hydrodynamic equations describing the state of the system.

A similar symmetry breaking has been observed in a system in presence of a gravitational force acting in the $`x`$ direction. In this case, the system is unbounded for $`x>0`$, and the partition has a hole at a certain height. Again, an asymmetry in the number of particles in the two compartments develops if a control parameter, dependent on the amplitude of the vibration and the degree of inelasticity, is larger than a critical value .

Symmetry breaking in the direction parallel to the vibrating wall has also been observed in systems without any partition of the container. Sunthar and Kumaran  have reported molecular dynamics simulation results showing the presence of convection rolls and phase separation into coexisting dense and dilute regions in a granular system in presence of gravity. The phase separation takes place on the surface of the vibrating wall and, as already indicated, in the direction perpendicular to the energy flux. No theoretical explanation for this phenomenon is provided in ref. , although the effect of the different parameters controlling the behavior of the system is discussed in detail, on the basis of the simulation results. For a two-dimensional closed system in absence of gravity, a transversal continuous spontaneous symmetry breaking has also been predicted . Not at all surprisingly, the gradients are now sharper next to the elastic wall, opposite to the energy source. This is consistent with the positions that must have the holes in systems with a separating wall in order to observe the symmetry breaking, with and without gravitational field acting on the system. The work by Livne et al.  is restricted to the nearly elastic limit and it is based on a numerical marginal stability analysis of the hydrodynamic equations. The predictions of this analysis are compared with numerical solutions of the own hydrodynamic equations with the appropriate boundary conditions.

In this paper, the bifurcation predicted in ref.  will be considered again. There are several reasons for that. Firstly, hydrodynamic equations derived from the Boltzmann equation for smooth inelastic hard disks and valid, in principle, for arbitrary inelasticity will be used, then somehow extending the previous results. Secondly, instead of a numerical analysis of the stability of the solutions of the hydrodynamic equations, an analytical study, based on the WKB approximation, will be presented here. One of the main advantages of this approach is that the dimensionless control parameter governing the bifurcation phenomenon is clearly identified. A third motivation for the present work is to report simulation results, both from molecular dynamics and also from Monte Carlo simulation of the Boltzmann equation, showing the existence of the predicted transition. Since these simulation techniques do not contain any externally introduced hydrodynamic concept, they provide a direct proof of the existence of the continuous symmetry breaking, and a test of the theoretical predictions. Attention will be also paid to the form of the hydrodynamic profiles beyond the bifurcation. This leads to a deeper understanding of the development of the instability. In any case, it is clear that the work in ref.  opened the way to more systematic investigations of this instability, like the one in this paper.

The plan of the paper is as follows. In Sec. II, the hydrodynamic description of the one-dimensional state of a low density vibrofluidized granular gas will be briefly summarized. The analytical expressions of the hydrodynamic profiles are given. This state is the starting point for the marginal stability analysis developed in Sec. III. Linearization of the hydrodynamic equations around the one-dimensional state leads to a second order linear differential equation. The WKB solution of the closed problem posed by this equation and the corresponding boundary conditions is built up. This requires to consider three different cases, depending on the value of the parameters characterizing the system. From the WKB solution, the marginal stability curve follows easily. Simulation results are presented and compared with the theoretical predictions in Sec. IV, where an order parameter characterizing the transition is defined. A good agreement is found. The last section contains some final remarks, as well as the comparison of the results derived here with those in ref.  in the common range of applicability, namely nearly elastic collisions and very low density.

## II Basic equations and the reference state

For a steady state without macroscopic flows, the balance equations for a two-dimensional gas of smooth inelastic hard disks of mass $`m`$ and diameter $`\sigma `$ have the form

$$𝖯=0,$$
(1)
$$𝐪+nT\zeta =0,$$
(2)
where $`n`$ and $`T`$ are the number density and granular temperature (with Boltzmann’s constant set equal to unity), respectively. In the low density limit, for a gas described by the inelastic Boltzmann equation, and to lowest order in the gradients (Navier-Stokes order), the pressure tensor $`𝖯`$ and heat flux $`𝐪`$, for the steady state under consideration are given by

$$𝖯=p𝖨,$$
(3)
$$𝐪=\kappa T\mu n,$$
(4)
$`𝖨`$ being the unit tensor, $`p=nT`$ the hydrodynamic pressure, $`\kappa `$ the heat conductivity, and $`\mu `$ a transport coefficient that has no analogue in the elastic case. These transport coefficients are proportional to the elastic heat conductivity $`\kappa _0(T)`$,

$$\kappa =\kappa ^{}(\alpha )\kappa _0(T),\mu =\mu ^{}(\alpha )\frac{T\kappa _0(T)}{n},$$
(5)
$$\kappa _0(T)=\frac{2}{\sigma }\left(\frac{T}{\pi m}\right)^{1/2}.$$
(6)
Finally, the cooling rate in the same approximation is related with the elastic shear viscosity $`\eta _0`$ by

$$\zeta =\zeta ^{}(\alpha )\frac{p}{\eta _0},$$
(7)
$$\eta _0(T)=\frac{1}{2\sigma }\left(\frac{mT}{\pi }\right)^{1/2}.$$
(8)
The functions $`\kappa ^{}`$, $`\mu ^{}`$, and $`\zeta ^{}`$ only depend on the constant coefficient of normal restitution $`\alpha `$ characterizing the inelasticity of collisions. Their explicit expressions are given in refs.  and , and will be not reproduced here. By using the above expressions we get from Eqs. (1) and (2),

$$\frac{p}{x}=\frac{p}{y}=0,$$
(9)
$$\left(\kappa ^{}\mu ^{}\right)\left[\frac{}{x}\left(\kappa _0\frac{T}{x}\right)+\frac{}{y}\left(\kappa _0\frac{T}{y}\right)\right]\zeta ^{}\frac{p^2}{\eta _0}=0.$$
(10)
Next, we specify the boundary conditions. We consider that the system has $`N`$ particles enclosed in a rectangular box with dimensions $`L_x`$ and $`L_y`$. The wall located at $`x=0`$ is vibrating and, therefore, supplies energy to the system. This energy is needed in order to keep and sustain a fluidized steady state. The other three walls, located at $`x=L_x`$, $`y=0`$, and $`y=L_y`$, respectively, are at rest. For the sake of simplicity, collisions of particles with all four walls are assumed to be elastic. Then, the mathematical boundary conditions to be imposed are:

$$\left(\frac{T}{x}\right)_{x=L_x}=\left(\frac{T}{y}\right)_{y=0}=\left(\frac{T}{y}\right)_{y=L_y}=0,$$
(11)
$$\left[\kappa ^{}(\alpha )\mu ^{}(\alpha )\right]\left[\kappa _0(T)\frac{T}{x}\right]_{x=0}=Q.$$
(12)
Equations (11) express that the heat flux must vanish at the inmobile walls, while Eq. (12) is the energy balance at the vibrating wall. The quantity $`Q`$ is the rate of energy input through this wall per unit of length. Its calculation in terms of the parameters defining the motion of the wall has been addressed in several works. Here we will consider the simplest possibility, namely that the wall moves in a sawtooth manner with velocity $`v_b`$ . This is a good approximation to more realistic motions, as long as the characteristic frequency of vibration of the wall be much larger than the collision rate of the gas molecules in its vicinity. In addition, it will be assumed that the amplitude of the vibration is much smaller than the mean free path of the particles of the gas next to it, so collective motions in the system are not being generated. Under the above conditions, it is

$$Q=pv_b.$$
(13)
The closed mathematical problem defined by Eqs. (9)–(13) admits a $`y`$-independent solution $`T=T_R(x)`$ that has been discussed in detail in ref. . The existence of this one-dimensional solution had been previously noticed by Grossman et al. . In the following, we will refer to this solution as the reference state and its properties will be characterized with a subindex $`R`$. The reference temperature profile is

$$T_R(\xi _x)=T_{0,R}\left[\frac{\mathrm{cosh}(\xi _x^{}\xi _x)}{\mathrm{cosh}\xi _x^{}}\right]^2,$$
(14)
where $`\xi _x`$ is a dimensionless scaled length defined by

$$\xi _x=\sqrt{a(\alpha )}_0^x\frac{dx}{\lambda _R(x)},$$
(15)
with $`\lambda _R(x)`$ being the local mean free path,

$$\lambda _R(x)=\frac{1}{2\sqrt{2}\sigma n_R(x)},$$
(16)
and

$$a(\alpha )=\frac{\pi \zeta ^{}}{16(\kappa ^{}\mu ^{})}.$$
(17)
Moreover, $`\xi _x^{}`$ is the value of $`\xi _x`$ for $`x=L_x`$, i.e.,

$$\xi _x^{}=2\sqrt{2a(\alpha )}\frac{\sigma N}{L_y}.$$
(18)
Therefore, $`\xi _x^{}`$ is proportional to the number of monolayers of particles perpendicular to the $`x`$-axis at rest, $`\sigma N/L_y`$. Finally, the uniform pressure of the reference state is

$$p_R=\frac{T_{0,R}(2\xi _x^{}+\mathrm{sinh}2\xi _x^{})}{8\sqrt{2a(\alpha )}\sigma L_x\mathrm{cosh}^2\xi _x^{}},$$
(19)
and the temperature $`T_{0,R}`$ of the gas next to the vibrating wall is given by

$$T_{0,R}=\left(\frac{e}{\mathrm{tanh}\xi _x^{}}\right)^2,e(\alpha )=\left(\frac{2am}{\pi }\right)^{1/2}\frac{v_b}{\zeta ^{}}.$$
(20)
The simplicity of the above results is a consequence of the limiting kind of motion of the vibrating wall considered, whose only only effect is to transfer energy to the grains, without inducing any periodic motion in the system. In some previous studies , a “thermal” wall, instead of a vibrating one, has been considered at $`x=0`$. By definition, particles which collide with a thermal wall leave it with the velocity distribution corresponding to the temperature of the wall. Although this kind of walls are far from reality for granular systems, their consideration might be useful for comparison purposes. The only change to be made in the above discussion in order to apply it to a system with a thermal wall, is to replace the boundary condition given in Eq. (12) by the requirement that the temperature $`T_{0,R}`$ has the value determined by the wall. As a consequence, the expressions of the hydrodynamic profiles remain the same, while Eq. (20) does not apply in this case.

## III Marginal stability analysis

Our aim now is to investigate whether the system described in the previous Section exhibits another steady state, in addition to the reference one. Then, we introduce a perturbation $`\delta T(x,y)`$ by

$$T(x,y)=T_R(x)+\delta T(x,y),$$
(21)
with $`\delta T(x,y)T_R(x)`$, and search for a solution of this form to Eqs. (9)-(10) with the boundary conditions (11) and (12). Substitution of Eq. (21) into Eq. (10) and linearization in $`\delta T`$ yields

$$\left[\frac{^2}{\delta \xi _x^2}+2\frac{1}{2}\left(\frac{\mathrm{ln}T_R}{\xi _x}\right)^2+\left(\frac{\lambda _R}{\overline{\lambda }}\right)^2\frac{^2}{\xi _y^2}\right]\delta T=\frac{4T_R}{p_R}\delta p,$$
(22)
where $`\delta p`$ has been defined by $`p=p_R+\delta p`$, $`\xi _x`$ is the dimensionless scale introduced in Eq. (15), and

$$\xi _y=\sqrt{a(\alpha )}\frac{y}{\overline{\lambda }}.$$
(23)
Here $`\overline{\lambda }`$ is the average mean free path,

$$\overline{\lambda }=\frac{1}{2\sqrt{2}\sigma \overline{n}},$$
(24)
$`\overline{n}=N/L_xL_y`$. Next, we consider factorized solutions of Eq. (22),

$$\delta T(\xi _x,\xi _y)=\varphi (\xi _x)\phi (\xi _y).$$
(25)
The conditions at the boundaries to be verified by the functions $`\varphi `$ and $`\phi `$ are

$$\left(\frac{\varphi }{\xi _x}\right)_{\xi _x=\xi _x^{}}=\left(\frac{\phi }{\xi _y}\right)_{\xi _y=0}=\left(\frac{\phi }{\xi _y}\right)_{\xi _y=\xi _y^{}}=0,$$
(26)
$$\frac{1}{2}\left(\frac{\mathrm{ln}T_R}{\xi _x}\right)_{\xi _x=0}\varphi (0)+\left(\frac{\varphi }{\xi _x}\right)_{\xi _x=0}=0,$$
(27)
where

$$\xi _y^{}=2\sqrt{2a(\alpha )}\frac{\sigma N}{L_x}.$$
(28)
Moreover, Eqs. (11) and (12) also imply that it must be $`\delta p=0`$, i.e. the pressure is not changed by the small perturbation. When Eq. (25) is substituted into Eq. (22), use of separation of variables leads to the equations

$$\frac{1}{\phi (\xi _y)}\frac{d^2\phi (\xi _y)}{d\xi _y^2}=k^2,$$
(29)
$$\left[\frac{d^2}{d\xi _x^2}+2\frac{1}{2}\left(\frac{d\mathrm{ln}T_R}{d\xi _x}\right)^2\left(\frac{k\overline{n}T_R}{p_R}\right)^2\right]\varphi (\xi _x)=0,$$
(30)
where $`k`$ is the constant of separation. The solution of Eq. (29) verifying the corresponding boundary conditions in Eq. (26) is

$$\phi =A\mathrm{cos}k\xi _y,$$
(31)
with $`A`$ an arbitrary constant. Moreover, the values of $`k`$ are restricted to $`k=\pi q/\xi _y^{}`$, $`q`$ being an integer. It is important to realize that the low density limit does not imply either $`\xi _x^{}1`$ or $`\xi _y^{}1`$. In fact, both quantities can be large in a very dilute system. An additional restriction to be required to $`\delta T(x,y)`$ is that the total number of particles in the system, $`N`$, is conserved. It is easily seen that this condition is equivalent to

$$_0^{\xi _x^{}}𝑑\xi _x_0^{\xi _y^{}}𝑑\xi _y\frac{\delta T(\xi _x,\xi _y)}{T_R(\xi _x)}=0.$$
(32)
Equation (25) with $`\phi (\xi _y)`$ given by Eq. (31) guarantees that this equality is verified.

An equation having a structure similar to Eq. (30) was obtained in ref.  for a dense system in the limit of nearly elastic collisions, by employing approximate constitutive relations introduced by Grossman et al. . While the equation in  was solved using numerical techniques, here we will use a Wentzel, Kramers, Brillouin (WKB) approximation  to investigate the possible solutions . First we use Eqs. (14) and (19) to rewrite Eq. (30) as

$$\frac{d^2\varphi (\xi _x)}{d\xi _x^2}+f(\xi _x)\varphi (\xi _x)=0,$$
(33)
with

$$f(\xi _x)=\frac{2}{\mathrm{cosh}^2(\xi _x^{}\xi _x)}\frac{16k^2\xi _x^2\mathrm{cosh}^4(\xi _x^{}\xi _x)}{(2\xi _x^{}+\mathrm{sinh}2\xi _x^{})^2},$$
(34)
which is a monotonic increasing function of $`\xi _x`$ in all the interval $`0\xi _x\xi _x^{}`$.

To construct the WKB exponential approximation of Eq. (33) it is necessary to consider three different ranges of parameters, that will be discussed separately in the following.
a) The function $`f(\xi _x)`$ is positive everywhere in the system. This is equivalent to $`f(0)>0`$ or

$$k\frac{2\xi _x^{}+\mathrm{sinh}2\xi _x^{}}{2\sqrt{2}\xi _x^{}\mathrm{cosh}^3\xi _x^{}}.$$
(35)
Then, the WKB solution is oscillatory,

$$\varphi (\xi _x)=\frac{c_1}{\sqrt{g(\xi _x)}}\mathrm{exp}\left[i_0^{\xi _x}𝑑\xi _x^{}g(\xi _x^{})\right]+\frac{c_2}{\sqrt{g(\xi _x)}}\mathrm{exp}\left[i_0^{\xi _x}𝑑\xi _x^{}g(\xi _x^{})\right],$$
(36)
where $`c_1`$ and $`c_2`$ are constants, and

$$g(\xi _x)=\sqrt{f(\xi _x)}.$$
(37)
Imposing the boundary condition at $`\xi _x=\xi _x^{}`$ leads to

$$\varphi (\xi _x)=\frac{c}{\sqrt{g(\xi _x)}}\mathrm{cos}_{\xi _x}^{\xi _x^{}}𝑑\xi _x^{}g(\xi _x^{}),$$
(38)
with $`c`$ another constant. When the boundary condition at $`\xi _x=0`$, Eq. (27), is also required, the consistency condition

$$\mathrm{tan}_0^{\xi _x^{}}𝑑\xi _xg(\xi _x)=\frac{3\mathrm{tanh}\xi _x^{}}{g^3(0)\mathrm{cosh}^2\xi _x^{}}$$
(39)
follows. The above equation determines the possible values of the parameters for which a WKB solution of the differential equation (33) exists in the region under study.

If a thermal wall is considered at $`x=0`$, the boundary condition (27) must be replaced by $`\varphi (0)=0`$, and instead of Eq. (39) it is found:

$$\mathrm{cos}_0^{\xi _x^{}}𝑑\xi _xg(\xi _x)=0,$$
(40)
i.e.,

$$_0^{\xi _x^{}}𝑑\xi _xg(\xi _x)=\frac{(2q+1)\pi }{2},$$
(41)
where $`q`$ is an arbitrary integer.
b) The function $`f(\xi _x)`$ is negative everywhere in the system. This is the case if $`f(\xi _x^{})0`$. Therefore, the region of parameters being considered is defined by

$$k\frac{2\xi _x^{}+\mathrm{sinh}2\xi _x^{}}{2\sqrt{2}\xi _x^{}}.$$
(42)
The exponential WKB approximation in this case reads

$$\varphi (\xi _x)=\frac{b_1}{\sqrt{h(\xi _x)}}\mathrm{exp}\left[_0^{\xi _x}𝑑\xi _x^{}h(\xi _x^{})\right]+\frac{b_2}{\sqrt{h(\xi _x)}}\mathrm{exp}\left[_0^{\xi _x}𝑑\xi _x^{}h(\xi _x^{})\right],$$
(43)
with

$$h(\xi _x)=\sqrt{f(\xi _x)},$$
(44)
and $`b_1`$ and $`b_2`$ arbitrary constants. Imposing the boundary conditions (26) and (27) leads to the relationship

$$h(0)^3\mathrm{tanh}_0^{\xi _x^{}}𝑑\xi _xh(\xi _x)=\frac{3\mathrm{tanh}\xi _x^{}}{\mathrm{cosh}^2\xi _x^{}}.$$
(45)
For a thermal wall at $`x=0`$, the above equation is substituted by

$$\mathrm{cosh}_0^{\xi _x^{}}𝑑\xi _xh(\xi _x)=0,$$
(46)
and, therefore, there is no WKB solution in this range of parameters.
c) The function $`f(\xi _x)`$ changes sign in the interval $`0<\xi _x<\xi _x^{}`$. It must be $`f(0)<0`$ and $`f(\xi _x^{})>0`$. This corresponds to the $`k`$ interval

$$\frac{2\xi _x^{}+\mathrm{sinh}2\xi _x^{}}{2\sqrt{2}\xi _x^{}\mathrm{cosh}^32\xi _x^{}}k\frac{2\xi _x^{}+\mathrm{sinh}2\xi _x^{}}{2\sqrt{2}\xi _x^{}}.$$
(47)
Since $`f(\xi _x)`$ exhibits in this case a zero in the integration range, we have to consider separately the regions $`\xi _x<a`$ and $`\xi _x>a`$, where $`a`$ is the turning point, i.e. $`f(a)=0`$. In the former region, the WKB solution is given by an expression of the form (43), while in the latter it has the form (36). A global solution is constructed by matching both WKB approximations through the connection formulas, expressing the connection between the oscillatory and exponential behaviors. Since it is easily seen that the turning point is a simple (first-order) zero, a standard application of the theory suffices . Then, imposing the boundary conditions gives the following equation to be verified by the solutions in this range:

$$\mathrm{tan}\left(\theta \frac{\pi }{4}\right)\frac{3\mathrm{tanh}\xi _x^{}h^3(0)\mathrm{cosh}^2\xi _x^{}}{3\mathrm{tanh}\xi _x^{}+h^3(0)\mathrm{cosh}^2\xi _x^{}}=\frac{1}{2}\mathrm{exp}\left[2_0^a𝑑\xi _xh(\xi _x)\right],$$
(48)
where

$$\theta =_a^{\xi _x^{}}𝑑\xi _xg(\xi _x).$$
(49)
The result for a thermal wall at $`x=0`$ is given in this case by

$$\mathrm{tan}\left(\theta \frac{\pi }{4}\right)=\frac{1}{2}\mathrm{exp}\left[2_0^a𝑑\xi _xh(\xi _x)\right].$$
(50)
Once we have derived the equations determining the solutions to the boundary value problem, the strategy to build up the marginal stability curve is as follows. Given a value of $`\xi _x^{}`$, the first question is to see whether there is a bifurcation from the one-dimensional solution, i.e. whether Eq. (30) has a solution. In case the answer to this question is affirmative, the smallest value of $`\xi _y^{}`$ for which there is solution gives the point of the marginal stability curve corresponding to that value of $`\xi _x^{}`$. For larger values of $`\xi _y^{}`$, the one-dimensional solution is not stable. Suppose a solution of Eq. (33) exists for a given wavenumber $`k`$. This value is compatible with many (infinite) values of $`\xi _y^{}`$, the smallest one corresponding to the choice $`q=1`$ in the relationship between the possible values of $`k`$ and $`\xi _y^{}`$ (see below Eq. (31)). Moreover, the larger the value of $`k`$, the smaller the value of $`\xi _y^{}`$ associated. The conclusion is that the stability curve is determined by the largest value of $`k`$ for which the problem has solution. From a practical point of view, we have to start by looking for solutions belonging to case b discussed above. The values $`\xi _y^{}=\pi /k`$ obtained from the solutions of Eq. (45) are the critical values, i.e., they define points on the marginal stability curve. Nevertheless, there is no solution for every value of $`\xi _x^{}`$ belonging to the $`k`$-interval defining case b. More precisely, there is no solution at all for a thermal wall, while it is a simple matter to show that the existence of solution in this region for a vibrating wall requires that

$$\left[h^{(1)}(0)\right]^3\mathrm{tanh}_0^{\xi _x^{}}𝑑\xi _xh^{(1)}(\xi _x)\frac{3\mathrm{tanh}\xi _x^{}}{\mathrm{cosh}^2\xi _x^{}},$$
(51)
where

$$h^{(1)}(\xi _x)=\underset{kk_1}{lim}h(\xi _x).$$
(52)
Equation $`(\text{51})`$ gives $`\xi _x^{}0.555`$. Therefore, solutions for larger values of $`\xi _x^{}`$ must be found, if they exist, in a different region of values of $`k`$, namely in that considered in case c. With this procedure, it is an easy task to find the largest value of $`k`$ (smallest value of $`\xi _y^{}`$) for which the eigenvalue problem defined by Eq. (30) has solution for each value of $`\xi _x^{}`$, in the WKB approximation. The marginal stability curve for a vibrating wall obtained in this way is given in Fig. 1, where we have plotted the critical values $`\mathrm{\Delta }_c`$ of the asymmetry parameter $`\mathrm{\Delta }L_y/L_x=\xi _y^{}/\xi _x^{}`$, as a function of $`\xi _x^{}`$. The solid line is the WKB solution corresponding to the so-called case b, while the dashed line is from case c. Note that the solutions from both regions of parameters match rather smoothly at $`\xi _x^{}0.555`$. Moreover, the critical asymmetry $`\mathrm{\Delta }_c`$ is a monotonic decreasing function of $`\xi _x^{}`$, growing very fast when $`\xi _x^{}`$ tends to zero, as it was expected since the transversal inhomogeneities must disappear in the elastic limit.

For a thermal wall, we have obtained that there is no solution in region b. Therefore, we must search for solutions belonging to the range c of parameters. Analysis of Eq. (50) leads to the result that the equation has no solution for $`\xi _x^{}1.181`$. Moreover, Eq. (41) has no solution in that region either, implying that, in the WKB approximation, the transition to the transversally inhomogeneous state requires, in the case of system driven by a thermal wall, that the inelasticity and the number of monolayers at rest be not small. The dotted line in Fig. 1 shows the marginal stability curve for a thermal wall. In the limit of large $`\xi _x^{}`$, the curve overlaps with the one for a vibrated system.

## IV Simulation results

To test the above theoretical results and also to investigate the nature of the predicted symmetry breaking, two more fundamental descriptions of the system, via the direct simulation Monte Carlo (DSMC) method and molecular dynamics (MD) simulation, respectively, have been considered. Although both are based on a dynamical simulation of the system, their nature is rather different. The DSMC method provides an algorithm to obtain numerical solutions of the Boltzmann equation for given initial and boundary conditions . Therefore, it relies on the validity of a kinetic theory description, namely the one given by the Boltzmann equation. On the other hand, no hydrodynamic approximations are introduced into the description, so that the validity of a hydrodynamic level of description, as provided for instance by the Navier-Stokes equations, is not taken for granted in this approach.

The MD simulations follow the motion of the particles of the system as a sequence of free motions and binary collisions , i.e. by direct application of Newton’s equations of motion. Therefore, MD provides the more basic description of the evolution of the system. In this context, it is important to stress that although the DSMC method also uses “particles” at a formal level, they are not real particles, but fictitious ones, which are employed to mimic the ideal dynamics described by the Boltzmann equation. It is precisely the ideal nature of the particles in the DSMC method what allows a very high numerical accuracy, since the number of particles used do not affect at all the physical state being simulated. In particular, the above number is not related with the actual density of the system.

In the following, we will report simulation results obtained for the system studied in the previous sections, restricting ourselves to the vibrating wall. We have used the two methods, DSMC and MD, because they complement one another. For instance, comparison of the results obtained by both methods provides a test of the validity of the inelastic Boltzmann equation to describe a system under the physical conditions used in the MD simulations. Since the one-dimensional state has been discussed in detail elsewhere , attention will be focused here on the behavior of the system in the vicinity of the bifurcation.

In the MD simulations, the density plays of course a relevant role. As the interest here is in the low density limit, small values of the surface fraction $`\nu =N\pi \sigma ^2/4L_xLy`$, typically of the order of $`10^2`$, have been used. For fixed given values of $`\nu `$, $`\alpha `$ and $`L_x`$, a set of simulations have been run corresponding to different widths $`L_y`$ of the system, starting from a small enough value. This means that each set of simulations corresponds to the same value of $`\xi _x^{}`$, as defined in Eq. (28), while they differ in the asymmetry parameter $`\mathrm{\Delta }`$. The simulations of the Boltzmann equation have been carried out with a similar systematic, the main difference being that the density plays no role in them. The other parameter needed to specify the simulations is the velocity of the vibrating wall $`v_b`$. We have verified that, in agreement with the theory developed in the previous sections, the formation of transversal inhomogeneities is not altered by modifying this velocity, as long as it is large enough to fluidize the complete system.

All the simulations started from a spatially homogeneous configuration, with a Gaussian velocity distribution, corresponding to an arbitrary temperature. The simulation is then followed until the system reaches a steady state, in which the several monitored properties of the system (mean kinetic energy, density fluctuations, and hydrodynamic profiles) become time independent. Once the system is in the steady state, all statistical averages of interest are accumulated. For the purpose here, the density and temperature profiles provide the relevant information. Let us describe what is observed in a set of simulations, namely we are going to present DSMC data from systems with $`\alpha =0.95`$ and $`L_x=10\overline{\lambda }`$. This corresponds to $`\xi _x^{}=1.015`$. In Fig. 2 the steady two-dimensional density profile for $`L_y=20\overline{\lambda }`$ ($`\mathrm{\Delta }=2`$) is shown in a three-dimensional plot. It is seen that no gradients in the $`y`$-direction are present, i.e. the system is in the one-dimensional reference state. The transversal homogeneity appears even clearer in Fig. 3, where the density is plotted as a function of $`y`$ for several fixed values of $`x`$.

When the width $`L_y`$ is increased keeping all the other parameters fixed, a critical value shows up such that gradients in the $`y`$-direction spontaneously develop in the system for larger widths. An example of this is given in Fig. 4, where the density surface for the same parameters as in Fig. 2, except that now $`L_y=26.5\overline{\lambda }`$, is plotted. A gradient in the direction parallel to the vibrating wall is clearly identified, becoming more pronounced next to the opposite wall, as illustrated in Fig. 5. The density gradients in this case are relatively small, being the maximum variation of the density in the $`y`$ direction of the order of $`5\%`$. In fact, in all the simulations, both by DSMC and MD, it has been found that when increasing $`L_y`$, a continuous transition from the one-dimensional state to a state with weak inhomogeneities in the transversal direction occurs. Moreover, the transversal density profile near the transition exhibits, as in the case of Fig. 5, a wavelength equal to twice the width of the system. In the language used in Sec. III, what is observed is the perturbation with $`q=1`$ ($`k=\pi /\xi _y^{}`$), consistently with our theoretical analysis.

If the width of the system is increased further on, the transversal inhomogeneities grow very fast and, of course, the results from the linear marginal stability analysis do not apply. Simulations show that the density becomes very large in one of the corners of the system, away from the vibrating wall. The rapid increase of the gradients and the large value of the density in a localized region of the system, lead to conclude that in this regime a cluster of particles is formed . Of course, for such a region of parameters the hydrodynamic profiles obtained from DSMC and MD are quite different, as the former considers the particles as points while the latter assigns them a finite diameter, implying that the density is bounded by the close-packing value.

Once the appearance of transversal inhomogeneities has been observed by visual inspection, it is desirable to have a “quantitative” criterium to establish whether the system is or is not transversally homogeneous. This is equivalent to identifying an order parameter to characterize the transition. Since there may be gradients in both directions, the identification of such a parameter is not at all a trivial task. The one we have chosen is defined as follows. First, we introduce the dimensionless quantity $`\rho _x(y)`$ by

$$\rho _x(y)=\frac{n(x,y)}{L_y^1_0^{L_y}𝑑yn(x,y)}1.$$
(53)
If the system is transversally homogeneous, $`\rho _x(y)`$ is independent of both $`x`$ and $`y`$, and equal to zero. When transversal inhomogeneities are present, it depends of course on $`y`$ but, as it can be guessed from Figs. 4 and 5, it also depends on $`x`$. Nevertheless, the simulation data indicate that this dependence is rather weak, at least near the transition. As an example, in Fig. 6 the function $`\rho _x(y)`$ has been plotted for the same system as in Figs. 4 and 5. The different lines correspond to four different values of $`x`$, equally separated, in the interval $`[3L_x/4,L_x]`$. This is the region where the transversal gradients are larger. From the figure it follows that the $`x`$ dependence is practically scaled out in the definition of $`\rho _x(y)`$. It must be mentioned, however, that if the whole range of variation of $`x`$ is considered, some dependence on $`x`$ shows up. In any case, the departure from zero of the average value of $`\rho _x(y)`$, $`\overline{\rho }(y)`$, over a certain $`x`$ interval, next to the vibrating wall and not too large to avoid $`x`$ dependence, provides a good criterium to distinguish transversally inhomogeneous systems from homogeneous ones. The results to be discussed in the following have been obtained using the interval $`[3L_x/4,L_x]`$.

The above discussion, the theoretical analysis, and the numerical results, as those illustrated in Fig. 6, suggest that a good order parameter may be given by the absolute value of the first Fourier component, $`|f_1|`$, of $`\overline{\rho }(y)`$. In fact, we have verified that the absolute value of the Fourier transformed of this quantity exhibits an abrupt maximum for the first component when transversal gradients begin to build up in the system. The behavior of $`|f_1|`$ as a function of the asymmetry $`\mathrm{\Delta }`$ in the vicinity of the transition point, is shown in Fig. 7 for the same parameters considered in Figs. 2-6. First of all, it must be noted that $`|f_1|`$ varies in a continuous way through the transition, although it grows very fast when one goes into the inhomogeneous region. This is the typical behavior of the order parameter of a nonequilibrium second order phase transition . The continuous character of $`|f_1|`$ introduces some arbitrariness in the determination of the transition point $`\mathrm{\Delta }_c`$. We have made the choice, somewhat arbitrarily but consistently, that the transition takes place when $`|f_1|`$ becomes an order of magnitude larger than its typical value in the reference state, which is determined by the noise level. For the systems used in DSMC, this latter value is of the order of $`10^3`$, so that a system has been considered as transversally inhomogeneous when $`|f_1|10^2`$, which implies deviations from homogeneity of the order of $`1\%`$. This leads in the case of Fig. 6 to an estimation $`\mathrm{\Delta }_c=2.5\pm 0.1`$ for the critical asymmetry.

By changing the initial parameters of the system and repeating the above procedure, $`\mathrm{\Delta }_c`$ has been computed for different values of $`\xi _x^{}`$. The results from DSMC are represented by the filled symbols in Fig. 1. The agreement between the theoretical predictions and the simulation data is rather good, although a systematic deviation appears, larger the smaller $`\xi _x^{}`$. When valuing the comparison, it must be taken into account that the simulations only provide an upper bound for $`\mathrm{\Delta }_c`$. When the system is very close to the transition point, the time required to go from the transversally homogeneous state to the inhomogeneous one may be too large to observe the transition during the simulation time. In any case, it is fair to say that the hydrodynamic equations and the WKB approximation provide an accurate description of what is observed in the simulations.

The results discussed up to this point were obtained from DSMC. Just to illustrate how MD simulations lead to an equivalent scenario, we present next some results for a set of MD simulations with $`\alpha =0.925`$, $`\nu =10^2`$, and $`L_x=100\sigma `$. For these values, it is $`\xi _x^{}=0.281`$. In Fig. 8 a three-dimensional plot of the density profile is shown for $`\mathrm{\Delta }=6.4`$. The system exhibits gradients in the transversal direction, increasing again as we move away from the vibrating wall. However, the density gradients are very small, being the maximum deviation from homogeneity of the order of $`10\%`$. That means that the system is close to the transition, probably in the region where a linear approximation around the steady state still provides an accurate description.

When the asymmetry is increased further on, gradients in the perpendicular direction become sharper, and a state with a very peaked density is observed. This is illustrated in Fig. 9 for $`\mathrm{\Delta }=7`$. In Fig. 10, the quantity $`\overline{\rho }`$ is shown for different simulations belonging to the set we are considering, i.e. they only differ in the value of $`\mathrm{\Delta }`$. The continuous transition from the reference state to the transversally asymmetric one is clearly observed, as well as the dramatic increase of the transversal gradients when the system gets well inside the unstable region. In conclusion, MD results are in full qualitative agreement with the DSMC ones. Even more, the critical values of the asymmetry parameters obtained from MD are in very good quantitative agreement with those following from DSMC, as seen in Fig. 1, where they have been represented by the open symbols. This provides a proof of the validity of the kinetic theory description as given by the Boltzmann equation for dilute inelastic gases.

It is worth to emphasize that points in Fig. 1 were obtained from DSMC and MD simulations by changing the values of $`\alpha `$ and $`\sigma N/L_y`$ to sample different values of $`\xi _x^{}`$. The fact that $`\mathrm{\Delta }_c`$ obtained in this way varies smoothly with $`\xi _x^{}`$ supports the theoretical prediction that the dependence on the different parameters of $`\mathrm{\Delta }_c`$ occurs through it. This has also been ratified by considering two sets of simulations having different values of both $`\alpha `$ and $`\sigma N/L_y`$, but leading to the same value of $`\xi _x^{}`$. The same critical asymmetry parameter was found, supporting the scaling predicted by the theory.

## V Discussion and concluding remarks

In this paper, the spontaneous transversal symmetry breaking in a vibrated granular fluid predicted by Livne et al.  has been further investigated for low density systems. It has been shown that the transition takes places both in MD simulations and in systems described by the Boltzmann equation. Moreover, there is a reasonable good agreement between the theoretical predictions, following from a marginal stability analysis of the hydrodynamic description of the system, and the results from MD and also from the numerical solution of the Boltzmann equation by the DSMC method. This refers to the values of the critical asymmetry as a function of the control parameter, and also to the shape of the hydrodynamic profiles in the vicinity of the symmetry breaking.

To characterize the transition, an order parameter quantifying the initial set up of transversal inhomogeneities has been introduced. In terms of this parameter, the transition presents the features of a second order nonequilibrium phase transition. In this context, it is worth mentioning that no sub-critical bifurcations have been observed in the simulations. Moreover, for states well inside the instability curve, the density profile exhibits a characteristic nonlinear $`\lambda /2`$ shape. On the other hand, in ref. , nonlinear two-dimensional states inside the lineal stability region were found at high densities from the numerical solutions of the hydrodynamic equation. Of course, there is not any contradiction in this, since our analysis has been restricted to low density gases. In ref.  an analytical expression for the marginal stability curve in a certain limit is given. In our notation, the limit considered is $`\xi _x^{}1`$, and the expression reads $`\mathrm{\Delta }1.6/\xi _x^2`$, where we have neglected subleading terms in the density. An asymptotic analysis of the WKB results in this paper, namely of Eq. (45) leads to $`\mathrm{\Delta }_c1.63/\xi _x^{}`$, i.e. a qualitatively different behavior. Given the asymptotic character of the region where these expressions are derived, it is hard to discriminate between the two results from the simulation data.

The work reported here shows once again the generality of spontaneous symmetry breaking phenomena in granular fluids. They occur in isolated as well as in driven granular systems, in dense and dilute flows, with and without gravitational field acting on the particles, i.e., they appear as quite a ubiquitous effect. Interestingly, all the indications point out that they always have a collective origin, which is fully captured by a hydrodynamic description.

Although it may be thought that the spontaneous symmetry breaking discussed in this paper must be closely related with the clustering instability exhibited by a freely evolving granular gas , we believe this relationship deserves more work. The information we have up to now about each of the two phenomena is not the same. We know the mechanism responsible for the clustering instability (the growth of the shear mode relative to the granular temperature), but not the final state attracting the system. On the other hand, the existence of a transversally inhomogeneous steady state has been established for vibrated systems, but the detailed mechanism responsible for the development of instabilities from the reference state is not known.

## VI Acknowledgments

This research was partially supported by the Dirección General de Investigación Cinetífica y Técnica (Spain) through Grant No. PB98-1124.


format: text

text:
# Interlayer tunneling spectroscopy of Bi2Sr2CaCu2O8+δ: a look from inside on the doping phase diagram of high 𝑇_𝑐 superconductors.

## Abstract

A systematic, doping dependent interlayer tunneling spectroscopy of Bi2212 high $`T_c`$ superconductor is presented. An improved resolution made it possible to simultaneously trace the superconducting gap (SG) and the normal state pseudo-gap (PG) in a close vicinity of $`T_c`$ and to analyze closing of the PG at $`T^{}`$. The obtained doping phase diagram exhibits a critical doping point for appearance of the PG and a characteristic crossing of the SG and the PG close to the optimal doping. This points towards coexistence of two different and competing order parameters in Bi2212. Experimental data indicate that the SG can form a combined (large) gap with the PG at $`T<T_c`$ and that the interlayer tunneling becomes progressively incoherent with decreasing doping.

PACS numbers: 74.25.-q, 74.50.+r, 74.72.Hs, 74.80.Dm

Observation of an energy gap in the electronic density of states (DOS) had a decisive role in understanding of low $`T_c`$ superconductivity . However, fifteen years after discovery of high $`T_c`$ superconductors (HTSC), there is still no consensus about HTSC energy gap. Several experiments revealed different energy scales in HTSC . One of those, a normal state pseudo-gap (PG), persists at $`T>T_c`$ . The origin of the PG is an intriguing open question, which is crucial for understanding HTSC. Currently, the scientific community is divided, believing either in superconducting or nonsuperconducting origins of the PG. The resolution can be provided by a doping phase diagram, both because Oxygen-doping is the most critical HTSC parameter (HTSC can be altered from a metal to an antiferromagnetic insulator by decreasing O-content) and because distinctly different diagrams are expected for different scenarios . In a superconducting scenario, the PG represents the pairing energy, which can be finite at $`T>T_c`$ in a strong coupling case. The smaller gap represents the energy required for maintenance of a long range coherence  at $`T<T_c`$. Those two energies should merge in the overdoped (OD) region, as the weak coupling limit is approached. If, on the contrary, the PG appears abruptly at some critical doping point $`p_c`$ and the PG crosses the superconducting gap (SG) at the phase diagram, it would correspond to a non-superconducting PG , which develops in the underdoped (UD) region at the expense of the SG.

The present state of confusion requires further studies using advanced experimental techniques. One of those is an interlayer tunneling spectroscopy, which is unique in it’s ability to measure properties inside HTSC single crystals. This method is specific to strongly anisotropic HTSC, like Bi2Sr2CaCu2O8+δ (Bi2212), in which mobile charge carriers are localized in double CuO2 layers, while the transverse ($`c`$axis) transport is due to interlayer tunneling . Interlayer tunneling has become a powerful tool for studying both electron and phonon DOS of HTSC. It has several important advantages compared to surface tunneling techniques: (i) it probes bulk properties and is insensitive to surface deterioration or surface states; (ii) the current direction is well defined; (iii) the tunnel barrier is atomically perfect and has no extrinsic scattering centers; (iv) mesa structures are mechanically stable and can sustain high bias in a wide range of temperatures ($`T`$) and magnetic fields ($`H`$).

Here I present a systematic O-doping dependent interlayer tunneling study of Bi2212. The spectroscopic resolution was improved by decreasing in-plane mesa sizes, thus avoiding stacking faults and self-heating in the mesas . This way it was possible to trace the SG and the PG at $`TT_c`$ and analyze ”closing” of the PG at a characteristic temperature $`T^{}`$. The obtained doping phase diagram exhibits a critical doping point for appearance of the PG and a characteristic crossing of the SG and the PG close to the optimal doping. This points towards coexistence of two different, competing order parameters in HTSC. Experimental data indicate that the SG can either form a combined gap with the PG or remain uncombined at $`T<T_c`$ and that the interlayer tunneling is predominantly coherent in OD samples, but becomes progressively incoherent with decreasing doping.

Small mesa structures, with areas $`A=1030\mu m^2`$, containing $`N=512`$ intrinsic junctions, were made on top of Bi2212 single crystals by a microfabrication technique . The fabrication was highly reproducible: all mesas on the same crystal exhibited similar behavior, independent of $`A`$ and $`N`$. UD crystals were prepared by annealing in vacuum at 600$`{}_{}{}^{}C`$.

Fig. 1 shows current-voltage characteristics (IVC’s) per junction at 4.2K for different doping. A characteristic knee in IVC’s is clearly seen, followed by a normal resistance branch $`R_N`$. The knee is strongly suppressed both by $`T`$  and $`H`$ , while $`R_N`$ is almost $`T,H`$ independent. Such behavior is typical for SIS-type tunnel junctions, in which the knee occurs at a sum-gap voltage $`2\mathrm{\Delta }_{SG}/e`$, where $`\mathrm{\Delta }_{SG}`$ is the maximum SG. Multiple branches at low bias correspond to one-by-one switching

of junctions from a supercurrent to a quasiparticle (QP) branch. QP branches carry important information: (i) the maximum spacing between QP branches $`\delta V_{QP}`$ is an additional parameter for estimation of the SG; (ii) the extent of QP branches along the vertical axis in Fig. 1 represents the $`I_cR_N`$ product per junction, which is a critical parameters of a Josephson junction.

Figs. 2 and 3 show tunneling conductance $`\sigma =dI/dV`$ curves for slightly OD and UD samples, respectively. Below $`T_c`$ a sharp peak, corresponding to the knee in IVC’s, is seen. The peak voltage, $`V_{peak}`$, decreases as $`TT_c`$. Above $`T_c`$ the peak disappears, but a distinctly different dip-and-hump structure remains, representing the persisting PG . $`T`$ dependencies of the peak (large symbols) and the hump (small symbols + lines) voltages for four samples with different doping are shown in Fig. 4.

For OD samples, $`V_{peak}`$ can be clearly traced up to $`T_c`$ and $`V_{peak}0`$ at $`T_c`$, see Figs. 2 a) and 4. Note that $`V_{peak}(TT_c)`$ is substantially larger than the hump voltage $`V_{hump}(T_c)`$ in the OD mesa. At $`150K`$, $`V_{hump}`$ starts to decrease and vanishes at $`T^{}200K`$ (see Figs. 2 b) and 4). Interestingly, IVC’s are nonlinear even above $`T^{}`$, see Fig. 2 b) and $`\sigma (V)`$ has an inverted parabola shape, which might indicate the presence of van-Hove singularity close to Fermi level in slightly OD samples . Details of the PG closing at $`T^{}`$ are important for understanding the origin of the PG. At the first glance $`V_{hump}(TT^{})`$ resembles a BCS-like dependence, typical for a phase transition due to an onset of charge or spin density waves . However, a different perspective opens when the parabolic background at $`T>T^{}`$ is subtracted, see Fig. 2 c). In such a plot the PG simply ”fills-in” at $`T^{}\mathrm{\Delta }_{PG}`$ without a significant change in $`V_{hump}`$. This may indicate that there is a smooth

crossover rather than a true phase transition at $`T^{}`$.

The behaviour of the SG in UD samples at $`TT_c`$ is one of the most important and yet controversial issues . For UD samples the peak is much weaker than for OD samples even at low $`T`$, cf. Figs. 2 a) and 3 a), and it rapidly smears out with increasing $`T`$. The contrast of the small peak can be increased by subtracting the background PG dip-and-hump at $`T>T_c`$, as shown in Fig. 3 c). Thus $`V_{peak}`$ can be located at $`TT_c`$.

UD samples showed two distinct types of behavior, which I refer to as ”small” and ”large” gap cases, cf. Figs. 1 c) and d). d$`I`$/d$`V`$ curves for the large gap case are shown in Fig. 3. The behavior of large and small gaps is different: (i) for large gaps, the dip-and-hump is strongly enhanced at the expense of the peak, e.g., in Fig. 3 the PG dip-and-hump is clearly recognizable even at low $`T`$. For the large (UD84.4) and small (UD85) gap samples, in Fig. 4, ratios of hump to dip conductances $`\sigma (V_{hump})/\sigma (0)`$ at 100 K are $``$ 5.2 and 1.8, while $`\sigma (V_{peak})R_N`$ at 4.2K is $``$ 1.6 and 10, respectively. (ii) For small gaps $`V_{peak}0`$ at $`T_c`$ and decreases with UD together with $`T_c`$, while for large gaps $`V_{peak}`$ remains finite at $`T_c`$ (even though it drops considerably at $`T_c`$) see Figs. 3 c) and 4, and both peak and hump voltages increase with underdoping despite the decrease of $`T_c`$. (iii) Noticeably, apart from gap magnitudes, other parameters, such as $`\rho _c`$, $`J_c`$ and $`I_cR_N`$ are similar, implying that the tunneling barrier is not affected.

The observed differences can be explained by the following scenarios for formation of small and large gaps,

shown schematically in insets a) and b) of Fig. 4:

The small gap is developed on top of a modest suppression of the DOS at Fermi level, i.e., when there is no true gap at $`T_c`$, which might interfere with the opening SG. Therefore, the peak in d$`I`$/d$`V`$ represents the bare (uncombined) SG, which vanishes at $`T_c`$, while the dip-and-hump represent a ”normal” background, which is hindered by the growing SG. Such behavior was observed for OD, optimally doped and UD samples with the small gap, see Figs. 2 a) and 4.

On the other hand, in the large gap case the SG is developed on top of a true gap $`\mathrm{\Delta }_0`$, see inset b) in Fig. 4. Indeed, from Fig. 3 b) it is seen that the PG dip-and-hump flatten with increasing $`T`$ in a state conserving manner, characteristic for a ”true” energy gap in DOS, and $`\sigma (V)`$ curves intersect in one point, indicating approximately constant value of the PG in the measured $`T`$ range. Below $`T_c`$ this causes formation of the combined ($`\mathrm{\Delta }_0`$ and $`\mathrm{\Delta }_{SG}`$) large gap. In agreement with this assumption: (i) the large gap does not vanish, but approaches $`\mathrm{\Delta }_0`$ at $`T_c`$, see Figs. 3a) and 4. (ii) The peak completely disappears at $`T_c`$ but does not transform into the hump because $`eV_{hump}>\mathrm{\Delta }_0`$, see the dashed line in inset b). (iii) The volume of the peak (superfluid density) is small because it builds up from an initially suppressed DOS. (iv) The opening of the SG at $`T<T_c`$ shifts all DOS features, including the hump, as shown in inset b)

of Fig. 4. The correlated shift of both the peak and the hump with $`T`$ for UD84.4 sample, as shown in Figs. 3 and 4, is a strong argument in favor of the combined scenario of the large gap. Similarly, uncorrelated $`T`$ dependent peak and $`T`$ independent hump in the OD93 sample, see Figs. 2 and 4, suggests that the small gap represents the uncombined SG. Interestingly, if we take $`V_{peak}(4.2K)V_{peak}(T_c)`$ as a measure of the SG part of the combined gap, it will coincide with the small gap for a similar doping, as shown by arrows in Fig. 5 d). A systematic increase of the hump energy with decreasing $`T`$, observed by ARPES , would have been consistent with the combined scenario of the large gap if not for the lack of correlated $`T`$dependence of the coherence peak.

Fig. 5 shows O-doping dependencies of: a) $`T_c`$, dashed line represents the empirical expression, used for estimation of $`p`$; b) the critical current density, $`J_c`$, and the $`I_cR_N`$ product per junction; c) the tunneling resistivity at large bias $`\rho _c=R_NA/(Ns)`$. The $`I_cR_N`$ is an important parameters of a Josephson junction. As Bi2212 is likely to be a d-wave superconductor, the $`I_cR_N`$ depends both on $`\mathrm{\Delta }_{SG}`$ and the coherence (in-plane momentum conservation) of $`c`$axis tunneling (another highly debated issue in HTSC). The $`I_cR_N`$ is maximum $`\mathrm{\Delta }_{SG}/e`$ for coherent, and zero for completely incoherent tunneling. For OD mesas $`I_cR_N10`$mV is a considerable fraction $`0.6`$ of $`\mathrm{\Delta }_{SG}/e`$, indicating predominantly coherent nature of the interlayer tunneling. With underdoping,

the $`I_cR_N`$ decreases dramatically at a much faster rate than $`\mathrm{\Delta }_{SG}`$. This indicates that the interlayer tunneling becomes progressively incoherent in UD Bi2212.

Fig. 5 d) shows the obtained doping phase diagram of Bi2212. Here I plot $`\frac{1}{2}V_{peak}(4.2K)`$ $`\mathrm{\Delta }_{SG}/e`$, $`\frac{1}{2}V_{hump}(100K>T_c)`$ $`\mathrm{\Delta }_{PG}/e`$ and $`\delta V_{QP}(4.2K)`$, It is seen that the small gap ($``$ SG) shows a similar tendency as $`T_c`$ and decreases both on OD and UD sides. This is also supported by a correlated behavior of $`\delta V_{QP}`$. In contrast, the large gap ($``$ PG) increases approximately linearly with underdoping, as shown by the solid line. The PG and the SG lines cross at about the optimal doping, $`p=0.16`$. On the OD side the PG becomes considerably less than the SG and shows a clear tendency to vanish at the critical doping point, $`p_c0.19`$. This speaks in favor of a nonsuperconducting origin of the PG , consistent with earlier observations of different $`T`$  and $`H`$  dependencies of the SG and the PG. Within such a scenario, a suppression of superconductivity (decrease of $`T_c`$, $`\mathrm{\Delta }_{SG}`$, the superfluid density, etc.) in UD HTSC is caused by appearance of the competing order parameter (PG), e.g., due to strengthening of antiferromagnetic correlations and formation of spin density waves. Note that a similar phase diagram, attributed to competition between superconducting and antiferromagnetic orders, was reported for heavy fermion superconductors.

At present, the reason for appearance of either small or large gaps in UD samples is unclear. However, it is not due to irreproducibility of fabrication (all mesas on the same crystal show the same behavior) or macroscopic defects (regular QP branches are observed in both cases). Presumably, the ambiguity is connected with a microscopic inhomoginiety of UD crystals. The presence of ambiguity obscures identification of the genuine HTSC behavior in the UD region. However, there is no ambiguity for overdoped and optimally doped samples. Therefore, conclusions that there is a critical doping point in HTSC phase diagram and that the SG and the PG cross rather than merge near the optimal doping are robust.

In summary, O-doping dependence of Bi2212 was studied using high resolution interlayer tunneling spectroscopy. We were able to simultaneously trace the superconducting gap and the $`c`$axis pseudo-gap at $`TT_c`$ and analyze ”closing” of the PG at $`T^{}`$. The obtained doping phase diagram exhibits a critical doping point for appearance of the PG and a characteristic crossing of the SG and the PG close to the optimal doping, indicating a competing nature of two coexisting order parameters in HTSC. In UD samples, the SG can either form a combined gap with the PG or remain uncombined at $`T<T_c`$, but the bare SG vanishes at $`TT_c`$ for all studied doping levels. Analysis of $`I_cR_N`$ vs. $`\mathrm{\Delta }_{SG}`$ indicates that the interlayer tunneling is predominantly coherent in OD, but becomes progressively incoherent in UD samples.


text:
# Electron-Transport Properties of Na Nanowires under Applied Bias Voltages

## Abstract

We present first-principles calculations on electron transport through Na nanowires at finite bias voltages. The nanowire exhibits a nonlinear current-voltage characteristic and negative differential conductance. The latter is explained by the drastic suppression of the transmission peaks which is attributed to the electron transportability of the negatively biased plinth attached to the end of the nanowire. In addition, the finding that a voltage drop preferentially occurs on the negatively biased side of the nanowire is discussed in relation to the electronic structure and conduction.

preprint: APS/123-QED

Currently, quantized conductances of atomic nanowires with a few atoms width and nanometer length, which are in steps of $`G_0=2e^2/h`$ ($`e`$: the electron charge, $`h`$: Planck’s constant), have attracted much attention because of both fundamental and technological importance. A number of exciting experimental and theoretical studies on the conductance quantization for various metal nanowires have been reportedJMRuitenbeek00 ; JMKrans95 ; HOhnishi98 ; STsukamoto01 ; NKobayashi00 ; ANakamura99 ; NDLang97 ; NDLang95 ; RNBarnett97 ; JCCuevas98 ; HSSim01 . Above all, the studies for Na nanowires are of particular interest, because the nature of electron conduction through nanowires is expected to emerge best for Na nanowires having the simplest electronic structure where monovalent electrons predominate. Elaborate investigations on electron-transport properties of Na nanowires are, however, not many yet. A few experiments have been reportedJMKrans95 , where the variation of the conductance was measured during the withdrawal of the Na film to show the stairlike behavior of the quantized conductance. Theoretically, there have been some first-principles calculations fully self-consistently performed using the model incorporating a couple of semi-infinite electrodes so as to evaluate electron flow through Na nanowiresSTsukamoto01 ; NKobayashi00 ; ANakamura99 ; NDLang97 . But, these theoretical studies are limited to zero bias voltageANakamura99 ; NKobayashi00 ; NDLang97  or the low bias regimeSTsukamoto01 . The analysis and discussion of electron transport through the Na nanowire should be extended to a higher bias voltage in order to gain insight into a unique conduction property possible under a strong electric field.

In this letter, we present a theoretical analysis of electron transport through the Na nanowire under the application of a sizable bias voltage, based on the density functional theory within the local density approximation. All the calculations determining the electronic structure and transport properties are self-consistently evaluated by the method making use of the Lippmann-Schwinger equation proposed by LangNDLang95 . The norm-conserving pseudopotentialsNTroullier91  are adopted to describe the electron-ion interaction, and the real-space timesaving double-grid techniqueTOno99  is used by virtue of an efficient computation of the nonlocal parts of the pseudopotentials. We will demonstrate that the Na nanowire exhibits an unusual nonlinear current-voltage (I-V) characteristic which leads to negative differential conductance behavior. Notice that the negative differential conductance is expected to be an essential property to realize fast switching in future electronic devices. This paper is probably the first report describing that a simple nanowire in which all atoms are tightly binding together possesses the negative feature of the differential conductancecomment2 . We will also discuss how an applied bias voltage gives rise to the change of the effective potential along the nanowire (voltage drop) and of the charge density distribution.

According to Refs. NKobayashi00,  and STsukamoto01, , we here adopt the model composed of three parts for the nanowire system: a straight atomic nanowire, two plinths, and a couple of semi-infinite jellium electrodes sandwiching them (see Fig. 6 for illustration). The straight nanowire is made up of $`N`$ ($`N`$=1$``$5) Na atoms equispaced at 7.0 a.u., which is equal to the nearest neighbor atomic distance of the Na crystal. Each plinth is modeled after the Na(001) ideal surface, i.e., it is formed by four Na atoms, one at each corner of a square with the side length of 8.1 a.u. These plinths are attached to the two ends of the straight nanowire and act as buffer layers to moderate the unphysical effect due to an oversimplified electrode model of jellium. The electron density of the jellium electrodes is taken equal to the average valence-electron density of the Na crystal ($`r_s=4.0`$ a.u.)NKobayashi00 . The distance between the jellium edge and the atoms forming the plinth is set to be 4.0 a.u. In the present calculation, we do not consider any structural relaxation.

In the numerical calculations, we employ the supercell on which periodic boundary conditions are imposed in the directions perpendicular to the nanowire ($`x`$ and $`y`$ directions) and nonperiodic boundary condition in the direction parallel to it ($`z`$ direction). Therefore, wave functions are expressed by the Laue representation. The side length of the supercell in the $`x`$ and $`y`$ directions is 24.0 a.u. and the number of the plane waves in these directions is taken to be $`24^2`$, which corresponds to a cutoff energy of 7.4 Ry. We use only the $`\mathrm{\Gamma }`$ point in the $`xy`$ Brillouin zone because of a large supercell adopted here. The grid size in the $`z`$ direction is chosen to be 0.5 a.u.

Before moving on to the discussion of electron-transport properties for the Na nanowire system at a finite bias voltage, we briefly report the results obtained at the zero bias limit. In Fig. 1(a), the calculated conductance is drawn as a function of the number of atoms $`N`$ in the straight nanowire: in the case of $`N=1`$, i.e., the case of the closed packed structure of Na atoms, the conductance takes the value $`2.8G_0`$, while in the cases of $`N=`$2$``$5, it is close to the quantized value of $`1G_0`$NKobayashi00 . These features of the conductance are in good agreement with the experimental results obtained by Krans et al., where the conductance histogram shows the three distinct points: a peak at slightly lower than the $`3G_0`$ quantized value, a peak at $`1G_0`$ and a missing $`2G_0`$ peakJMKrans95 . In Fig. 1(a), an oscillatory behavior in the conductance is observed for $`N=25`$, i.e., the even-atom-number nanowire systems exhibit slightly lower conductance than the quantized value of $`1G_0`$, and the odd-atom-number systems show fully quantized conductance $`1G_0`$. This oscillation can be understood in terms of the local density of states (LDOS) around the nanowire at the Fermi level: as seen in Fig. 1(b), the LDOS for the even-atom-number (odd-atom-number) nanowire has an off-resonance (on-resonance) transport property at the Fermi levelNDLang98 . Based on the Friedel sum rule, Sim et al. recently predicted the occurrence of the even-odd oscillation for the Na nanowire sandwiched by the pyramidal-shape tips HSSim01 . Our present results verify their prediction of the even-odd oscillation.

We now discuss how electron-transport properties of the nanowire are affected when a sizable bias voltage is applied. The range of bias voltages is 0.0$``$3.5V and a long nanowire of $`N=5`$ is adopted. The current per supercell is simply expressed as the integration of the transmission over the voltage window, which represents the energy range between the Fermi levels of the two electrodes. Hereafter, the respective electrodes are called left and right electrodes, and the bias voltage is applied with the left electrode negative and the Fermi level of the right electrode kept fixed. Figure 2 shows the I-V characteristic of the system and its differential conductance. Here, the differential conductance is defined as the derivative of the current with respect to the bias voltage. It is clearly seen that the system exhibits the nonlinear I-V characteristic and in particular, the differential conductance shows unusual behavior near zero bias voltage and around 2.5 V; the differential conductance rapidly decreases in the range of 0.0$``$0.5 V, deviating from the quantized value $`1G_0`$ observed at the zero bias limit, and moreover, it is negative at the bias voltages of 2.3$``$2.7 V with a minimum at 2.5 V.

We first examine why the reduction of the differential conductance is induced at the low bias voltages of 0.0$``$0.5 V. The transmission and the LDOS around the straight nanowire of electrons incident from the left electrode are depicted as a function of the electron energy measured from the Fermi level of the right electrode in Figs. 3(a) and (b), respectively. The right edges of the voltage windows (shaded areas) in these figures correspond to the Fermi level of the left electrode for the respective bias voltages. When the bias voltage is boosted from 0.0 V to 0.5 V, a dip in the transmission curve, which is indicated by D in Fig. 3(a), appears in the voltage window, causing the reduction in the differential conductance. This transmission dip D is brought about by the off-peak of the LDOS marked D in Fig. 3(b)

Next, we discuss the negative feature of the differential conductance at around 2.5 V seen in Fig. 2. In the transmission in Fig. 3(a), one can see that the leftmost peak clearly observed at the bias voltage of 2.0 V is suddenly suppressed at 2.5 V (A$``$A’), and that one more peak at 2.5 V is also suppressed when the bias voltage is raised to 3.0 V (B$``$B’)comment1 . As a consequence of the suppression, the integration of the transmission over the voltage window decreases upon raising the bias voltage from 2.0 V to 3.0 V, which leads to the negative differential conductance. In contrast, such a disappearance of the transmission peak is no longer observed at the bias voltage of 3.5 V. It should be noticed that the suppression of the transmission peaks is in accordance with the vanishing of the corresponding peaks in the LDOS around the straight nanowire of electrons incident from the left electrode, marked by A$``$A’ and B$``$B’ in Fig. 3(b). This vanishing implies that some electrons coming from the left electrode with specific energies cannot reach the nanowire when the bias voltage exceeds a certain threshold.

In order to investigate why these transmission peaks and LDOS peaks around the straight nanowire suddenly disappear as the bias voltage is raised, we focus on the electronic structure around the left plinth. Figure 4 shows the LDOS around the left plinth of electrons incident from the left electrode at the bias voltages of 0.0, 2.0, 2.5 and 3.0 V, where each curve shifts to the higher energy side according to the applied bias voltage. This shift corresponds to the linear rising of the effective potential around the left plinth upon applying bias voltage, which will be discussed later for Fig. 5. We see in Fig. 4 that the leftmost LDOS peaks appear at 0.0, 0.5 and 1.0 eV at the bias voltages of 2.0, 2.5 and 3.0 V, respectively. It is now common knowledge that the energy of the leftmost LDOS peak corresponds to the first onset energy of the transmissionNKobayashi00 ; STsukamoto01 ; JCCuevas98 . Consequently, when the bias voltage of 2.0 V is applied, electrons with energy lower than 0.0 eV are refused to go through the left plinth, while those with energy above 0.0 eV readily enter the plinth and the nanowire to yield the transmission peak A in Fig. 3(a). In the case of the bias voltage of 2.5 V, the threshold energy for electrons to pass through the left plinth is 0.5 eV, therefore, electrons with the energy of 0.0 eV are prohibited from entering the left plinth. Thus, as the bias voltage is boosted from 2.0 V to 2.5 V, the suppression of the transmission peak at 0.0 eV, A$``$A’, is brought about. This mechanism of the suppression of the transmission peak also holds in the case of B$``$B’ in Fig. 3(a).

Finally, we explore how the applied bias voltages drop along the nanowire. Figure 5 shows the differences of the local effective potential along the axis of the nanowire at some significant bias voltages from that at zero bias voltage, $`V_{eff}(V_B)V_{eff}(0.0)`$, at $`V_B=`$2.0, 2.5, 3.0 and 3.5 V. In Fig. 5, one can see that the intensive voltage drop takes place mainly around the negatively biased end of the nanowire (atoms numbered 1 and 2 in Fig. 5)MBrandbyge99 , while the voltage drop is nominal near the two plinths and on the positively biased side (atoms 3$``$5). The localization of the voltage drop is understood as follows. Electrons abundant in proportion to the applied bias voltage rush into the left end of the nanowire from the left electrode through the plinth. The congestion of electron flow occurs around the region between atoms 1 and 2, where the narrowing of the path along which electrons pass takes place due to the geometrical transition from the wider electrode to the narrower nanowire. Although only electrons with some specific momenta can enter the conduction channel and reach the opposite electrode without any backscattering, the remaining electrons suffer from a considerable amount of reflection. Thus, the resistance for electron flow localizes around the region between atoms 1 and 2, which results in the intensive voltage drop there. This mechanism of electron flow is analogous to the vena contracta effect known in hydrodynamics. The absence of a sizable voltage drop near the two plinths and on the positively biased side of the nanowire is associated with the change in the charge density distribution upon applying bias voltage. Figure 6 shows the profiles of the difference of the charge densities, $`\rho (V_B)\rho (0.0)`$, at $`V_B=`$2.5 and 3.5 V. On the right of the atoms forming the left plinth and the left of atom 2, one observes considerable increases of the charge density. These increases reflect the screening effect which keeps the electric potential of the negatively biased jellium electrode and plinth equal to the negative electrode potential applied to the left jellium. On the other hand, the most significant decreases of electrons occur on the right of atom 2 so that the electric potential of the positively biased side of the nanowire system (atoms 3$``$5, the right plinth and jellium electrode) becomes close to the positive electrode potential assigned to the right jellium. As a consequence of the negative electric potential on the left of atom 1 and the positive one on the right of atom 2, the localization of the voltage drop between atoms 1 and 2 is naturally realized.

In conclusion, a first-principles calculation method using the Lippmann-Schwinger equation was implemented to elucidate electron-transport properties of Na nanowires at finite bias voltages. Employing a nanowire model including the plinths as well as the semi-infinite jellium electrodes, we obtained the nonlinear I-V characteristic and the unusual behavior of the differential conductance. We demonstrated that the negative feature of the differential conductance is explained by the suppression of the transmission peaks due to the electron transportability of the negatively biased plinth. It was also found that the voltage drop preferentially occurs around the negatively biased end of the nanowire, and this finding was discussed in relation to the electronic structure and conduction of the nanowire.

This work was supported by a Grand-in-Aid for COE Reseach (No. 08CE2004) from the Ministry of Education, Culture, Sports, Science and Technology.


text:
# Secular Evolution in Barred Galaxies: Observations

## 1. Introduction

I had originally intended to give an overview of the observational evidence relating to the mechanisms by which bars evolve over time. However, this overview has been expertly provided by a number of other authors in these proceedings. Bureau (2002), for example, gives a thorough treatment of the buckling of bars to form boxy bulges, and how this may lead to their dissolution. Similarly, Das et al. (2002) discuss the role that central mass concentrations, such as massive black holes, may play in weakening and ultimately destroying bars. Debattista (2002) describes the way that bars might be expected to slow down through dynamical friction, while Gerssen (2002) presents the observational evidence that not much slowing has occurred. This apparent contradiction suggests that the mass distribution in barred galaxies cannot be very centrally concentrated, or that some other mechanism destroys the bars before they have had the opportunity to slow down. Finally, Athanassoula (2002), through her sophisticated N-body simulations, discusses much of the theory behind these observations. This wide sweep of articles leaves me very little to say by way of background review.

I will therefore limit this paper to discussing how these observations, and the accompanying theory, might be brought together in a common quantitative framework, and to presenting a couple of preliminary applications of this procedure. Section 2. discusses the quantification of bar strength, and how it can be combined with other quantitative measures of morphology to define a parameter space directly analogous to that used in the previous qualitative classification of Hubble’s tuning fork. Section 3. describes an application of this parameterization in studying the evolution of bar pattern speeds, and Section 4. looks at direct evidence for bar evolution in the Hubble Deep Fields. Finally, Section 5. speculates on the next steps in this field.

## 2. Quantifying Bar Strength and Morphology

Historically, morphological features such as the presence or absence of a bar in a galaxy have been assessed qualitatively by inspection of photographic plates \[e.g. Sandage & Tammann 1981 (RSA), de Vaucouleurs et al. 1991 (RC3)\]. Although this approach has a strong romantic appeal, it has a number of shortcomings in modern studies of galaxy morphology. First, such human inspection unavoidably has a subjective element, which makes comparison between different analyses rather difficult. Second, projects such as the Sloan Digital Sky Survey are now becoming so large that it is not practical to inspect each galaxy image individually. Third, such classifications do not make full use of the quantitative data that is now available from CCD images. Although this under-utilization may not matter when classifying high-quality images of nearby galaxies, it becomes vital when analyzing small noisy images of very distant galaxies, where every photon counts. Finally, an objective quantitative measure of a property such as bar strength is vital if one wishes to correlate galaxies’ “barriness” with their other properties to seek insight into the way that bars form and evolve.

A number of methods have been advocated for measuring bar strength in a quantitative manner. The simplest involves just fitting isophotes to a galaxy at different radii (e.g. Martin 1995). The intrinsic shape of the galaxy’s isophotes can then be measured by deprojecting it to face on. This deprojection requires that one assume that the galaxy can be treated as a thin disk, and that its inclination can be estimated from the apparent shape of its outer isophotes, which are taken to be intrinsically round. A more robust variant on this technique (Abraham et al. 1999) involves measuring the apparent ellipticity as a function of radius by calculating the second moments of the light distribution as a function of surface brightness threshold. Whichever approach is used, some measure of the bar flattening, such as its intrinsic ellipticity or axis ratio, then provides the measure of bar strength. A more physically-motivated approach has been advocated by Block et al. (2002), who seek to measure the gravitational influence of a galaxy’s non-axisymmetric component relative to the axisymmetric forces. Each of these measures is optimized for different purposes: for example, the moment-based technique is designed to deal with faint distant galaxies where the signal-to-noise ratio is low, whereas the gravitational force method is intended to tie the observations as closely as possible to the underlying physics where higher quality data exist. All the methods have to make some assumption about the shapes of the galaxies surveyed in order to derive their three-dimensional properties from two-dimensional images. In many cases, these assumptions will be violated – treating a galaxy as a thin disk when it contains a significant bulge component is clearly incorrect. However, even where the assumptions are invalid, the technique will still provide a measure of barriness. As long as the measure is compared to the same quantity derived from other samples that have been treated in the same way, one can still draw meaningful conclusions from any differences between the statistics from the various samples.

It seems very likely that bars cannot be fully described by any such single parameter – for example, is a small highly-distorted central isophote a stronger or weaker bar than a much larger slightly-distorted central feature? Ultimately, one probably needs at least two parameters to capture the essence of a bar, one describing the degree to which it distorts the galaxy’s isophotes, while the other specifies the fraction of the galaxy’s light that lies in the bar. However, even a single parameter does provide an objective and repeatable measure of barriness which for the first time allows a quantitative discussion of bar properties.

One simple use of such a measure of bar strength has been in an investigation of whether Hubble’s (1936) tuning fork classification of galaxies represents a physical phenomenon, or whether it is just an aesthetic idealization. Hubble’s scheme arranges the galaxies from early type (big bulge, weak tightly-wound spiral structure) to late type (small bulge, strong open spiral structure) along a horizontal axis, bifurcating vertically into barred and unbarred galaxies. Thus, it can be thought of as a two-dimensional parameter space, with a measure of lateness along the $`x`$ axis and a measure of barriness on the $`y`$ axis. Figure 1 shows a quantitative realization of this “Hubble Space,” where Abraham et al.’s (1994) central concentration index, $`C`$, is used to measure the position of a galaxy along the Hubble sequence, and Abraham & Merrifield’s (2001) measure of bar strength, $`f_{\mathrm{bar}}`$, delineates the $`y`$ axis. The picture is a little complicated here, partly because the high values of $`C`$ correspond to galaxies with large central concentrations of light, which translate into the big bulges of early-type galaxies, so the axis runs the opposite way to Hubble’s traditional early-to-late sequence. Further, the parameters are somewhat correlated: the appearance of a bar will be weakened in a galaxy with a strong central bulge, whereas even a rather small bar shows up plainly in a late-type galaxy with a weak bulge, so there is a general trend from top left to bottom right in Fig. 1. However, despite these complications, there does appear to be at least some sign of Hubble’s tuning fork in this picture, as the dotted line in the figure indicates. The barred galaxies, in particular, seem to form a remarkably tight sequence, well separated from the main group of unbarred systems. A simple statistical test suggests that the deficit of galaxies around values of $`f_{\mathrm{bar}}0.2`$, $`C0.4`$ is significant at the 95% level, although such a posteriori calculations should not be over-interpreted.

If the gap between barred and unbarred galaxies is real, it has important consequences for the secular evolution of bars. It could mean that barred and unbarred galaxies are fundamentally distinct entities, with no evolution between the two types of system. However, these systems are remarkably similar in all their other properties. For example, one might have thought that the presence or absence of a bar could be related to how strongly the galaxy is dominated by dark matter, which will alter the effectiveness of the classical bar instability (Ostriker & Peebles 1973). The Tully-Fisher relation essentially plots the correlation between luminous mass as quantified by the absolute magnitude against total mass as quantified by the rotation speed, so one might expect the relation to be different for barred and unbarred galaxies. However, the mean Tully Fisher relation is identical for barred and unbarred systems (Debattista & Sellwood 2000). There is a somewhat larger scatter in the relation for barred galaxies, but this can be attributed to the difficulty in measuring the circular rotation speed of a barred galaxy where the stars and gas do not follow circular orbits (Franx & de Zeuuw 1992). Given that barred and unbarred galaxies seem otherwise identical, and that there are the various credible mechanisms mentioned in the Introduction by which bars might be destroyed, it seems unlikely that barred and unbarred galaxies are fundamentally distinct. More plausibly, the absence of “semi-barred” galaxies can be taken to show that galaxies evolve rapidly between the two extreme states – the absence of intermediate galaxies is analogous to the Hertzsprung gap in the color-magnitude diagram for stars, which represents a region through which stars evolve rapidly.

This result conflicts with the RC3 classifications, which give a large number of galaxies the “SAB” designation, indicating that they have intermediate bar strength. This conflict can be resolved by looking at where the SAB galaxies are located in Fig. 1. They are found to lie scattered around $`C0.3`$, $`f_{\mathrm{bar}}0.15`$; thus, they are the late-type galaxies on the unbarred branch of Hubble’s tuning fork. Their values of $`f_{\mathrm{bar}}`$ are, indeed, intermediate between those of the unbarred and strongly barred galaxies, but this is simple because of the tilt in the tuning fork. In essence, the SAB systems are unbarred galaxies that reveal minor bar-like features at their centers; these fairly insignificant features would have been hidden in galaxies with bigger bulges.

One concern with the result presented in Fig. 1 is that the morphologies have been derived from optical imaging, where dust and localized star formation might disguise the underlying galaxy morphologies (Block et al. 1999). Further, it is not clear how the galaxies in the Frei Catalog were selected. It is therefore possible that the bifurcation in Fig. 1 represents the natural tendency to pick the extremes of a continuous distribution to represent the archetypes of barred and unbarred systems. To address these concerns, Whyte et al. (2002) have analyzed H-band images from the Ohio State University Bright Spiral Galaxy Survey. These galaxies form a relatively complete sub-sample from the RC3, so should not have any bias toward artificial bimodality. Nonetheless, as Figure 2 shows, there is still a significant dip in the distribution of bar strengths at $`f_{\mathrm{bar}}0.2`$. Since this plot does not have the smearing effect of the tilt in the tuning fork taken out, the true bimodality of the distribution is likely to be even stronger. Thus, the picture remains that any evolution from barred to unbarred or vice-versa must happen on a rapid timescale compared to the ages of the galaxies.

## 3. Evolution in Morphology and Pattern Speed

A second application of this quantitative approach to morphology is provided by studies of bar pattern speeds. A variety of techniques have now been employed for measuring bar pattern speeds (e.g. Gerssen 2002). The results all indicate that the bar patterns are rotating at such a rate that their ends lie close to the point at which the pattern rotates as rapidly as the galaxy, the co-rotation radius $`r_{\mathrm{cr}}`$. A bar can be conveniently parameterized by the co-rotation radius divided by the length of its major axis, $`=r_{\mathrm{cr}}/a`$, so the observations imply that $`1`$ for all galaxies, which is the smallest value possible for any self-consistent bar pattern. This result is somewhat surprising because, as described in Debattista (2002), a mechanism exists by which bars evolve to ever lower pattern speeds until they essentially grind to a halt, so one might expect to find bars with a wide range of pattern speeds and values of $``$ much greater than unity.

Some clue as to what is going on can be gleaned by looking at the distribution in Hubble Space of galaxies with measured bar pattern speed. As Figure 3 shows, these galaxies describe the expected fairly tight sequence which matches that already derived for barred galaxies. However, there are indications that the pattern speeds vary systematically across the width of the sequence, with somewhat slower pattern speeds being associated with weaker bars. This result would suggest that barred galaxies follow evolutionary tracks that lead them downward in this figure as their bars start to slow. The absence of any slower bars suggest that any such galaxies have entered the gap region where systems evolve rapidly via some other mechanism that destroys the bar.

Clearly, however, larger samples of bars with measured pattern speeds are required if this result is to be confirmed. It will also be highly instructive to calculate the Hubble Space parameters for the bars that form and slow down in N-body simulations (Debattista 2002) in order to follow their evolutionary tracks in this space, to see how well they match up to the observations.

## 4. Direct Observations of Bar Evolution

With the more robust measures of bar strength, such as the moment-based techniques, it is possible to characterize this property even in small faint galaxies, such as those imaged in the Hubble Deep Fields (HDFs). This possibility is particularly intriguing as it allows one to investigate the secular evolution of barred galaxies over a cosmological timescale by examining the properties of bars in the population as a function of redshift.

An analysis of the HDFs, summarized in Fig. 4, led Abraham et al. (1999) to conclude that there has been strong evolution in the barred galaxy population to the extent that there are almost no barred systems beyond a redshift of $`z0.7`$. This result quantified the qualitative impression that van den Bergh et al. (1996) had obtained from visual inspection of the first HDF data. This observation lends itself to the simple physical interpretation that at these early epochs disks were either too dynamically hot or too low in mass to undergo the conventional bar instability.

There are, however, several possible snags with this analysis. The first issue is bandshifting: all the data analyzed came from I-band images, but at significant redshifts the observed I-band light will have been emitted further to the blue, so one is not comparing like with like when looking at galaxies at different redshifts. However, the effects here are fairly modest: at the redshifts where there is an absence of bars, the I-band images correspond to rest-frame B-band emission. Since B-band images of nearby galaxies reveal plenty of barred systems, the absence of bars beyond $`z0.7`$ cannot simply be attributed to the observations being made in waveband where bars are not detectable.

The second problem is that the more distant galaxies tend to have lower signal-to-noise ratios, which could prevent their bars being detected. One subtler effect is that cosmological dimming may reduce the surface brightness of a galaxy’s disk to an undetectable level, so a face-on barred galaxy at high redshift may be misidentified as an edge-on system. Since edge-on systems cannot have their bar strengths measured, they are excluded from the analysis, which might explain where the barred systems have gone. However, Abraham et al.’s (1999) experiments involving artificially redshifting and degrading the images of the closer galaxies showed that their bars would have remained detectable even at the highest redshifts observed, and that the disks remain bright enough for there to be no confusion over inclination.

It is still possible that barred galaxies exist at high redshifts: if, for example, these systems have much weaker disks than those in the nearby Universe, then they could have been misidentified as edge-on systems. One is, however, forced to the conclusion that if there are barred galaxies beyond $`z0.7`$, then they are sufficiently different in structure from those in the nearby Universe that they cannot be detected in the HDF data.

## 5. Conclusions and Future Work

With the advent of largescale CCD surveys of galaxies, both nearby and in the more distant Universe, the study of galaxy morphology has matured from the qualitative to the quantitative. Robust techniques have been developed for automatically measuring properties of galaxies such as their bulge-to-disk ratios and bar strengths. Perhaps unsurprisingly, these quantitative measures correspond quite closely to their qualitative predecessors, to the extent that one can reproduce Hubble’s tuning fork classification. The bifurcation between barred and unbarred galaxies has significant implications for investigations of secular evolution in bars, since it implies that any evolution between these two end states must occur on a timescale that is short compared to the galaxies’ cosmological lifetimes. Further evidence for evolution in bar strength on timescales that are less than the cosmological have come from studies of faint galaxies in the HDFs, which indicate that bar formation is a relatively recent phenomenon.

An important next step will be to use this framework to compare the observations with the predictions of N-body simulations. It is now possible to calculate morphological parameters to follow the properties of N-body simulations, such as those described by Athanassoula (2002), as the systems undergo secular evolution. Further, cosmological simulations like those presented by Navarro (2002) are rapidly advancing to the point where they can reliably follow the formation and evolution of individual disk galaxies, so we can also follow the variations in morphological parameters in a cosmological framework. Through such analysis, galaxy observations such as those presented here can be directly tied to N-body simulations to test our understanding of the formation and evolution of these systems, just as the observed sequences of stars in the color-magnitude diagram can be combined with stellar evolutionary tracks in order to understand how stars form and evolve.

One other area where there is still much room for analysis lies in the morphological description of bars. As described in Section 2., it seems unlikely that the essence of a bar can be captured in a single parameter, but it is by no means clear what combination of parameters does provide the optimal description, nor what the physical differences are that give rise to the range of parameters observed. A fuller description of morphological features such as bars will enable us to bring more sophistication to our exploration of the mechanisms by which these features evolve, and hopefully answer many of the outstanding questions.

## References

Abraham, R.G. & Merrifield, M.R. 2000, AJ, 120, 2835

Abraham, R.G., Merrifield, M.R., Ellis, R.S., Tanvir, N.R. & Brinchmann, J. 1999, MNRAS, 308, 569

Abraham, R.G., Valdes, F., Yee, H.K.C. & van den Bergh, S. 1994, ApJ, 432, 75

Athanassoula, L. 2002, this volume

Block, D.L., Puerari, I., Frogel, J.A., Eskridge, P.B., Stockton, A. & Fuchs, B. 1999, Ap&SS, 269, 423

Block, D.L., Puerari, I., Knapen, J.H., Elmegreen, B.G., Buta, R., Stedman, S., Elmegreen, D.M. 2002, this volume

Bureau, M. 2002, this volume

Das, M., Teuben, P.J., Vogel, S.N., Regan, M.W., Sheth, K. & Harris, A. 2002, this volume

Debattista, V. 2002, this volume

Debattista, V. & Sellwood, J. 2000, ApJ, 543, 704

de Vaucouleurs, G., de Vaucouleurs, A., Corwin, H.G., Buta, R.J., Paturel, G. & Fouque, P. 1991, Third Reference Catalog of Bright Galaxies (New York: Springer-Verlag)

Elmegreen, B. 1996, in ASP Conf. Ser. Vol. 91, Barred Galaxies, ed. R. Buta, D.A. Crocker & B.G. Elmegreen (San Francisco: ASP), 197

Franx, M. & de Zeuuw, T. 1992, ApJ, 392, L47

Frei, Z., Guhathakurta, P., Gunn, J.E. & Tyson, J.A. 1996, AJ, 111, 174

Gerssen, J. 2002, this volume

Hubble, E.P. 1936, The Realm of the Nebulae (New Haven: Yale University Press)

Martin, P. 1995, AJ, 109, 2428

Navarro, J. 2002, this volume

Ostriker, J.P. & Peebles, P.J.E. 1973, ApJ, 186, 467

Sandage, A. & Tammann, G.A. 1981, Revised Shapley-Ames Catalog of Bright Galaxies (Washington: Carnegie Institute)

van den Bergh S., Abraham R.G., Ellis R.S., Tanvir N.R. & Santiago B.X. 1996, AJ, 112, 359

Whyte, L.F., Abraham, R.G., Merrifield, M.R., Eskridge, P.B. & Frogel, J.A. 2002, in preparation


text: # Untitled Document