text:
# On the Stability of Knots in Excitable Media

\[

To appear in Phys. Rev. E.

## Abstract

Through extensive numerical simulations we investigate the evolution of knotted and linked vortices in the FitzHugh-Nagumo model. On medium time scales, of the order of a hundred times the vortex rotation period, knots simultaneously translate and precess with very little change of shape. However, on long time scales we find that knots evolve in a more complicated manner, with particular arcs expanding and contracting, producing substantial variations in the total length. The topology of a knot is preserved during the evolution and after several thousand vortex rotation periods the knot appears to approach an asymptotic state. Furthermore, this asymptotic state is dependent upon the initial conditions and suggests that, even within a given topology, a host of meta-stable configurations exist, rather than a unique stable solution. We discuss a possible mechanism for the observed evolution, associated with the impact of higher frequency wavefronts emanating from parts of the knot which are more twisted than the expanding arcs.

\]

There are a wide variety of naturally occurring excitable media which possess spiral wave vortices. Examples include chemical concentrations in the Belousov-Zhabotinsky reaction  and electrical depolarization waves in cardiac tissue . The vortices in this last example are of particular significance since they are believed to play a vital role in ventricular fibrillation and hence sudden cardiac death . Both these systems, and many others, have a common mathematical description in terms of nonlinear partial differential equations of reaction-diffusion type. In the case of cardiac tissue the simplest continuous mathematical model is the FitzHugh-Nagumo equation, and it is the investigation of dynamical three-dimensional solutions of this equation which is the topic of this letter.

More than twenty years ago it was conjectured that stable, or at least persistent, three-dimensional solutions (termed ‘organizing centres’) might exist in excitable media in which two-dimensional vortices are embedded into three-dimensional space in such a way that they form knotted (or linked) vortex strings. The anatomy of these objects was clarified in terms of the topology of iso-concentration surfaces bordered by vortex strings . The hope was that the non-trivial topology of a configuration, perhaps aided by a short-range repulsive force between vortex cores, or by an effect of phase twist along the vortex string, might provide a barrier to its decay . However, others  argued against this optimism with the view that curvature, tension and reconnection processes would ultimately lead to the collapse and extinction of all knots. A framework was proposed  for thinking about vortex string dynamics in the limiting case of slight curvature and twist, but attempts to verify it were successful only in the strict limit of no twist . Ultimately, to address the fundamental issue of the existence of stable knots one must turn to numerical methods. About a decade ago a number of preliminary numerical investigations were performed  which suggested that certain knotted (and linked) configurations were stable, having a soliton-like behaviour in which the knot moved through the medium as a rotating rigid body with a constant shape. However, due to computational constraints, such simulations were limited to time scales which never exceeded about one hundred vortex rotation periods (often substantially less) and used very symmetric initial conditions.

In this letter we present the results of extensive numerical simulations of a duration well beyond a thousand times the vortex period, and using perturbed asymmetric initial conditions. We investigate several knots and links and conclude that all appear to be meta-stable in the sense that small perturbations produce dramatic changes in the evolution over time scales of the order of thousands of vortex rotation periods. The evolution is quite exotic, and very far from the simple curvature and tension-induced collapse suggested previously. In all cases the topology of the knot (or link) is preserved during the evolution, as we observe no reconnection events. Rather than a simple uniform contraction of the knot, which might be expected as a result of tension, we find that a particular arc of the knot both expands and contracts. After substantial variations in its total length the knot eventually approaches a steady state. However, this state does not appear to be unique and suggests that, even within a given topology, a host of meta-stable configurations exist. We discuss a possible mechanism for the observed evolution, associated with the impact of higher frequency wavefronts emanating from parts of the knot which are more twisted than the expanding arcs.

The FitzHugh-Nagumo equations are given by

$$\frac{u}{t}=\frac{(uu^3/3v)}{ϵ}+^2u,\frac{v}{t}=ϵ(u+\beta \gamma v)$$
(1)
where $`u(t,𝐱)`$ and $`v(t,𝐱)`$ are both real fields with $`u`$ the electric potential and $`v`$ the recovery variable associated with membrane channel conductivity. We take the constants appearing in equation (1) to have the values $`ϵ=0.3,\beta =0.7,\gamma =0.5.`$ This choice of constants is non-generic and is motivated by our aim of trying to find stable knots. This set of values has a number of special properties which might be conducive to knot stability, such as the lack of meander of a two-dimensional vortex and the stability of an untwisted vortex ring in three dimensions. See  for a description of the properties of a two-dimensional vortex as a function of the parameters $`ϵ,\beta ,\gamma .`$

In two space dimensions the FitzHugh-Nagumo equations with these parameter values has plane wave solutions which travel at a speed $`c=1.9`$ and rotating vortex solutions (often called spiral waves) with a period $`T_0=11.2.`$ The vortex solution has $`u`$ and $`v`$ wavefronts in the form of an involute spiral with a wavelength $`\lambda _0=cT_0=21.3.`$ Geometrically this means that all lines which are perpendicular to the level curves of the field $`u`$ are tangent to a small circle of diameter $`\lambda _0/\pi .`$ This circle represents the vortex core and is the region in space in which the gradients of the $`u`$ and $`v`$ fields differ substantially from being parallel. For later use it is convenient to define the quantity

$$\mathrm{\Phi }=|u\times v|$$
(2)
which is highly localized at the vortex core.

We solve equations (1) in three-space dimensions using an explicit finite difference scheme which is accurate to second order in the spatial derivatives and first order in the time derivative. Although this scheme appears very simplistic, it appears that the nature of these equations is such that more sophisticated or higher order algorithms do not lead to substantial gains in efficiency or accuracy, although this can be achieved if one is willing to modify the FitzHugh-Nagumo equations to a form designed specifically for the applicability of a more efficient numerical approach . For our simulations we use a grid containing $`201^3`$ points and a lattice spacing $`\delta x=0.5,`$ so that our spatial coordinates are confined to the range $`50x_i50.`$ The time step used is $`\delta t=0.02.`$ In the $`x_1`$ and $`x_2`$ directions we apply Neumann boundary conditions and in the $`x_3`$ direction the boundary conditions are periodic. The selection of the $`x_3`$ direction as periodic is because we shall orient our knots so that they initially translate as rigid bodies moving in the $`x_3`$ direction, and we do not wish to impede their motion.

We create initial conditions which form knotted vortex strings by making use of complex curves as described in  and is similar to the approach used in  for the study of knotted topological solitons. Recall that a knot may be written as the intersection of a complex curve $`𝒞`$ with the unit 3-sphere $`S^3.`$ Here $`S^3`$ should be thought of as compactified three-dimensional Euclidean space, with the explicit coordinates given by stereographic projection

$$Z_0=\frac{2(x_1+ix_2)}{1+r^2},Z_1=\frac{r^21+2ix_3}{1+r^2}$$
(3)
where $`r`$ is the Euclidean distance from the origin and $`Z_0`$ and $`Z_1`$ are two complex coordinates satisfying $`|Z_0|^2+|Z_1|^2=1`$ and hence parameterize $`S^3.`$ With this identification the knot is the one-dimensional locus in space of the complex curve $`𝒞`$ with coordinates $`Z_0`$ and $`Z_1.`$ As an example, to represent the $`(m:n)`$ torus knot we take $`𝒞=Z_1^mZ_0^n,`$ where for later convenience we have identified $`𝒞`$ with its zero set. If $`𝒞`$ has $`p`$ factors then it describes an object with $`p`$ components and hence this formalism can also be used to describe disconnected knots as well as links.

For a given knot (or link) we create initial conditions for the fields $`u`$ and $`v`$ from the associated curve $`𝒞`$ through the prescription

$$u=\mathrm{\Lambda }_1\text{Real}(𝒞)+u_{},v=\mathrm{\Lambda }_2\text{Imag}(𝒞)+v_{}.$$
(4)
Here $`\mathrm{\Lambda }_1`$ and $`\mathrm{\Lambda }_2`$ are two real constants (taken to be $`\mathrm{\Lambda }_1=2,\mathrm{\Lambda }_2=1`$) which are used to scale the initial conditions so that they cover the range of the excitation-recovery loop in $`(u,v)`$ space associated with the ODE part of equation (1). The constants $`u_{}=v_{}=0.4`$ are the values which can be roughly attributed to the vortex core.

The simplest non-trivial knot is the trefoil knot, given by the curve $`𝒞=Z_1^2Z_0^3.`$ We use this in the above prescription to obtain our initial conditions. In fig. 1 we plot the isosurface $`\mathrm{\Phi }=0.01`$, which indicates the core of the vortex string, at the times $`t=500,5000,10000,40000.`$ Note that at time $`t=0`$ this isosurface vanishes, since the initial conditions do not produce the vortex string itself but only seed the field configuration which will form into a vortex string after a time scale of the order of ten vortex periods. In fig. 1a the symmetric trefoil knot has clearly formed. In fact the knot forms at a much earlier time, but it is slightly larger and quickly shrinks to this size. The knot moves in the $`x_3`$ direction (towards the back left-hand-side of the box in the figure), with little change of shape, and at a speed of approximately $`c/80,`$ where $`c`$ denotes the wavefront speed given earlier. The knot also rotates around the $`x_3`$-axis with a period of around $`160T_0,`$ where $`T_0`$ is the vortex period given above. Fig. 1b reveals that the knot has drifted slightly away from the $`x_3`$-axis and one of the three lobes has expanded in comparison to the other two. Although the initial configuration has a cyclic $`C_3`$ symmetry, the cubic grid, and more importantly its boundary, breaks this symmetry and allows an asymmetric instability to develop. The larger lobe continues to expand, fig. 1c, and now the knot no longer simply translates in the $`x_3`$-direction, but rotates and follows a complicated path in space. Though still preserving the topology of a trefoil, the knot is now better viewed as a large expanding ring with a small knot tied in it. Eventually the expansion of the arc stops and a contraction begins. By fig. 1d the knot has regained its more symmetric form and has a similar length as before the expansion, but now it appears to be an asymptotic state. This can be seen by computing the length as a function of time, which is displayed in fig. 2.

To understand a possible mechanism responsible for the expansion of one arc of the knot and the subsequent contraction to a steady state we need to recall two facts. First, analytical and numerical work shows that a straight and uniformly twisted vortex line has a period which is slightly less than that of the two-dimensional vortex, or equivalently the untwisted vortex line . Here twist refers to the variation of the phase in the complex $`(u,v)`$-plane as one moves along the vortex string. Second, it is known that for a system of two vortices in which the vortices have different periods (for example, as arises in a model with spatially varying parameters) then the collision interface, which is the point at which the spiral wavefronts from the two vortices meet and annihilate, gradually moves towards the vortex with the larger period. In the absence of dispersion the collision interface moves at a speed $`\widehat{c}=c|T_1T_2|/(T_1+T_2)`$ where $`T_1`$ and $`T_2`$ are the periods of the two vortices. Eventually the collision interface reaches the core of the larger period vortex and it gets slapped away by the higher frequency wavefronts emanating from the shorter period vortex . Combining these two facts we see that a reasonable explanation for the expanding arc is that the more knotted part has a greater local twist rate than at least some part of the large expanding ring, so its period is slightly less and this results in its higher frequency wavefronts slapping away the ring and producing its expansion. This slapping mechanism is discussed in  in the context of stabilizing a knot against contraction. To check this hypothesis we have examined the collision interface by taking a slice through the configuration in which the expanding arc and most parts of the knot pass almost perpendicularly through the selected plane. This reveals that the wavefront produced by the tightly knotted cores impacts almost on top of the core of the expanding arc, in agreement with our hypothesis for slapping induced expansion. The details will be presented elsewhere .

The simplest example, in which the above issues regarding stability can be investigated, is to study two rings linked once. Consider the complex curve $`𝒞=Z_1^2Z_0^2\mu Z_0,`$ where $`\mu `$ is a real parameter. If $`\mu =0`$ then this curve is associated with two identical rings which each contain one full twist and are linked once. This configuration has a $`C_2`$ symmetry corresponding to a rotation by $`180^{}`$ around the $`x_3`$-axis. The link formed from this initial condition moves along the $`x_3`$-axis as a rotating rigid structure and shows no sign of instability even upto $`t=20000.`$ The reason this example differs from the trefoil knot in this respect is that a $`C_2`$ symmetry is clearly more compatible with the cubic lattice (and boundary) of the numerical grid than the $`C_3`$ symmetry of the trefoil knot. For this example we therefore require an explicit perturbation to test the stability of this link. This is achieved by setting $`\mu `$ to be non-zero in the above curve, which distorts one of the rings, making it larger than the other and hence breaking the $`C_2`$ symmetry. The results of a numerical evolution with $`\mu =0.5`$ are displayed in fig. 3, where we plot the vortex cores ($`\mathrm{\Phi }=0.01`$ isosurface) at times $`t=200,6200,10200,15200.`$ The larger of the two rings initially contracts but this is followed by an expansion which yields an asymptotic state in which the larger ring has a length similar to that in the perturbed initial condition. In fig. 4 we plot the lengths of the small ring (bottom curve) and the large ring (top curve) as a function of time. From this figure it can be seen that the length of the small ring remains almost constant and an asymptotic state has been reached which is certainly very different from the unperturbed solution ($`\mu =0`$) in which both rings have an equal length. By examination of the collision interface we again verify that the wavefronts from the small ring impact on the vortex core of the large ring. Moreover, an examination of the twist along each of the rings reveals that the small ring has a roughly constant positive twist along its length, but the large ring has a substantial variation in its twist rate, containing regions of negative twist even though the total twist along its length sums to one full turn in the positive direction. The fact that such a highly non-trivial distribution of twist occurs in an apparently asymptotic state is further evidence that a variety of meta-stable configurations exist in which the relative spatial distribution of the strings is in equilibrium under the action of several complicated forces in which the rate of twisting plays a vital role.

To summarize, we have found a novel dynamical behaviour of knotted vortex strings in the FitzHugh-Nagumo model with parameter values chosen to minimize any knot instabilities. It would be interesting to determine if our results are generic for the FitzHugh-Nagumo model with other parameter values and also for other excitable media. In fact, there is already some evidence for this in the initial expansion of a trefoil knot in a medium with equal diffusion of both reactants , but this example was regarded as an unexplained peculiarity at the time and simulations could not be performed for the length of time required to observe the full expansion and approach to an asymptotic state that we have described in this letter. It would certainly be worthwhile performing extensive numerical investigations, over very long time scales, on a variety of equations modelling different excitable media.

The construction of knotted vortex strings in laboratory experiments on excitable media would be of significant interest, though it is unlikely that the full evolution described in this letter could be studied in this setting since the typical lifetime of vortices in current experiments is limited to less than a hundred vortex periods.

Finally, the interaction and scattering of two initially well separated knots is also worthy of investigation.


format: text

text:
# Anomalous high-temperature Hall effect on the triangular lattice in NaxCoO2

## Abstract

The Hall coefficient $`R_H`$ of $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ ($`x=0.68`$) behaves anomalously at high temperatures ($`T`$). From 200 to 500 K, $`R_H`$ increases linearly with $`T`$ to 8 times the expected Drude value, with no sign of saturation. Together with the thermopower $`Q`$, the behavior of $`R_H`$ provides firm evidence for strong correlation. We discuss the effect of hopping on a triangular lattice and compare $`R_H`$ with a recent prediction by Kumar and Shastry.

The layered cobaltates are known to exhibit fairly high conductivities for perovskites. The discovery by Terasaki, Sasago, and Uchinokura Terasaki  of enhanced thermopower in the layered Na-doped cobaltate NaxCoO2 provided the first hint that these materials may harbor unusual electronic properties stemming from strong Coulomb interaction. Recently, interest in the cobaltates has escalated with the report of Takada *et al.* Takada  that superconductivity occurs below $``$5 K in the oxyhydrate $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2\mathrm{yH}_2\mathrm{O}`$. In addition, Wang *et al.* Wang  reported that the thermopower $`Q`$ in $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ ($`x=0.68`$) is completely suppressed by a 10-Tesla longitudinal magnetic field. The field suppression of $`Q`$, constituting direct evidence for a large spin-entropy term in the Seebeck coefficient, implies that a strong-correlation picture is necessary to describe the electronic transport properties. Superconductivity in $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2\mathrm{yH}_2\mathrm{O}`$ has been confirmed by several groups super .

Motivated by the observed field-suppression of $`Q`$, we have measured in detail the Hall effect in $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ ($`x`$ 0.68), and found that the dependence of the Hall coefficient $`R_H`$ on temperature $`T`$ is anomalous. Above 200 K, the Hall coefficient $`R_H`$ displays a steep, $`T`$-linear increase with no signs of saturation up to our highest $`T`$ (500 K). The linear increase continues through a previously unreported weak transition at $`T_D`$ = 430 K. In addition to the anomalous high-temperature behavior, the Hall angle and resistivity are also highly unusual at low-$`T`$.

Several groups have proposed phase diagrams for the layered cobaltates based on the resonating valence bond (RVB) theory applied to the triangular lattice Baskaran ; Shastry ; Lee . Independent of our Hall measurements, Kumar and Shastry Shastry  recently predicted that, on the triangular lattice, the high-frequency Hall coefficient $`R_H(\omega )`$ should be linear in $`T`$. We compare our experiment with this interesting prediction.

$`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ is comprised of layers of edge-sharing tilted octahedra. Each octahedron is made up of a Co ion surrounded by six O atoms at the vertices. Within each CoO2 layer, the Co ions occupy the sites of a triangular lattice. Doping is accomplished by Na ions which partially occupy the sites of a triangular lattice between the CoO2 layers. With the Na content fixed at $`x`$, the ratio of Co3+ to Co4+ is $`x:(1x)`$. As common in cobaltates, the splitting between the $`e_g`$ and $`t_{2g}`$ states in each octahedron is so large that the Co ions are always in their low-spin states. A further trigonal distortion splits the highest lying $`t_{2g}`$ state (called $`A_g`$) from the remaining 2 $`E_g`$ states. In the Co4+ ion, the $`A_g`$ state is occupied by one electron, while the $`E_g`$ states are completely full. According to band-structure calculation Singh , the Fermi level (for $`x=0.5`$) falls near the top of the band derived from the $`t_{2g}`$ states. At $`x`$ = 0 (all Co ions of formal valence Co4+), the band formed from the $`A_g`$ state is half-filled. Within the strong-interaction picture favored by the $`Q`$-suppression experiment Wang , this half-filled band is a Mott insulator. Increasing the Na content $`x`$ is equivalent to doping away from the half-filled Mott limit by adding electrons.

In our crystals, the Na content is determined Wang  by inductive coupled plasma analysis to be $`x`$ = 0.68, which implies that the ratio of $`\mathrm{Co}^{3+}:\mathrm{Co}^{4+}`$ is $`2:1`$. This is confirmed by a plot of the inverse in-plane susceptibility $`\chi ^1`$ which displays a Curie-Weiss behavior with a slope consistent with this ratio of Co4+ ions with $`S=\frac{1}{2}`$ (the Weiss temperature is $`\theta _W`$55 K) Wang .

As shown in Fig. 1, the in-plane resistivity varies linearly from 2 to 100 K with a slope $`d\rho /dT=2.45\mu \mathrm{\Omega }`$cm/K about 5 times larger than in optimally-doped $`\mathrm{YBa}_2\mathrm{Cu}_3\mathrm{O}_7`$. At high $`T`$, $`\rho `$ reveals a weakly hysteretic transition at $`T_D`$ = 430 K. Because the transition does not impact qualitatively the overall behavior of $`\rho `$, $`Q`$ and $`R_H`$, we associate it with an entropically driven order-disorder transition of the Na ions. The ‘metallic’ profile of $`\rho `$ suggestive of a large Fermi Surface (FS), juxtaposed with a large Curie-Weiss susceptibility Wang ; Ray  indicative of fluctuating local moments with antiferromagnetic (AF) coupling, presents an interesting dichotomy if we assume that the same electrons are responsible for transport and magnetization. In fact, the metallic profile of $`\rho `$ masks a set of anomalous transport properties which we describe next.

As noted Terasaki , the thermopower $`Q`$ is about 10 times larger than predicted by the Fermi-liquid expression $`Q=(\pi ^2/3)(k_B/e)(T/T_F)`$ with a Fermi temperature $`T_F10^4`$ K. A clue to this enhancement derives from the sensitivity of $`Q`$ to a longitudinal magnetic field $`𝐇`$. Wang *et al.* Wang  found that, below 20 K, $`Q`$ is increasingly suppressed by $`𝐇`$, until at 2 K, it is driven to zero by a field of 10 T. The suppression provides direct evidence that spin entropy contributes a dominant fraction of the enhanced $`Q`$. The curve of $`Q`$ vs. $`H`$ fits closely to a simple model of non-interacting spins with $`s=\frac{1}{2}`$ and a Lande $`g`$-factor of 2.2 $`\pm `$ 0.1 (Ref. Wang ). The plot in Fig. 1 shows that $`Q`$ increases with nominally constant slope $`dQ/dT`$ from 200 K to about 450 K, above which the slope increases slighly. The transition at $`T_D`$ is barely resolved as a weak peak.

In view of the unusual $`Q`$, it is natural to investigate other transport properties. The Hall coefficient was previously reported Terasaki  to be negative and only weakly $`T`$ dependent in the restricted range of $`T`$ between 4 and 130 K. When we extended the measurements above 200 K, we found the suprisingly strong $`T`$ dependence shown in Fig. 2. Near 200 K, $`R_H`$ changes sign and begins a steep increase that shows no sign of saturation up to our highest $`T`$ 500 K, apart from a weak anomaly at the transition at $`T_D`$. The slope of $`R_H`$ between 200 and 500 K is nearly independent of $`T`$ away from the anomaly at $`T_D`$. The ratio $`2:1`$ of the formal valence charge implies that the carrier density equals 8.65$`\times 10^{21}\mathrm{cm}^3`$, and a Drude value for the Hall coefficient of $`|R_H^0|=0.71\times 10^9\mathrm{m}^3/C`$. Here, we see that $`R_H`$ at 500 K is $``$8 times larger than $`R_H^0`$.

Because this Hall profile is so starkly different from the standard $`R_H`$ profile in conventional metals, we have adopted several precautions and checks. We used two methods that complement each others strengths. Method 1 is the standard practice of sweeping $`H`$ and measuring the Hall voltage $`V_H`$ at fixed $`T`$. At each $`T`$, $`H`$ is slowly ramped from -14 T to +14 T and then back to -14 T to yield 2 measurements of $`R_H`$. Method 2 is the technique of Sample *et al.* Sample  in which $`H`$ is fixed while $`T`$ is varied very slowly. Each measurement of $`R_H`$ is made by electronically toggling between the voltage and current leads to extract the antisymmetric component of the resistivity tensor. Method 1 allows checks for non-linearities in $`V_H`$ vs. $`H`$, but is error-prone if $`T`$ drifts significantly during the field ramp. Method 2 provides a high-density trace of $`R_H`$ vs. $`T`$ and is robust against contact noise, but is unreliable if $`R_H`$ varies with $`H`$. With method 2, $`R_H`$ is seen to have an anomaly at the transition $`T_D`$, but resumes its $`T`$-linear trend several degrees above $`T_D`$. Method 1 reveals that $`R_H`$ is strictly independent of $`H`$ up to 14 T above 20 K. Below 20 K, $`V_H`$ reveals a slight curvature vs. $`H`$. As shown in Fig. 2, the 2 methods agree well over the whole range of $`T`$.

In conventional metals, the empirical rule that $`R_H`$ must saturate to a constant above a ‘Hall scattering’ temperature $`\mathrm{\Theta }_H`$ is nearly universally obeyed in non-ferromagnetic metals. Here $`\mathrm{\Theta }_H`$ is the temperature above which the average momentum $`𝐪`$ of the available scattering excitation exceeds the largest Fermi Surface (FS) caliper. For phonon scattering, $`\mathrm{\Theta }_Hs\mathrm{\Theta }_D`$, where $`s`$ = 0.2-0.3 and $`\mathrm{\Theta }_D`$ is the Debye temperature. \[In $`\mathrm{Na}_{0.5}\mathrm{CoO}_2`$, $`\mathrm{\Theta }_D`$ = 380 K (Ref. Ando ), so we would have $`\mathrm{\Theta }_H`$ 100 K if phonon scattering were predominant.\] The basis for this rule is clearest in the two-dimensional (2D) FS, where the Hall conductivity $`\sigma _H`$ is proportional to the area $`𝒜_{\mathrm{}}`$ swept out by the mean-free-path $`\stackrel{}{\mathrm{}}(𝐤)`$ as the wavevector $`𝐤`$ goes around the FS Ong . For $`T>\mathrm{\Theta }_H`$, an electron can always encounter an excitation with large enough $`|𝐪|`$ to scatter it to an arbitrary point on the FS. The $`𝐤`$ dependence of $`\mathrm{}(𝐤)`$ then becomes independent of $`T`$ apart from an overall scale factor. Hence the *geometric shape* of $`𝒜_{\mathrm{}}`$ is $`T`$-independent, and its magnitude is rigorously proportional to $`\mathrm{}_{max}^2`$, where $`2\mathrm{}_{max}`$ is the maximum caliper of the swept area. As the conductivity $`\sigma \mathrm{}_{max}`$, this immediately implies that $`R_H`$ is independent of $`T`$. In the published literature of the Hall effect in metals that are non-ferromagnets Hurd , this empirical rule is rigorously obeyed (even in systems with both hole and electron FS pockets and in $`3d`$ metals with FS ‘monsters’). The only known exception is the cuprate family. The strong violation of this empirical rule provided an early, conspicuous clue Chien  that the charge transport in the cuprates is anomalous. In Fig. 2, the violation is just as striking as in the cuprates, but has the opposite trend.

A different perspective of the anomalous Hall effect is provided by the Hall angle $`\mathrm{tan}\theta _H=\sigma _H/\sigma `$, which is displayed in Fig. 3. In conventional metals, $`\mathrm{tan}\theta _H`$ is proportional to $`\mathrm{}_{max}`$, and must therefore decrease as $`T`$ increases. Here, we see just the opposite trend. $`\mathrm{tan}\theta _H`$ *increases* with $`T`$ right through the zero-crossing at 200 K as well as the disorder transition at $`T_D`$.

The sharp change in slope near 25 K in $`\mathrm{tan}\theta `$ vs. $`T`$ implies that the Hall response crosses over to a distinct regime at low $`T`$. As shown in the inset to Fig. 3, the magnitude of $`R_H`$ increases steeply below 25 K and saturates to a constant only below 3 K. However, $`\rho `$ is rigorously $`T`$-linear from 80 K down to $`2`$ K as noted above (see Fig. 4). These profiles of $`\rho `$ and $`R_H`$ are reminiscent of the corresponding curves in the cuprates (except there the curves extend to 500 K or higher and $`R_H`$ is positive). Following the analysis of Chien *et al.* Chien  in the cuprates, we find that $`\mathrm{cot}\theta _H`$ is well-fitted by the expression $`a+bT^n`$ with $`n=1.5`$ (dashed line in Fig. 4). This implies that the strong $`T`$ dependence of $`|R_H|`$ below 25 K results from distinct power-law exponents for the Hall-angle lifetime ($`n=1.5`$) and the transport lifetime obtained from $`\rho `$ $`(n^{}=1)`$. Interestingly, the temperature scale over which this holds is about 10-20 times smaller than in the cuprates. This analysis suggests that, below 25 K, the carriers in $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ enter into a strong-correlation regime that involves the same ‘strange-metal’ physics as observed in the cuprates. The reduction of the Hall-angle exponent from 2 to 1.5 also matches well the observed reduction from optimally doped to overdoped cuprates.

We return to the high-temperature behavior of $`R_H`$. The Hall effect problem in cuprates has motivated many theoretical investigations of $`R_H`$ in the $`t`$-$`J`$ model. Recently, Kumar and Shastry (KS) applied to $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ a series expansion technique SSS  to calculate the high-frequency limit of the Hall coefficient $`R_H^{}`$ (the small parameter is $`\beta t`$ where $`t`$ is the hopping amplitude and $`\beta =1/k_BT`$). For a triangular lattice, KS obtain Shastry ; SSS

$$R_H^{}=\frac{a^2d}{8|e|}\frac{k_BT}{t}\frac{(1+\delta )}{\delta (1\delta )},$$
(1)
where $`a`$ is the Co-Co bond length, $`d`$ the interlayer spacing, and $`\delta `$ the doping away from half-filling.

In Fig. 2, $`R_H`$ is $`T`$-linear between 200 and 500 K. If we fit this linear portion to Eq. 1, we find that $`t`$ = 25 K, consistent with $`\beta t<1`$. However, it is not obvious *a priori* that the prediction for $`R_H^{}`$ can be compared meaningfully with the dc data in Fig. 2 (moreover, the theory does not predict the weaker dependence on $`T`$ below 100 K). It is clearly important to verify theoretically whether the domain of validity of Eq. 1 extends to the dc regime in the $`t`$-$`J`$ model.

The singular behavior of $`R_H`$ is intimately related to charge transport on the triangular lattice. It is helpful to recall electrons hopping on a lattice (Fig. 2, inset). The Hall current arises from the Peierls phase factor $`\mathrm{e}^{i\alpha }`$ acquired around the minimal closed loop, where $`\alpha =e\varphi /\mathrm{}`$ is proportional to the flux $`\varphi `$ piercing the loop. For the 3-site problem (site indices 1,2,3), Holstein Holstein ; Friedman  computed the transition rate for hopping from site 1 to 3 in a field, and found that the direct, first-order hopping amplitude $`\widehat{t}_{31}`$ interferes with the second-order amplitude $`\widehat{t}_{32}\widehat{t}_{21}`$ (with $`\widehat{t}_{ij}`$ the complex hopping amplitude between sites $`i`$, $`j`$ in a field), to yield a giant enhancement of the Hall mobility (Fig. 2, inset) Holstein , so that $`\sigma _H\widehat{t}_{13}\widehat{t}_{32}\widehat{t}_{21}t^3\mathrm{e}^{\mathrm{i}\alpha }`$. However, this non-interacting, hopping model obviously does not describe the physics here.

To get a $`T`$-linear $`R_H`$, we need to consider strong correlation. In the $`t`$-$`J`$ model calculation of KS, transport around a loop on the triangular lattice gives likewise $`\sigma _H(\beta t)^3\mathrm{sin}\alpha `$ Shastry . Since $`\sigma (\beta t)^2`$, we have $`R_H(\beta t)^1`$, which implies a $`T`$-linear increase. In the strong-correlation problem, the singular behavior in Eq. 1 is a consequence of transport on the triangular lattice.

The accumulated evidence increasingly support the strong-interaction picture in the layered cobaltate. Initially, observations of an enhanced $`Q`$ Terasaki  and a Curie-Weiss susceptibility Ray  in the presence of a ‘metallic’ $`\rho `$ provided early hints of strong correlation. Later, direct evidence for dominant spin entropy contribution to the thermopower came from the effect of $`H`$ on $`Q`$ Wang . The steep increase in $`R_H`$ at high temperatures, as well as the cuprate-like behavior of $`R_H`$ and $`\rho `$ below 25 K, reported here provide further evidence for strange-metal behavior. Both the field suppression of $`Q`$ and the Hall profile in Fig. 2 provide striking examples of how strong interaction can change *qualitatively* electronic transport behavior from conventional behavior even at temperatures up to 500 K. Hopefully, a systematic search may uncover more examples of $`T`$-linear $`R_H`$ in layered perovskites with triangular lattices.

We thank B. S. Shastry for enlightening discussions, and acknowledge support from a MRSEC grant (DMR 0213706) from the National Science Foundation.


text:
# Universal threshold enhancement

## Abstract

By assuming certain analytic properties of the propagator, it is shown that universal features of the spectral function including threshold enhancement arise if a pole describing a particle at high temperature approaches in the complex energy plane the threshold position of its two-body decay with the variation of $`T`$. The case is considered, when one can disregard any other decay processes. The quality of the proposed description is demonstrated by comparing it with the detailed large $`N`$ solution of the linear $`\sigma `$ model around the pole-threshold coincidence.

sigma meson, spectral function, finite temperature, linear sigma model

The near-threshold two-body decay of a resonance is a well-studied problem in scattering theory walecka . In a medium with tunable parameters (e.g. temperature, density etc.) new features may arise yoko . It has been found hatsuda ; kuni  that a hallmark of partial chiral symmetry restoration in a medium is the enhancement of the scalar-isoscalar spectral function near the two-pion threshold. In order to elucidate further the origins and the nature of this phenomenon, in this note we formulate a general approach to it. On the basis of this investigation some features of threshold enhancement appear to be universal, showing up whenever the mass of a stable particle can be driven by tuning some parameters of its environment to pass the threshold of a two-body decay into stable final particles. Following our general discussion we shall show evidence for this situation in the linear $`\sigma `$ model by studying the propagator of the $`\sigma `$ particle, neglecting the finite in-medium lifetime of the pions.

Let us start our general discussion by assuming the existence of a branch point in the complex $`\sqrt{p^2}`$ rest frame energy plane on the real energy axis for the propagator $`G_\sigma (p^2,|𝐩|),(p^2=p_0^2𝐩^2)`$, which corresponds to the $`\sigma \alpha _1+\alpha _2`$ hypothetical decay. (In this part of our discussion $`\sigma `$ refers to an arbitrary particle. The extra dependence on the 3-momentum is a medium effect.) We can restrict the discussion to the halfplane $`\mathrm{Re}\sqrt{p^2}>0`$, due to the symmetric behavior of the propagator. At the branch point a square root singularity, proportional to $`\delta \sqrt{1p^2/M^2}`$, shows up in the self-energy contribution to the propagator, where $`M`$ denotes the temperature dependent position of the branch point. The analytic (Riemann-sheet) structure of the $`\sqrt{p^2}`$-plane is specified by requiring that $`G_\sigma `$ is analytically continued through the part of the real axis where $`p^2>M^2`$, i.e. where the complex threshold factor $`\delta `$ is purely imaginary. Furthermore the plane we are working on is bounded by a cut along the section $`p^2<M^2`$ of the real axis. Along the cut $`\delta `$ is real positive (negative) on the upper (lower) halfplane. Conventionally, one associates the upper halfplane with the first, the lower one with the second Riemann-sheet. We shall call the upper edge of the real axis reached from the first Riemann sheet as physical, and the lower edge of the cut as unphysical.

The goal of our analysis is to investigate the threshold behavior of the spectral function in the neighborhood of the branch point in a small interval around some temperature $`T^{}`$, to be specified below. We first restrict the discussion to the case $`|𝐩|=0`$, which makes our analysis easier to follow. There are now two small quantities characterizing this situation: $`ϵ(TT^{})/T^{}`$ and $`\delta `$. \[It is important to note, that it is $`M(T)`$ and not $`M(T^{})`$, which appears in the definition of $`\delta `$.\] Our basic assumption is that the inverse propagator $`G_\sigma ^1(p_0)`$ can be expanded in terms of these two independent variables. Then to quadratic order one obtains

$`G_\sigma ^1(p_0)`$ $``$ $`a_1ϵ+\stackrel{~}{a}_2ϵ^2+b_1\delta +b_2\delta ^2+c\delta ϵ`$ (1)
$``$ $`{\displaystyle \frac{1}{1+cϵ/b_1}}[a_1ϵ+a_2ϵ^2+b_1\delta +b_2\delta ^2],`$
where $`a_2=\stackrel{~}{a}_2+a_1c/b_1`$ was defined, after factoring out a temperature dependent finite constant.

We investigate the behavior of the spectral function, which can be derived from (1) by studying the pole structure starting at high temperature. In the concrete case of the $`\sigma `$-model the tendency to restore the chiral symmetry with increasing temperature reduces the $`\sigma `$-mass below $`2m_\pi (T).`$ For such temperature $`\sigma `$ is necessarily stable if the pions are assumed to be stable.

When decreasing the temperature the location of the $`\sigma `$-pole moves along the physical real axis to the two-body decay threshold from below. Let us assume that this pole reaches the branch point at $`T=T^{}`$. This is the temperature $`T^{}`$ around which the behavior of the spectral function will be investigated.

The stability of the $`O(N)`$ multiplet enforced by the chiral symmetry is now postulated, without reference to any particular symmetry, for the general discussion. We assume that for $`ϵ>0`$ in addition to the approximate stability of the particles $`\alpha _1`$ and $`\alpha _2`$, $`\sigma `$ is stable. This requires a single pole, $`\delta _+`$, for $`G_\sigma `$ on the physical real axis below the threshold $`(1>\delta >0)`$, corresponding to a stable state. The requirement, that no other positive root should be present is fulfilled only if $`a_1<0,b_1>0,b_2>0`$. This choice implies the existence of a root $`\stackrel{~}{\delta }<0`$ on the unphysical real axis for both signs of $`ϵ.`$ This pole stays away from the threshold even when $`ϵ0`$, see Fig. 1.

When the temperature is decreased to the region $`ϵ<0`$, $`G_\sigma ^1`$ has no roots on the physical real axis. The pole, which did correspond to the stable particle moves continuously over the threshold to the unphysical real axis. In this temperature range this solution is denoted by $`\delta _{}`$. This means that for $`ϵ<0`$ one has two poles $`(\delta _{},\stackrel{~}{\delta })`$ on the unphysical real axis for sufficiently small $`|ϵ|`$.

In this approximation the spectral function $`\rho (ϵ,\delta )=\mathrm{Im}G_\sigma (p_0)/\pi `$ has the following parameterization above the threshold $`p_0^2>M^2`$:

$$\rho (p_0,ϵ)=\frac{\left[1+\frac{c}{b_1}ϵ\right](b_2\pi )^1\sqrt{\frac{p_0^2}{M^2}1}(\stackrel{~}{\delta }+\delta _\pm )}{(\frac{p_0^2}{M^2}1\delta _\pm \stackrel{~}{\delta })^2+(\frac{p_0^2}{M^2}1)(\stackrel{~}{\delta }+\delta _\pm )^2}$$
(2)
(the expression depends on $`ϵ`$ also indirectly, through the $`\stackrel{~}{\delta }\delta _\pm `$ term). The pole terms associated with $`\delta _+`$ and $`\delta _{}`$ in the respective temperature ranges $`ϵ>0`$ and $`ϵ<0`$ can be uniquely separated in $`G_\sigma `$. The approximate spectral function based exclusively on this singular “pole” contribution, above the threshold $`M`$ reads:

$$\rho _{\mathrm{sing}}=\frac{(1+cϵ/b_1)}{b_2(\delta _\pm \stackrel{~}{\delta })\pi }\times \frac{|\delta |}{|\delta |^2+\delta _\pm ^2}.$$
(3)
The maximal threshold enhancement clearly occurs for $`ϵ=0`$, when $`\rho _{\mathrm{sing}}1/|\delta |`$, which has been found also in Ref. hatsu . Since $`\delta _+(ϵ)`$ and $`\delta _{}(ϵ)`$ in the respective temperature ranges are equal to linear order in $`ϵ`$, the mirror variation to $`𝒪(ϵ)`$ of $`\rho _{\mathrm{sing}}`$ around $`T^{}`$ is deformed only by the temperature dependent first factor on the right hand side of Eq. (3).

The threshold enhancement sets in in a universal way, reflecting the motion of $`\delta _+(\delta _{})`$ along the physical (unphysical) real axis towards the position of the threshold. Concerning universality based on the expansion given in Eq. (1) three further remarks are in order: i) the linear dependence of $`\delta _\pm `$ on $`ϵ`$ is universal, which determines the leading term of the enhancement according to the second factor in Eq. (3), ii) when looking at the first correction to this leading behavior the one-pole description is not sufficient, $`\stackrel{~}{\delta }`$ plays an important role (if the expansion is pushed to higher orders, further roots would appear), iii) it is obvious that the expansion given in Eq. (1) is independent of the low temperature quasiparticle spectra.

Another characteristic temperature value may be defined, corresponding to the temperature where $`\delta _{}`$ and $`\stackrel{~}{\delta }`$ become degenerate on the unphysical real axis, when the temperature is gradually decreased further. It is determined by the smaller negative root of the equation $`b_1^24b_2(a_1ϵ_0+a_2ϵ_0^2)=0.`$ For smaller temperature values $`(ϵ<ϵ_0)`$ the two poles become complex. One of the poles moves into the second Riemann sheet and can be interpreted as a resonance of finite lifetime if the temperature becomes sufficiently low. The other pole is its “mirror” located on the continuation of the second Riemann sheet through the unphysical real axis into the positive halfplane, and does not have any physical meaning. This picture fully agrees with the qualitative pole trajectory found in nagyN2 . It could happen that the validity of the quadratic expansion of $`G_\sigma ^1`$ breaks down already before the degeneracy of the two poles on the unphysical real axis occurs and $`\delta _{}`$ may stay on this axis until the temperature is decreased to zero. This alternative seems to be realized in the calculation presented in hidaka . It would be very interesting to see the presence of such a pole in the low temperature representation of the pion-pion scattering amplitude based on chiral perturbation theory dobado .

The extension of the investigation to the case of nonzero spatial momentum $`|𝐩|/M`$, regarded as a new small parameter can be readily done, by assuming the analyticity of $`G_\sigma `$ also on this quantity. Now the threshold parameter $`\delta `$ takes its original Lorentz invariant form: $`\delta =\sqrt{1p^2/M^2}`$, where $`p^2`$ is the squared four-momentum. At finite temperature, in the propagator also a separate dependence on the spatial momentum appears. Therefore the expansion Eq. (1) will contain also terms of the form $`(𝐩^2/M^2)^lϵ^k\delta ^m`$ (note, that $`ϵ`$ remains as it was introduced for $`|𝐩|=0`$!). To lowest order (single pole approximation) one keeps terms linear in $`ϵ`$ and $`\delta `$ and the coefficient $`d`$ of a new term $`d\times |𝐩|^2/M^2`$ has to be computed.

The coefficient $`d`$ determines the shift in the value of the temperature for which the maximal threshold enhancement of $`\rho _{sing}`$ occurs at small finite $`|𝐩|`$: $`T^{}(|𝐩|)=T^{}[1+\frac{d}{4|a_1|}\frac{𝐩^2}{m_\pi ^2}]`$. Since the term proportional to $`𝐩^2`$ is real even after the continuation above the threshold, $`|𝐩|`$ and $`ϵ`$ play similar roles in the spectral function. For fixed $`T`$, it is $`d`$, which controls the location and height of the maximum of $`\rho _{sing}`$ when $`|𝐩|`$ is increasing. For instance, one finds $`\delta _{max}|a_1ϵ+d|𝐩|^2/M^2|/b_1`$.

Below, we illustrate the validity of the simplified representation discussed above for the scalar-isoscalar spectral function of the linear $`\sigma `$ model, obtained to leading order in an expansion with respect to the inverse of the number of the Goldstone bosons (large $`N`$ approximation). For this approximate solution one checks explicitly the existence of the power series of $`G_\sigma ^1`$ in $`\delta ,`$ $`|𝐩|`$ and $`ϵ`$. Furthermore, the ranges of validity in $`|ϵ|`$ of the single and the two pole descriptions of the spectral function can be analyzed by comparing Eq. (2) and Eq. (3) with the complete expression of the spectral function in the model.

The first goal in the second part is to sketch the way the $`\sigma `$ propagator given by Eq. (19) of Ref. nagyN2  can be cast for $`p_0<M2m_\pi (T)`$ and $`|𝐩|=0`$ in the form presented in Eq. (1) of the present paper. Employing an obvious notation, the renormalized couplings $`(m_R^2,\lambda _R,h)`$ of this model were fixed at $`T=0`$ at values which ensure $`m_\pi (0)=140\mathrm{MeV},f_\pi =\sqrt{N}\mathrm{\Phi }(T=0)=93\mathrm{MeV},`$ $`m_\sigma =3.95f_\pi ,`$ $`m_\sigma /\mathrm{\Gamma }_\sigma 1.4`$. The temperature dependence of the vacuum expectation value of the $`O(N)`$ field, $`\mathrm{\Phi }(T)`$, was determined by the equation of state (its renormalized form appears in Eq. (11) of Ref. nagyN2 ):

$$m^2+\frac{\lambda }{6}\left[\mathrm{\Phi }^2(T)+_k\left[n(\omega _k)+\frac{1}{2}\right]\frac{1}{\omega _k}\right]\frac{h}{\mathrm{\Phi }(T)}=0,$$
(4)
where $`n`$ is the Bose-Einstein factor and $`\omega _k^2=𝐤^2+m_\pi ^2`$. The expansion of $`m_\pi (T)`$ to quadratic order in $`ϵ`$ is provided through the combination of the relation $`m_\pi ^2(T)=h/\mathrm{\Phi }(T)`$ with the expansion of $`\mathrm{\Phi }(T)`$ around $`\mathrm{\Phi }(T^{})`$. $`\mathrm{\Phi }(T^{})`$ is found from Eq. (31) of Ref. nagyN2 .

The selfenergy of $`\sigma `$ is determined to leading large $`N`$ order by a sum of the chain of pion bubbles, characterized by the function $`b^>(p_0)`$ as shown by Eq. (5) of Ref. nagyN2  (the index $`>`$ refers to the first Riemann sheet). The expression of the bubble above the threshold given in Ref. nagyN2  is written, after continuing to the physical real axis below the threshold and some convenient integral transformations, in a form where the expansion in powers of $`\delta `$ and $`ϵ`$ actually starts from:

$`4\pi ^2b^>(p_0)={\displaystyle \frac{1}{4}}\mathrm{ln}{\displaystyle \frac{m_\pi ^2(T)}{M_0^2}}+{\displaystyle \frac{\delta \mathrm{arccos}(\delta )}{2\sqrt{1\delta ^2}}}`$
$`{\displaystyle \underset{\xi }{\overset{\mathrm{}}{}}}𝑑x{\displaystyle \frac{(x^2\xi ^2)^{\frac{1}{2}}}{e^x1}}+\delta {\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑x{\displaystyle \frac{[(1+x^2)\sqrt{1+\delta ^2x^2}]^1}{e^{\xi \sqrt{1+\delta ^2x^2}}1}},`$ (5)
where $`m_\pi (T)`$ and $`\xi =m_\pi (T)/T`$ carry the explicit $`T`$ dependence ($`M_0`$ is the normalization scale of the theory defined at $`T=0`$). One performs first the expansion in powers of $`\delta `$ to quadratic order. For the second integral in Eq. (5) one finds to linear order $`(\frac{\pi }{2}\delta )n(m_\pi (T))`$. Next, the remaining explicit pion mass dependence (also in $`\xi `$!) is expanded around $`T^{}`$ with the procedure sketched above. With these hints one easily establishes the approximate expression Eq. (1), obtaining numerically the following values for the coefficients: $`a_1=18.32`$, $`a_2=23.14`$, $`b_1=8.59`$, $`b_2=14.48`$, $`c=9.67`$.

In Fig. 2 the spectral function and its universal approximations are shown for three different temperatures. The curves labeled “exact” are calculated from the complete expression derived in Ref. nagyN2 . They are compared with the two approximate forms described above. All temperatures were chosen below $`T^{}`$ (above it the situation is qualitatively similar). For the highest value of $`ϵ`$ the pole $`\delta _{}`$ lies almost on the threshold, for its smallest value $`(ϵ=ϵ_00.077)`$ one has $`\delta _{}=\stackrel{~}{\delta }`$. When calculated from the exact expression one obtains: $`ϵ_0=0.13.`$ According to these curves the region where the one-pole description gives a good approximation is rather narrow in $`ϵ`$ and in particular in $`\delta .`$ The correction due to the second pole contribution is quite important for increasing frequency values. For $`ϵ=ϵ_0`$ the line corresponding to $`\rho _{sing}`$ was omitted, since its expression becomes meaningless. It is, however, remarkable how good approximation can be constructed using Eq. (1). The range of validity of the quadratic expansion of $`G_\sigma ^1`$ in $`ϵ`$ includes also the degeneracy temperature of the poles, persisting even for somewhat lower temperatures when both poles are complex and lie equidistantly from the threshold. The specific form of the spectral function, $`\rho |\delta |/(|\delta |^2+const)^2`$, characteristic at the degeneracy temperature is expected to occur even in the exact solution of the $`N=4`$ linear $`\sigma `$ model with physically relevant couplings.

For the case of the $`\sigma \pi `$ system, the coefficient $`d`$, which determines the influence of the explicit $`|𝐩|`$-dependence, comes uniquely from the temperature dependent part of the pion bubble, $`b_T^>b_T^>(p_0,|𝐩|)`$. Since the $`T=0`$ part of the expression of the bubble is Lorentz invariant, it remains unchanged when expressed with the complete expression of $`\delta `$. The expression of $`b_T^>`$ for $`𝐩0`$, see e.g. Refs. hidaka ; toni  can be cast in the form

$$b_T^>=\underset{1}{\overset{\mathrm{}}{}}\frac{dt}{8\pi ^2}\frac{n(tm_\pi )}{|𝐩|}\mathrm{ln}\left|\frac{\left[\frac{p^2}{2}+|𝐩|\sqrt{t^21}\right]^2p_0^2t^2}{\left[\frac{p^2}{2}|𝐩|\sqrt{t^21}\right]^2p_0^2t^2}\right|.$$
(6)
For the threshold behavior of the spectral function to the lowest (linear) order in $`(ϵ,\delta ,𝐩^2)`$ , one can use $`\delta =0`$ in the integral, that is for the quantities scaled by $`m_\pi (T)`$ one can put: $`p^2=4,`$ $`p_0^2=4+𝐩^2.`$ With a suitable transformation of the variable of integration we can ensure the smoothness of the integrand, obtaining:

$$b_T^>=\underset{0}{\overset{\mathrm{}}{}}\frac{dt}{4\pi ^2\xi }\mathrm{ln}\frac{1\mathrm{exp}\left[\xi \sqrt{1+\left(\frac{|𝐩|}{2}\mathrm{tanh}\frac{|𝐩|}{2}t\right)^2}\right]}{1\mathrm{exp}\left[\xi \sqrt{1+\left(\frac{|𝐩|}{2}\mathrm{coth}\frac{|𝐩|}{2}t\right)^2}\right]}.$$
(7)
Performing a partial integration one obtains the first two terms of the asymptotic expansion of the integral valid for small $`|𝐩|`$, relevant for our analysis. Next to the unchanged $`|𝐩|`$-independent term with the complete $`\delta `$ one can calculate the coefficient of $`𝐩^2`$ by performing the $`t`$-integral with the value of $`\xi `$ determined at $`T^{}`$. For the actual couplings $`d=1.65`$ is found.

The assumed expansion of the inverse propagator around the threshold in powers of $`ϵ`$ and $`\delta `$ resembles the starting point of the Landau theory of second order phase transitions. Like in the Landau theory, our purpose has not been to prove the analytic property of the relevant quantity, but to explore the consequences of its postulation. We conjecture the universality of the threshold behavior in any channel where a particle pole describing at high temperature a stable particle might pass through the threshold of its lowest two-body decay. Our discussion implies that for $`TT^{}`$ the pole dominating analytically the spectral function cannot correspond to any particle. Whether this remains valid down to $`T=0`$ is model dependent.

A possible generalization of Eq. (2), in which only the quadratic dependence of the inverse propagator on $`\delta `$ is assumed, can be easily written down:

$$\rho (ϵ,𝐩,p^2)=\frac{R_1(ϵ,𝐩)\sqrt{\frac{p^2}{M^2}1}}{[\frac{p^2}{M^2}1+R_2(ϵ,𝐩)]^2+R_3(ϵ,𝐩)[\frac{p^2}{M^2}1]}.$$
(8)
The functions $`R_i`$ are real, at $`T=0`$ they assume constant values. $`ϵ`$ is measured from the value $`T^{}`$ determined by the condition $`R_2(ϵ=0,0)=0`$. A sequence of reduced temperatures $`ϵ(𝐩)(T^{}(𝐩)T^{})/T^{}`$ can be defined as the solution of the equation $`R_2(ϵ(𝐩),𝐩)=0`$. Obviously $`T^{}(0)=T^{}`$. Asymptotically one expects the “scaling” behavior $`R_2[ϵϵ(𝐩)]^\alpha ,\alpha >0.`$ The treatment based on Eq. (1) leads to $`\alpha =1`$, but in principle a nontrivial non-integer scaling exponent can be envisaged. Note, that when there is no explicit symmetry breaking the branch point tends to the origin if $`|𝐩|`$ goes to zero and $`T^{}=T_c`$ independently of $`|𝐩|`$ nagyN1 .

The next step will be to discuss the effect of the finite pion lifetime. It shifts the branch-point into the lower halfplane. Around it there is still an enhancement (similar to the undamped case) but in the spectral function varying along the real axis the effect is smeared out partially, see ohtani . Either a complete two-loop (improved) perturbative analysis or the evaluation of the next-to-leading large N correction to $`G_\sigma `$ is needed to reliably assess the importance of the finite pion width.

This research has been supported by the research contract OTKA-T037689 of the Hungarian Research Fund.


text:
# Do Bulges of Early- and Late-type Spirals Have Different Morphology?

## 1 Introduction

Understanding bulges of galaxies is of pivotal importance for improving our overall theories for galaxy formation and evolution. Yet today, our general understanding of bulges is limited, and little is known about the formation scenario they have gone through. Up to the 1980’s, ideas about bulge formation were very much influenced by Baade’s population concept (e.g. Sandage 1986). Since in our Galaxy the stars in both the bulge and the stellar halo were found to be old, it was thought that the bulges and the stellar galaxy halo were a single entity, with a common $`r^{\frac{1}{4}}`$ surface brightness (SB) law (de Vaucouleurs 1948, 1959). They were thought to be formed early on in the lifetime of the galaxy, during a rapid collapse phase (Eggen, Lynden-Bell & Sandage 1962) before the formation of the disk. Their shapes were thought to be oblate spheroids (Kent, Dame & Fazio 1991). In recent years, however, a wealth of new data has altered our view of galactic bulges. Misalignments of the major axis with respect to the disk major axis indicate that bulges are probably not oblate (Bertola, Vietri & Zeilinger 1988). Not all stars in the bulge region are old (e.g. Böker et al. 2002). And the surface brightness profiles of many bulges do not follow the $`r^{\frac{1}{4}}`$ law - there are claims nowadays that there are no bulges for which the surface brightness profiles obey this law (Balcells et al. 2003). Bulges do seem to be rotationally flattened (Kormendy & Illingworth 1982, Davies & Illingworth 1983), with some bulges rotating even faster, almost as fast as disks (Kormendy 1993).

Before delving into the subject, it is good to first clearly define what we mean by bulge. In this paper, by bulge we mean the photometric inner component of the galaxy which appears to be superimposed onto an exponential disk. This definition is extremely simple, but is at the same time very robust. The definition is purely photometric, and is parametrical, but non-parametric studies by e.g. Kent (1986) show that the bulge obtained in this way is generally rounder than the disk, an indication that this definition is probably measuring something physical. Note that the definition does not depend on kinematics. When one measures the photometry of bulges defined in this way one finds that surface brightness profiles for galaxies of type earlier than Sc rise more steeply than those of late-types. Surface brightness profiles generally are well described by the Sérsic (1968) $`r^{\frac{1}{n}}`$ profile, for which the value of the shape parameter $`n`$ varies from values around 4 for bulges of early-type spirals, to 1 for bulges of late-type spirals (Andredakis, Peletier & Balcells 1995; Phillips et al. 1996; Moriondo, Giovanardi & Hunt 1998). Bulges also host a variety of central components including dust lanes, star forming rings and spiral structure reaching all the way to the centre (Zaritsky, Rix & Rieke 1993; Peletier et al. 1999; Carollo et al. 2002). Many bulges also contain central resolved sources, identified as star clusters (Carollo et al. 2002, Böker et al. 2001), which are fainter and less abundant in early-type spiral galaxies. Studying the nuclear properties of bulges is important, since they generally scale with the global galaxy properties (Faber et al. 1997). By studying the nuclear properties of bulges of early and late-type spirals, one might be able to provide the relative importance of dynamical effects in forming and/or maintaining the nuclear structure. It also provides a crucial test for formation scenarios of spheroidal stellar systems along the entire luminosity sequence.

The main theories about bulge formation include scenarios such as: primordial collapse where bulge and disk form sequentially (Eggen, Lynden-Bell & Sandage 1962); hierarchical galaxy formation (Kauffmann & White 1993); galaxy mergers or infall of satellite galaxies (Zinn 1985, Aguerri, Balcells & Peletier 2001); and secular evolution (Kormendy 1979, Norman, Sellwood & Hasan 1996, Hasan, Pfenniger & Norman 1998), where central bars, formed by disk instabilities, thicken gradually to form bulges. The theories are still of crude, and theoretical predictions up to now have not been good enough for observers to rigorously test these scenarios. Effects of flattening by rotation, or thickening of the bars have not directly been verified in simulations for large numbers of objects, mainly due to lack of fast computers. There are a number of crude predictions that one would naively draw from these formation models. If bulges were formed by secular evolution of bars, they would in general be flatter than if they were made in a primordial collapse. One might object against such simplified predictions, but since there is no alternative until robust results from simulations are available, we will assume that these predictions are correct.

In this paper, we investigate the ellipticity distribution of a large sample of early- and late-type spirals, to ultimately be able to infer information about the formation scenario of bulges. Given that in the infrared, the effects of dust extinction is considerably less than in optical, HST/NICMOS H-band data are ideal for this task. We discuss a sub-sample of late-type spirals consisting of HST-archival data, together with a control sample of early-type spirals consisting of HST-archival images of both Seyferts and non-Seyferts from Laine at al. (2002) (hereafter L02). Since the NICMOS field of view is not large enough to provide Bulge-Disk decompositions, we use WFPC2 optical images for determination of the bulge radius, the radius where the light of the disk starts dominating the galaxy light. Our samples are large enough to allow us to draw statistical conclusions using the KS-test. The outline of the paper is as follows: In Sect. 2 we introduce our data sets, and describe the reduction procedure. In Sect. 3 the analysis procedure is outlined, and Sect. 4 includes an extensive discussion of our results and corresponding implications.

## 2 The Data

We study a sample of archival HST/NICMOS (NIC2) F160W ($`\lambda _e1.6\mu `$ m, roughly corresponding to H-band) data. The instrument and its performance are described in Thompson et al. (1998) and the NICMOS instrument handbook (MacKenty et al. 1997). For collecting the data the well known large observational data sets including samples by Mulchaey ID=7330, Peletier ID=7450, Pogge ID=7867, and Stiavelli ID=7331, were searched and investigated. These four samples are the only large samples of bulges in the HST archive, and include galaxies covering a wide range of morphological types. Of these samples, all galaxies of type later than or equal to $`T_{RC3}=4`$ (corresponding to Sbc) are studied, and the earlier types are selected from the L02 sample. The \*_cal.FITS images were used as a starting point, after which additional data reduction steps were performed. Processed by the CALNICA pipeline (NICMOS Data Handbook), these images have been science calibrated, dark subtracted, flat fielded, and cosmic ray corrected in an automatic way. Visual inspection of the images showed a number of anomalies, and it was evident that some further reduction needed to be implemented.

Sky subtraction was done by subtracting the mean value of 10 randomly chosen outer regions of each image. Since in this study we are mainly interested in the central morphology, it is not important that the galaxies are much larger than the NIC2 frame. Masking bad pixels was done in two steps. First, only the most evident regions were masked by fitting planes interpolated from the boundary values (using our own software). This was followed by fitting ellipses to the isophotes (using Stsdas/Ellipse), and investigating the residual images, which showed the presence of additional undesired regions, such as very bright star forming knots, foreground stars, or very strong dust regions. After further masking of these regions, all the obvious artifacts were removed, and the images were ready for final ellipse fitting. The Ellipse task in the Isophote package in Stsdas was used on the masked images, this time providing more accurate results. The isophotes of each galaxy were fitted with a sequence of ellipses with different semi-major axis lengths based on the algorithm described by Jedrzejewski (1987). The final ellipses were then deprojected, using a two-dimensional deprojection procedure assuming that the galaxy was a thin disk, with inclination (given by the axis ratio in the outer parts) and PA inferred from the RC3 catalogue, (similar to the method used in L02). The deprojected fits are hereafter used in our analysis. Although removing the bad areas, before fitting ellipses, is believed to not seriously affect the fitting procedure (Carollo, Stiavelli & Mack 1998), we can be sure that the statistics are not affected by masking bad regions, since both our samples are treated in the same manner.

### 2.1 Sample Selection

We decided to split the sample of bulges at $`T_{RC3}=4`$, since our previous studies have shown that bulges of spirals with types earlier than that have several properties in common with ellipticals, such as surface brightness (SB) profile, age, age spread, and agreement with the fundamental plane of ellipticals (eg. de Jong 1996; Peletier et al. 1999; Falcón-Barroso, Peletier & Balcells 2002). Bulges of late-type spirals have been shown to be different. We require that the galaxies are not too inclined, as this would imply large projection effects, and hide important features. Absolute $`B`$ magnitude, and distance (or $`v_{hel}`$) requirements are imposed, in order to avoid large differences in spatial resolution, and/or discrepant prominence of small scale structures such as stellar clusters or dust lanes. Small scale features are more prominent for nearby galaxies, whereas for more distant galaxies, we would not be able to see much of the inner parts. The same criteria were used when selecting the early-type sample from L02, for which in addition we also rejected all the Seyfert 1 galaxies to avoid contamination of the central regions by the bright nucleus.

An initial sample of all bulges included 70 objects (all late-type), on which we applied further selection. Using parameters from the RC3 catalogue (de Vaucouleurs et al. 1991), we removed all the galaxies more inclined than $`63^{}`$ corresponding to the apparent outer, $`R_{25}`$ disk axis ratio of $`b/a0.45`$. Subsequently, the sample was restricted to a limited radial velocity range of 500 (km s-1) $``$ $`v_{hel}`$$``$ 4000 (km s-1), and absolute B-magnitude of $`M_B15`$ (all inferred from the RC3 catalogue). As we needed to estimate the bulge radius by performing a bulge-disk decomposition (see Sect. 2.2) on the corresponding WFPC2 images, we required availability of WFPC2 optical images on which ellipse fitting was possible. This rules out off-field galaxies (i.e. galaxies for which the central areas were only partly on the field), low SB galaxies and those with very strong star forming regions and/or strong dust features. Although the H-band images are believed to be unaffected by these effects, a few galaxies had to be removed for this purpose. Table 4 quantifies the number of rejected galaxies after application of each selection criterion. Finally, to obtain reliable statistics, it is of vital importance that the two samples are similarly distributed in all salient parameters. We compare our sample with the sample of L02 to make a matching control sample of early type spiral bulges, applying the same ranges in $`M_B`$, distance, and $`i`$, and since not many galaxies with $`0.4<b/a<0.6`$ are available in the sample of L02 we removed 5 galaxies from our late-type sample. Given that the very bright nucleus of type 1 Seyfert galaxies affects the central parts of the images, the Seyferts selected from the L02 sample are all of type 2. The final remaining list of the two subsamples are found in Tables 2 & 3, and the corresponding distributions are illustrated in Fig. 1.

### 2.2 Bulge-Disk Decomposition

A very important quantity for our study is the “bulge radius”. This is defined as the radius at which the SB of the exponential disk equals the surface brightness of the bulge. Successful fits to the bulges of spiral galaxies have been obtained using the Hubble law (Hubble 1930), King model (King 1966), de Vaucouleurs $`r^{\frac{1}{4}}`$ law (de Vaucouleurs 1948), by a generalised version of de Vaucouleurs’ law $`r^{\frac{1}{n}}`$ (Caon, Capaccioli & D’Onofrio 1993; Andredakis et al. 1995) and by an exponential function (Kent et al. 1991; Andredakis & Sanders 1994; Baggett, Baggett & Anderson 1998). Andredakis & Sanders (1994) and Andredakis et al. (1995) showed that for the late-type spirals of this paper an exponential fit of the bulge is entirely adequate, and much better than using an r$`^{\frac{1}{4}}`$ law.

We decomposed the galaxy SB profiles, fitting exponentials to both the bulge and the disk part. Despite the fact that potentially better fits can be obtained with a Sérsic $`r^{\frac{1}{n}}`$ law for the bulge, we decided to use exponential distributions, since these produce the most robust determinations for the bulge radius. As a test we have compared bulge radii obtained using a Sérsic + exponential distribution with bulge radii obtained with double-exponential distributions and found that both values are not very different in log-scale (see Table 4). We find that the measured bulge ellipticities do not change significantly when Sérsic bulge radii are used and our statistical conclusions are not affected. Since the field of view of

NICMOS is too small to be able to perform the Bulge-Disk decomposition on, we used optical WFPC2 images for this purpose. The optically obtained bulge radii do not differ much from those from the H-band images (Carollo et al. 2002). Images with the F555W, F547M, F606W, and F814W filters were used with the galaxy nucleus generally on the Planetary Camera chip. A manual reduction (bad pixel masking, sky subtraction etc.) as described previously was performed on the mosaiced images, followed by ellipse fitting. As a result of the high sensitivity of the optical images to dust extinction, the axis ratio of the successive ellipses was fixed (axis ratio and PA as used for the deprojection) throughout the galaxy. The resulting average radial SB profile is then decomposed into bulge and disk, and the bulge radius is obtained.

## 3 Analysis

The reduced data are now ready to be analysed. In Fig. 7 we plot the deprojected ellipticities and PAs, obtained from the ellipse fitting, as a function of radius (Fig. 7). We used the bulge radii, from the bulge-disk decomposition, to estimate the extent of the bulge region. In this region we quantified the bulge ellipticities in the following way: When the ellipticity profile within the bulge varies by less than 0.1, we define the bulge ellipticity as the mean value. Alternatively, when there is a rise or fall in the ellipticity of more than 0.1, we take the peak ellipticity. These two definitions provide the bulge ellipticity for the whole of our sample (see the histograms in bins of 0.1 in Fig. 2). As we are comparing the results with that of a control sample, any conclusion investigating whether the two datasets are statistically significant needs to be based on a statistical test of the two samples.

We used the Kolmogorov-Smirnov KS-test, which has the advantage of making no initial assumption about the distribution of data, i.e. it is non-parametric and distribution free (Press et al. 1992). For our deprojected ellipticity distributions we obtain a KS probability value of $`0.0160`$, meaning that the two distributions are different at the 98.4% confidence level.

To find out whether this result is real we applied a number of tests. We started by comparing the projected ellipticities of early and late type spirals. The corresponding distributions of the projected bulge ellipticities are given in Fig. 3. The two distributions are peaked at different values (e.g. there is a significant excess of late-type galaxies with bulge ellipticities in the $`0.30.4`$ bin, and a deficit for the $`0.20.3`$ bin). The bulges of late-type spirals are more flattened even when comparing the projected ellipticities, although in this case the KS-test leads to a probability value of $`0.0630`$, corresponding to 93.7% confidence level. The projected ellipticities are consistent with a single distribution, not necessarily implying the same after deprojection.

We continue by defining the bulge ellipticity in a different way, namely such that in all the cases we would take this value to be the average ellipticity of the bulge component of the galaxy. This provides a different distribution (see Fig. 4), with a KS-test result of $`0.3200`$. This implies that the ellipticity distribution of bulges in early and late type spirals is similar, and that the only difference between the two groups is the structure inside the bulge.

Investigating the residual images after fitting ellipses (Fig. 8) shows that in 12 cases the residual images are very smooth all the way to the very centre of the galaxy. This means that the galaxy light is very well described by the set of fitted ellipses This fact shows that the galaxy does not contain very much dust, and that the ellipticity analysis is a very reliable indicator of the distribution of the light. In the remaining 23 cases, the residual images show strong central features meaning that there are significant structures in the central parts. These are disks, bars, strong star formation regions, spiral structure, and/or dust features. In these cases the residual images also show clear indications of dusty bulges for which the inferred morphology is affected by the dust (e.g. the position of the centre changes as a function of radius).

## 4 Discussion

### 4.1 Possible Observational Bias

We have found that the distribution of bulge ellipticities is different between early and late-type spirals (see Fig. 2). The KS-test indicates however that the statistical significance of this result is small. Could it be that we made unjustified assumptions during the deprojection? We discuss here a number of observational biases which could cause this result and find out whether they might be significant.

The difference between the ellipticity distributions might be caused by the inferred 2D deprojection. This procedure assumes that the galaxy is two dimensional, which in reality is not the case. Depending on the thickness variations along the galaxy plane the “deprojected” bulge may seem more elongated as an artifact of the 2D deprojection. If the bulge is much thicker than the underlying disk it becomes more elongated, and consequently if bulges of late-type galaxies are intrinsically thicker than those of early-types, our statistics could be biased. For this purpose we analysed the intrinsic bulge thickness of a sample of 20 edge-on galaxies in the K-band and in the same Hubble-type range as our main sample (Guijarro Roman et al., in preparation). The sample consists of 14 early-type spiral galaxies and 6 late-type spirals, a subset of the sample of de Grijs (1998). For the 20 galaxies we derived the bulge scale heights by directly fitting exponentials to the vertical light profile, and the bulge scale length through B/D decomposition of the radial profile (as described in Sect. 2.2). The ratio of the bulge scale height ($`h_z`$) to the bulge scale length ($`h_r`$) directly measures the bulge axis ratio, and a comparison of these value tells whether the bulge of one galaxy is flatter or thicker than another. Fig. 5 shows that this ratio is not dependent on type, and that bulges of early-type spirals are not intrinsically thicker than bulges of late-type spirals, or vice versa. Hence, the 2D deprojection does not add any bias to our statistics for comparing the deprojected ellipticities.

It could be that the difference is caused by the fact that it is more difficult to detect elongated structures in big bulges of early-type spirals. We can investigate that by looking at the frequency of bars detected inside the bulge regions. Defining bars in the same way as in Knapen et al. (2000) we find no difference in the frequency of bars in the bulge region between early- and late-type galaxies. In our sample, 12 early-type and 13 late-type spirals indicate the presence of a bar within the bulge. These numbers are very similar, and no discrepancy is observed.

Instead of galaxy type we might also look at the bulge to total (B/T) flux ratio, given by

$$\frac{B}{T}=\frac{I_{0B}\times h_B^2}{I_{0B}\times h_B^2+I_{0D}\times h_D^2},$$
(1)
where $`I_{0B}`$ and $`h_B`$ denote the central intensity and the scale length of the bulge, respectively, and the index $`D`$ indicates the same parameters for the disk component.

We derive this parameter for our entire sample and divide the sample in two equal-sized subsamples according to this value, drawing the middle line. This results in two subsamples, one containing galaxies with B/T ratio $`<0.026`$, and one containing the rest of the sample. The bulge peak ellipticities are then compared (see Fig. 6), and the corresponding KS-test probability value is $`0.1975`$. As also seen from the histogram, one cannot say that the distributions are different. Although one might think that the B/T ratio is a good indicator of galaxy type, it is not a very good parameter for our sample here, since most B/T ratios are very small (less than 0.05), and in this range of B/T it is not a very good morphological type indicator.

### 4.2 Implications

We find that bulges of late-type galaxies are more elongated. What does this mean for the formation of bulges? Are we seeing some signs of secular evolution here? Could it be that secular evolution destroys nuclear bars, converting them slowly into bulges, which as a result are still somewhat more elongated?

Theoretical studies have shown that bars are capable of evolving in self dissolving mechanisms (Sellwood & Wilkinson 1993), i.e they undergo secular evolution (Toomre 1966, Hasan & Norman 1990; Pfenniger & Norman 1990). Secular evolution, however, is not always believed to destroy the bars. Merritt & Sellwood (1994) showed that the thickening of bars, through off plane bending instabilities, stops when the density in the mid-plane drops to a low enough value that the natural vertical frequency for a large fraction of particles drops below the forcing frequency below the global bend. There is also some strong theoretical evidence that long lived and strong bars survive loss of up to 2/3 of the angular momentum and reduced pattern speed by a factor of up to 5 (Sellwood & Debattista 1996).

Spiral activity has been shown to occur for a longer time in the outer disk when the initial bar ends at a radius well inside the outer edge of the disk. This spiral activity in the disk outside the bar is not affected by the bar, although the bar grows by trapping additional particles which are ready to lose angular momentum near the inner ends. The new particles added to the bar in this way still have too much angular momentum to sink deep into the bar, and are therefore added to the outer ends of the bar. As a result, the angular momentum of the bar content is rising although the bar itself is slowing down. Bar destruction and pattern speed change are also explained with other theories. Interaction with the halo (Weinberg 1985, Debattista & Sellwood 1996, Athanassoula 1996) is believed to be able to slow down a bar through dynamical friction. Mergers (Gerin, Combes & Athanassoula 1990; Athanassoula 2002; Heller & Shlosman 1994; Barnes & Hernquist 1992), or fuelling of an active nucleus by driving gas towards the centre, play a considerable role in destruction of bars resulting in a spheroidal bulge (Norman, Sellwood & Hasan 1996).

It is not obvious that secular evolution of bars leads to oblate structures. Secularly evolved bars may also end up as elongated structures, as bars are believed to be populated by orbits of greater eccentricity than those in the axisymmetric parts of the disk. One might think that bar formation, on the other hand, is expected to be less common in early-type systems (Kormendy 1979, 1993; Combes & Elmegreen 1993, Sellwood & Wilkinson 1993; Lütticke, Dettmar & Pohlen 2000), as there is less gas in early-type spirals, and there is also already a big bulge which stabilises the galaxy and prevents self-gravitating bar formation. If secular evolution is active one would therefore expect to find more elongated structures in late-type spirals, as we find observationally. This argument, however, might not be correct, since recent statistics on bars (e.g. Lütticke et al. 2000) shows that the bar fraction in early type spirals (Sa’s) is the same as in Sc’s. Also, Athanassoula (2002) has recently shown that massive bulges do not always suppress bar formation, when more accurate simulations with live-particles bulges/halos are done. In any case, in this paper we do not see a big difference between shapes of central features in early and late type spirals, so the process which makes elongated bulges in late-type spirals is not much more efficient than in early-type spirals.

To summarise, secular evolution seems to be able to make bulges (central disk concentrations) from bars, but it is not clear what their final shape is, and how efficient this process is. Yet so far it is not possible to confirm nor refute any of these theories, as very few observational studies have been done (Gerssen et al. 1999, Zimmer & Rand 2002), and only improved observations will be able to tell whether secular evolution of bulges is indeed possible.

## 5 Summary

Studying a sample of 70 galaxies in the H-band, we find that inner features in bulges of late-type spirals are more elongated than those in early-type spirals, although the statistical significance of this result is small ($`P_{KS}=0.0160`$). When we simply compare the average ellipticity in the bulge between both samples, we don’t find any difference. This probably indicates that bulges of later type spiral galaxies contain more elongated features like nucluear bars than bulges of early-type spirals. We have performed several tests to establish that this result is not due to observational effects: bars are visible just as easily in early-type bulges as in bulges of late type galaxies. Also, deprojection cannot cause the larger fraction of elongated features in late-type spirals.

The result could be explained if bulges of late-type spiral galaxies are formed primarily through secular evolution of bars, while this would not be the case for earlier-type bulges. Since however it is still unclear what the morphology of bulges that have been created through secular evolution, it is not possible at this stage to give more detailed conclusions. To better understand the formation process of bulges, it is important that we first obtain a better understanding of the process of secular evolution.

## Acknowledgement

We thank Michael Merrifield for very interesting and useful discussions. The data presented in this paper were obtained from the Multimission Archive at the Space Telescope Science Institute (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NAG5-7584 and by other grants and contracts. Much of the analysis was performed using Iraf which is distributed by NOAO.