text:
# Observation of electronic and atomic shell effects in gold nanowires

## Abstract

The formation of gold nanowires in vacuum at room temperature reveals a periodic spectrum of exceptionally stable diameters. This is identified as shell structure similar to that which was recently discovered for alkali metals at low temperatures. The gold nanowires present two competing ‘magic’ series of stable diameters, one governed by electronic structure and the other by the atomic packing.

preprint: 7

When the conductance of metallic nanowires can be described in terms of a finite number of quantum modes it is expected that the degree of filling of the quantum modes has a measurable effect on the total free energy of the nanowires van Ruitenbeek et al. (1997); Stafford et al. (1997); Yannouleas and Landman (1997). Indeed, Yanson et al. Yanson et al. (1999) observed that alkali nanowires, obtained by the mechanically controllable break junction (MCBJ) technique, exhibit exceptional stability for certain diameters. The series of stable diameters was shown to correspond directly to the well-known magic number series of alkali metal clusters produced in vapor jets (for a review see de Heer (1993)). For alkali metal clusters and nanowires the main features of this electronic shell effect can be understood in terms of the filling of simple free-electron modes inside a spherical or cylindrical cavity, respectively. The electronic level spectrum shows groups of levels bunched together and a cluster or nanowire is favored when such a shell is just filled. The principle of this shell structure is well known from atomic physics, giving rise to the periodic table of the elements, and from nuclear physics, explaining the periodic appearance of stable nuclei. Apart from this electronic shell structure alkali metal clusters and nanowires show, typically at larger diameters, a second series of stable diameters resulting from the atomic structure. Completed atomic layers stabilize the nanowires.

The alkali nanowires are very attractive from a fundamental point of view, since the electronic structure is well-described in terms of free and independent electrons, but they are not very suitable for possible applications. As suggested by research on metal clusters shell effects for nanowires should be observable for many other metals. As a first step beyond the alkali metals we have chosen gold, since the physics is expected to be very similar, it being to a good approximation still a free-electron metal.

Recently some experimental results on the formation and stability of gold nanowires have already been reported. By the use of an electron beam thinning technique in a transmission electron microscope under ultra-high vacuum at $`T=300`$ K Kondo and Takayanagi Kondo and Takayanagi (2000) observed gold nanowires that adopt a series of unusual helical configurations, which is an exceptional form of atomic shell structure. Using STM indentation experiments Medina et al. Medina et al. (2003); Díaz et al. (2003) reported the observation of atomic shell effects in gold nanowires. Here, we demonstrate clear evidence for electronic shell effects, with a cross-over to atomic shell structure at larger diameters.

The atoms need to have sufficient mobility in order to be able to explore many wire configurations and find the stable ones corresponding to shell filling. Therefore, in order to observe shell effects we need to raise the temperature to a sizeable fraction of the bulk melting temperature. For the alkali metals a temperature near 80 K was sufficient and the experiments were performed with a standard low-temperature MCBJ device with a local heating extension Yanson et al. (1999). For gold this is not sufficient and we have developed a new MCBJ instrument (Fig. 1) operating under ultra high vacuum (UHV) allowing the temperature to be varied between 70 K and 500 K. A notched gold wire is fixed with stainless steel clamps, bridging two bendable metal substrates. The wire is broken at the notch in a controlled way by mechanically applying pressure from below onto both substrates. Once broken the bending beams can be relaxed to bring the wire ends back into contact and atomic-sized contacts can be finely adjusted using a piezoelectric element. The base pressure is $`5\times 10^{10}`$ mbar.

Nanowires are obtained by indenting the two clean fracture surfaces and subsequently stretching the contact while monitoring the conductance of the wire. The conductance trace follows a sequence of plateaus separated by abrupt jumps. At the final stages the contact consists of just a few atoms and the conductance is changing in steps signaling the atomic structural rearrangements of the contact. For monovalent metals the conductance of a single-atom contact is close to 1 $`G_0=2e^2/h`$, and carried by a single conductance channel Scheer et al. (1998). The conductance is measured at constant bias voltage, recording the current through a current-voltage convertor and a 16 bit analog-to-digital convertor. A trace from 100 $`G_0`$ to 0 is recorded in about 10 s at a sampling rate of 10 ms.

At larger diameters the conductance jumps can be larger and each conductance trace is different depending on the evolution of atomic arrangements in the nanowire. In order to reveal underlaying variations in stability we perform a simple statistical analysis by adding all values of the digitized conductance traces into conductance histograms. Preferred nanowire structures are then seen as maxima in the conductance histograms. The wire radius $`R`$ corresponding to a given conductance value can be obtained from a semi-classical expression Torres et al. (1994); Höppler and Zwerger (1998):

$$GG_0\left[\left(\frac{k_\mathrm{F}R}{2}\right)^2\frac{k_\mathrm{F}R}{2}+\frac{1}{6}\right],$$
(1)
with $`k_\mathrm{F}`$ the Fermi wave vector.

Figure 2 shows a conductance histogram up to a conductance of 60 $`G_0`$ obtained for gold at room temperature (bold curve). The histogram is five-point smoothed by using a second order polynomial. The data show strong resemblance to the histogram for Na recorded at $`T=80`$ K by Yanson et al. Yanson et al. (1999), which is shown for comparison in the same graph. The high counts at the peaks result from stable wire diameters. The peaks between 7 and 25 $`G_0`$ agree very well. Above 25 $`G_0`$ only one or two more peaks are visible for Au, much less than for Na, and they may be shifted somewhat to lower conductance. The first few peaks below 7 $`G_0`$ do not correspond very well. The low-conductance histograms for Na have been shown before to differ significantly from those for Au Brandbyge et al. (1995); Krans et al. (1995). The direct correspondence in period and phase of the two series of peaks at higher $`G`$ is strong evidence for similar electronic shell structure in Au and Na.

For shell effects the peaks are expected to be equidistant when plotted as a function of the radius of the wire. This is verified in the inset to Fig. 2, where we plot $`k_\mathrm{F}R`$ obtained from Eq. (1) as a function of peak index. The peaks have been reproduced for different samples in six independent series of measurements on Au and in order to illustrate the degree of reproducibility the plot contains points from three measurements. The period $`\mathrm{\Delta }k_\mathrm{F}R`$ shows some variation between the individual measurements that averages as $`\mathrm{\Delta }k_\mathrm{F}R=1.02\pm 0.04`$. This is somewhat smaller than the slopes for electronic shell structure observed for the alkali metals Yanson (2001).

Often we observe another periodic structure in gold conductance histograms recorded under similar experimental conditions, as illustrated by two examples in the Fig. 3. Many histograms show a cross-over between the electronic shell period of Fig. 2 to a new shorter period. The peaks were identified using a software tool (symbols) and Fig. 4 shows the radii corresponding to the peaks in Fig. 3(a). The cross-over is at the second point, above which we find a period $`\mathrm{\Delta }k_\mathrm{F}R=0.40\pm 0.01`$. This structure has been reproduced in well over ten independent measurements and a second example is shown in Fig. 3(b) for comparison. It can be seen that the period $`\mathrm{\Delta }k_\mathrm{F}R`$ displays some variation, giving an overall average of $`\mathrm{\Delta }k_\mathrm{F}R=0.40\pm 0.03`$. Again, this second periodic structure agrees closely to the results obtained for alkali nanowires and by analogy we interpret it as an atomic shell effect.

Thus we attribute the peaks in the histograms in Figs. 2 and 3 to two different and independent effects: electronic and atomic shell effects. The electronic shell effect gives rise to peaks in the histograms because specific radii of the nanowire are more stable than average due to local minima of the electronic energy of the system. Electronic shell effects have been previously found for clusters of gold atoms in vapor jets Katakuse et al. (1985). The magic number series of stable clusters has been successfully explained by approximating the electronic structure of gold metal by homogeneous free-electron (jellium) models de Heer (1993); Brack (1993) confined by a spherically symmetric potential well. In the same spirit we model our nanowires as a free electron gas confined by a cylindrical potential well Stafford et al. (1997); Yannouleas and Landman (1997); Höppler and Zwerger (1998); Ogando et al. (2002). The magic diameters are obtained by a semiclassical approach using the Bohr-Sommerfeld quantization condition. Stable trajectories for the electron wavepackets are periodic orbits inside a cylinder, that can be classified as diametric, triangular, square, etc. Each of these trajectories gives rise to a periodic series of oscillations in the density of states as a function of $`k_\mathrm{F}R`$, that produce oscillations in the total energy of the system. The dominant period of the oscillation depends somewhat on the electron density and the shape of the potential well. We obtain a number of $`\mathrm{\Delta }k_\mathrm{F}R=1.21\pm 0.02`$ using an electron density corresponding to that for gold from the calculations by Ogando et al. Ogando et al. (2002), which agrees well with the periodicity of shell structure in gold clusters, $`\mathrm{\Delta }k_\mathrm{F}R=1.23\pm 0.01`$. The differences in the peaks observed for Na and Au at low $`G`$ can be largely explained by considering elliptically distorted wires and including differences in surface tension for the two metals Urban et al. (2003). However, we find a reproducible peak at 4 $`G_0`$ for gold, which does not appear in this analysis and may be due to quadrupolar distortions (Urban et al. (2003) and private communications).

The period obtained experimentally for gold nanowire magic radii, $`\mathrm{\Delta }k_\mathrm{F}R=1.02\pm 0.04`$, is slightly smaller than that expected from the semiclassical analysis. Deviations similar to the present results for Au have been observed for Li Yanson (2001). Such deviations have been attributed to backscattering of the electrons on defects, which is not included in the idealized semiclassical model. Defect scattering may also explain the observed variation of the slopes between experiments and the fact that the range of observed oscillations for Au is smaller than for Na, since scattering will lead to smearing of the peaks. In order to preserve a constant slope the correction should be proportional to the conductance $`G`$ itself. In other words, a certain fraction of the electrons is scattered back, independent of the wire radius. This would for instance be the case when the dominant mechanism is scattering by roughness on the nanowire surface, which seems plausible.

The Fermi surface of bulk gold has marked deviations from spherical symmetry. The nearly spherical main sheets of Fermi surface are connected by necks across the Brillouin zone boundaries. The necks may introduce an additional set of oscillations in the density of states. However, since the wavelength of the states in the neck is about six times larger than the main Fermi wavelength the period would be much larger than that due to the main states. Moreover, the number of states in the neck is small making the amplitude of the effect small. Thus we expect that it is a good approximation to consider gold as free electron metal having the same electronic shell effect as in alkali metals. This is supported by electronic structure calculations Opitz et al. (2002) for the quantum modes in nanowires of Na and Cu (which has a Fermi surface similar to Au).

For larger radius $`R`$ the configuration energy (atomic shell structure) becomes more important than the energy contribution due to the conduction electrons. The amplitude of the oscillations in the electronic free energy scale as $`1/R`$ Yannouleas et al. (1998), while the ones for the surface energy increase proportional to $`R`$, making the latter dominant at large radii. The cross-over point is seen to vary between experiments, which is likely due to differences in the local crystal orientation of the leads connecting the nanowire. Note that the histogram in Fig. 2 shows some additional fine structure that may be due to an admixture of atomic shell structure.

The  orientation has been shown to be most favorable to form nicely faceted long nanowires Rodrigues et al. (2000); Jagla and Tosatti (2001). At large diameters we expect the nanowires to order into densely packed wires, that minimize the surface energy. We assume that the bulk fcc packing of gold is preserved, which is supported by the observation of gold nanowires in transmission electron microscopy by Rodrigues et al. Rodrigues et al. (2000). The lowest energy surfaces for gold are perpendicular to $`[111]`$ and we can construct a densely packed wire along a $`[110]`$ axis with four $`(111)`$ and two $`(100)`$ facets, as also proposed by Ref. Jagla and Tosatti, 2001. Completing a full atomic shell gives rise to a rather long period of $`\mathrm{\Delta }k_\mathrm{F}R=2.85`$, but assuming stable wires are found when completing a single facet we obtain a six times smaller period $`\mathrm{\Delta }k_\mathrm{F}R=0.476`$ Yanson et al. (2001). The faceting structure of stable wires can be recognized in a Monte Carlo simulation of the thinning-down of a fcc wire by Jagla and Tosatti Jagla and Tosatti (2001). The experimental slope is about 20% smaller than expected, but this is consistent with a similar deviation observed for the electronic shell effect and is similarly attributed to defect scattering.

Our results and interpretation differ in several essential points from the work in Ref. Díaz et al., 2003, where the observed structure is entirely attributed to atomic shell structure. However, the well-developed peak numbers 3, 4 and 5 in Fig. 1 of Ref. Díaz et al., 2003 agree with the first three peaks of what we interpret as electronic shell structure in Fig. 2. We have no need to resort to a second-derivative analysis, as was done in the cited work to analyze the histograms at higher $`G`$, which entails the risk of identifying residual fluctuations in the background as peaks. Also, our peak structure shows the expected global decrease of the amplitude towards larger diameters, in contrast to the data by Medina et al. Medina et al. (2003); Díaz et al. (2003). We conclude that electronic shell structure is clearly observed in Au nanowires, with a cross over to atomic shells at larger diameters.

Apart from atomic shell structure based on conventional fcc lattice packing one should consider the ‘weird wire’ configurations that have been predicted Gülseren et al. (1998) and observed for gold in transmission electron microscopy Kondo and Takayanagi (2000). The helical atomic wire arrangements lead to a series of ‘shells’ for which we estimate the conductances shown by bars in Fig. 3. In view of the much more regular pattern observed we conclude that under conditions used in our experiment the regular lattice atomic shell structure is predominant. This may be a consequence of the path by which the nanowire evolves from a large fcc wire down. Making a structural cross-over to helical arrangement along the way presumably involves a relatively large barrier.

We thank R. van Egmond for valuable assistance with the experiments, and J. Bürki, H. Grabert, and D. Urban. This work is part of the research program of the “Stichting FOM,”and was supported by the EU TMR Network program DIENOW.


format: text

text:
# Possibility of a Metallic Field–Effect Transistor

## Abstract

We develop theoretical arguments that demonstrate the possibility of metallic field–effect transistors (METFET’s) in one–dimensional systems and particularly in armchair carbon nanotubes. A very inhomogeneous electric field, such as the field of a tunnelling tip, can penetrate the relatively weakly screened nanotubes and open an energy gap. As a consequence, an energy barrier forms that impedes electron flow and thus permits transistor action. This type of metallic field effect is advantageous because of the high conductance of the metallic tubes in the ON–state.

theory, metal FET, nanotube, nanodevice, electronic structure

Field–effect transistors in current use are semiconductor devices. The scaling trend to nanometer dimensions calls for ever higher doping and channel conductance of these devices hessbook . Ultimately one desires a conductance close to that of a metal if one wants to scale devices to the smallest possible size. However, a metallic conductance also prevents the penetration of the electric field except for extremely short distances; ordinarily too short to achieve device function.

We propose here an innovative approach to control electron transport in metallic one–dimensional (1D) systems by use of the inhomogeneous electric field of a low–dimensional highly localized ”modulating” gate (MG) such as nanometer tips or metal nano-interconnects. Use of a highly localized gate results in a strong enhancement of the electric field in a narrow region. Thus any depletion of charge is easier to achieve since it occurs only in extremely small volumes. In addition, the weak screening of electric potentials in 1D channels further enhances possible field effects.

The purpose of this letter is to show that metallic carbon nanotubes may be suitable for use as metallic field–effect transistors (METFET’s). Key to this novel metallic field effect is the opening of a band gap due to the breaking of the mirror symmetry by the MG. Note, that a region of spatially non-uniform electric field also represents a barrier that reflects electrons, while the homogeneous field of an extended gate changes predominantly the carrier density. These qualitative factors are discussed below in a more quantitative way and their potential for device applications is illustrated.

Metallic field effect in Armchair NT’s. – IBM researchers have enunciated the vision of combining metallic and semiconducting single–wall nanotubes (M–SWNT’s and S–SWNT’s) in circuits appenz-utah . The M–SWNT’s would serve as interconnects and the S–SWNT’s as active devices of the extremely small size (typically, the SWNT radius is $`0.7`$ nm) appenz2 ; appenz ; avouris-r ; avouris . Here we discuss the possibility to also use metallic tubes as transistors or switching devices, which is appealing because of their virtually ballistic conductance pseudos ; lundstr ; appenz2 .

The band structure of an armchair SWNT, as shown in Fig.1, has two subbands A and B (shown as bold green and blue curves) crossing at the Fermi level. Absence of the gap is ultimately related to the SWNT mirror symmetry damnjan , which forbids any mixing of these subbands. We propose to break this mirror symmetry by applying very non–uniform electric fields by use of a MG. Two possible mechanisms for gap opening in armchair SWNTs are: (1) by a direct mixing of electron and hole states of the A and B subbands, which happens if the external potential has an atomic scale corrugation, and (2) by the indirect subband mixing in a higher order of the perturbation theory, which essentially can be realized for any non–uniform potential. Both mechanisms require (i) a relatively high electric field (order of $`10^7`$ V/cm) and (ii) SWNT symmetry breaking.

Group theory proves that the direct matrix element is non–zero if the Fourier transform of the potential $`U_{q,\delta m}`$ is a full scalar with respect to the rotations of the symmetry group of the nanotube damnjan . This is fulfilled for $`\delta m=sn`$, where $`s=2,4,\mathrm{}`$ is an positive even integer and $`n`$ is the integer appearing in the notation \[n,n\] for armchair SWNT’s newyan . Then, the energy dispersion of two new subbands $`|A\pm B`$ in vicinity of the Fermi point reads as: $`E_\pm =\pm \sqrt{E_A(q)^2+U_{q,\delta m}^2},`$. Here $`E_A(q)=E_B(q)=E_B(q)`$ is the dispersion in A/B subband at zero gate voltage. At the Fermi point $`E_A(0)=0`$ and, thus, a gap opens:

$$E_g^{(sn)}=2|U_{q,sn}|$$
(1)
This is a direct mixing of the subbands A and B. The bandgap is linear in the applied gate potential (see Inset in Figure 1). A full calculation presenting a selfconsistent solution of joint Schroedinger and Poisson equations following the derivation of Refs. yanli ; jetpl  will be published elsewhere. Except for a natural condition that the field must not cause an electric breakdown, there exists no upper limit for the magnitude of the gap that opens in this case and we estimate the possibility of a gap of several eV.

Thus direct mixing of the crossing subbands allows to open a large gap. It does requires, however, a very high multipole moment of the potential in order to break the symmetry (we must at least have $`\delta m=2n`$) and this is difficult to achieve because the amplitude of any high multipole component will be typically much smaller than the gate voltage, $`U_{q,2n}<<V_g`$. Nevertheless, it may be possible to create such a potential by applying the gate voltage to a chemically modified surface of the SWNT. This approach is technologically complex.

Fortunately, there exists a more straightforward possibility to open a bandgap by virtually any kind of inhomogeneous electric field that breaks the mirror symmetry. A bandgap opens also by indirect mixing of the subbands. This bandgap opening is smaller because it arises from higher order perturbation terms (similar to a nonlinear Stark effect). In particular, any potential that contains both a uniform field and a higher order (quadrupole) component will lead to a gap which scales with the third power of the gate voltage. In what follows we consider an example of two metal tips at opposite sidewalls of the SWNT: one at $`+2V_g`$ and the other at $`V_g`$ bias. The dependence of the bandgap on $`V_g`$ is shown for  and  armchair SWNTs in the Inset of Fig.3. The gap $`E_g`$ first increases with $`V_g`$ and decreases beyond $`\pi \gamma /(en)`$ where $`\gamma 3`$ eV is a hopping integral. Perturbation theory yanli ; newyan ; nanotech  predicts a maximum gap that depends only on the size of the tube and scales in a universal way for nanotubes of arbitrary chirality:

$$E_g^{(c)}\frac{\mathrm{}v_Fb}{2R^2}\frac{\gamma }{n^2},$$
(2)
This was confirmed by numerical calculations in Refs.yanli ; nanotech . Here $`v_F=3\gamma b/2\mathrm{}`$ is the Fermi velocity for a M–SWNT, and $`b0.14`$ nm is the bond length. Even though the maximum gap is small in this case and thus has consequences for the conductance only at low temperatures, the gap opening should be much easier to achieve experimentally e.g. by use of the inhomogeneous electric field of a tunnelling tip (Fig.2). Also, ultra narrow leads rotkina , fabricated closely to the nanotube channel, special (electro–)chemical function groups at the tube sidewalls or inside the tube zharov  may be used as a MG. We notice that use of a dual gate (both local and backgate) may be beneficial for 1D METFET’s because of the uniform backgate controls the Fermi level (charge density) while the MG controls the conductivity of the channel.

Field effect modulation of metallic conductance. – Any opening of the gate induced semiconductor gap along the armchair metallic tube will create a potential barrier for the electrons and therefore modulate conduction. (Note that this modulation may be enhanced by Coulomb blockade ftft .) Inclusion of tunnelling is important for the calculation of the current because the effective mass of the electrons in SWNT’s is very small ($`0.06m_o`$, similar to the in–plane effective mass of graphite ieeenano ). At non zero temperature thermionic emission must be taken into account as well. Although the rate of tunnelling through the classically forbidden (gated) region is high, a significant metallic field effect can be achieved by increasing the width of the gated region, or by operating the device at low source–drain voltage, $`V_d`$, and low temperature, $`T`$. (The ON/OFF current ratio is controlled by the parameter $`E_g/V_d`$.) Because the transport in the armchair SWNT METFET is ballistic, its conductance in the ON state is limited only by reflections at the contact. We will use in what follows (mainly for normalization purposes) a maximum conductance of $`4G_o=2e^2/h`$ ftn .

For the OFF state, we estimate a semiclassical WKB transmission coefficient, $`T(E)`$, by assuming a uniform gap of width $`W`$ in the gated region. In order to extend the validity of our results to non–zero temperatures we calculate the total transmission of the channel by integrating the partial current, $`T(E)`$, multiplied by the difference of the Fermi distribution functions of the left and right electrodes.

Figure 3 shows typical IV curves (IVC) for a METFET using a  armchair M–SWNT of the diameter 0.7 nm at $`T=`$4K, $`W=50`$ nm. The upper (red) IVC corresponds to a zero gate voltage (no gap). The channel is fully open and the ON current is determined by injection from a contact and thus by the quantum conductance $`4G_o`$ ftn . With increasing gate voltage one observes a substantial decrease of the current due to the opening of the gap and a depletion of electrons in this region (blue curve). In the Right Inset of Fig.4, the METFET current is plotted vs. the $`V_g`$ for the dipole–quadrupole potential (two tip geometry) as described above. At low drain bias (blue curve) the METFET is switched OFF at $`V_g3`$ V.

The total transmission coefficient of the METFET, $`T(E)=G(E)/G_o`$, assuming transparent contacts, is plotted in the Left Inset of Fig.4. Several families of curves are given for different temperatures. In the upper family of curves (room temperature) the suppression of transmission through the gated region is clearly seen as a function of an increased gap: from 0 to 0.12 eV (from green to red, top to bottom). At lower temperatures $`T=`$77K and 4K this effect is naturally sharper.

The current is an exponential function of the gap (and, hence, of the gate voltage). Results are shown in Figure 4: by opening a gap of about 0.5 eV (which would certainly require sophisticated technology) one can decrease the current by five orders of magnitude.

In conclusion, we have proposed a novel type of electronic switching device based on carbon nanotubes: a one–dimensional metallic field–effect transistor (1D METFET). We have described armchair SWNT’s in detail because of their special symmetry that also gives rise to virtually ballistic transport. Our calculations demonstrate, at least principle, the possibility to open a band gap by application of inhomogeneous electric fields that are created by special gates (modulation gates) that in the simplest case resemble a tunnelling tip to break the SWNT symmetry. Assuming ballistic transport, we have calculated IV curves for the METFET’s with gate widths of the order of 15–50 nm and found significant modulation as well as reasonably large ON/OFF current ratios. We finally note that one can expect excellent scaling properties of the SWNT METFET because of the metallic conductivity enhanced by ballistic transport.

Authors are indebted to Mr. B. Grosser for help with image preparation and Ms. Y. Li for help in numerical computations. This work was supported by NSF Grant No. 9809520, by the Office of Naval Research grant NO0014-98-1-0604 and the Army Research Office grant DAAG55-09-1-0306. SVR acknowledges partial support of DoE Grant No. DE-FG02-01ER45932, and NSF grant No. ECS–0210495.


text:
# Cosmology and the Stellar Halo

## 1. Introduction

The LCDM model of structure formation predicts that the dark matter halos of galaxies like the Milky Way were formed by the continual accretion of lower-mass systems. While some fraction of the accreted dark halos are expected to survive as subhalos, most of the extended dark halo is built up from merging systems that are tidally destroyed (Zentner & Bullock 2003). Similarly, if the accreted systems contain dissipationless stars, e.g. in the form of a dwarf galaxy, then this process leads naturally to the formation of an extended stellar halo. The LCDM expectation for the formation of the stellar halo is reminiscent of the Searle & Zinn (1978) scenario of chaotic accretion, and has the positive feature that the expected accretion history is well-known theoretically, at least for the dark matter. In this contribution we describe efforts to model the formation of the stellar halo in the context LCDM cosmology using a set of idealized N-body simulations. We make predictions for the abundance of stellar halo substructure as a function of radius and argue that deep probes of the stellar halo offer a unique test of the hierarchical nature of structure formation on small scales.

Modeling the stellar halo from a cosmological context is non-trivial. Fully self-consistent cosmological simulations allow the dark matter dynamics to be followed in detail, but even simulations of millions of particles will only have a few thousand particles in the very largest satellites (e.g. Helmi, White, & Springel 2002). Moreover, such simulations are computationally intensive and hence prohibit examining more than a handful of halos. Analytic or semi-analytic approaches (e.g. Bullock et al. 2001, Johnston, Sackett & Bullock 2001) allow the production of many halos, but only approximate the dynamics.

In this study we take a hybrid approach, and use accretion histories generated semi-analytically within a cosmological model to inform us of the luminosity, mass and accretion times of satellites, as well as the instantaneous mass of the parent galaxy. An N-body simulation for each accreted satellite is then run within an analytic model for the parent. This approach allows a vast decrease in the necessary computing time, while at the same time allows us to follow the dynamics in detail. Our approach is limited in its ability to correctly follow major mergers, so our model results are most accurate for the outer halo ($`r>20`$kpc). The inner halo formed early on ($`t>8`$Gyr ago) during a period of rapid mass accretion when major mergers were common. In addition, the detailed phase-space structure of the inner halo was likely influenced by the formation process of the Galaxy itself, which is rather poorly understood. For this reason we focus mostly on the outer halo, which should have formed from more recently-accreted systems ($`t<8`$Gyr). Fortunately, the vast majority of expected stellar substructure comes from accretion events that occurred late, after the last $`>10\%`$ merger event.

We describe the method more fully in §2 and our results and conclusions in §3. In what follows we assume an LCDM cosmology with $`\mathrm{\Omega }_m=1\mathrm{\Omega }_\mathrm{\Lambda }=0.3`$, $`h=0.7`$, $`n=1`$, and $`\sigma _8=0.9`$. A more thorough presentation of this work will be given elsewhere (Johnston & Bullock; Bullock & Johnston, in preparation).

## 2. Methods

### 2.1. Cosmological Context and Semi-Analytic Star Formation

We set the mass accretion and subhalo accretion histories of our host halos using Extended-Press-Schechter (e.g. Lacey & Cole 1993) and employ the merger tree algorithm proposed by Somerville & Kolatt (1999). This allows us to generate mass accretion histories of our host halos as a function of time, as shown in Figure 1 for four example realizations (upper solid lines). We demand that the hosts at $`z=0`$ have a mass of $`1.4\times 10^{12}M_{}`$. The merger tree also provides a list of the masses, $`M_\mathrm{s}`$, and accretion redshifts, $`z_\mathrm{m}`$, of all subhalos that merge to form the host halo.

Assigning luminosities and stellar masses to each accreted dark matter halo is a fundamental step in any attempt to model the stellar halo in this context. We do so by modeling star formation in each system based on an estimate for its total mass in cold gas: $`\dot{M_{}}=M_{\mathrm{gas}}/t_{}`$, where $`t_{}`$ is the star formation timescale. For massive systems, the mass in cold gas accreted in a given time is proportional to the dark matter mass accreted: $`M_{\mathrm{gas}}=f_b(1f_{\mathrm{hot}})M_{\mathrm{dm}}`$. Here $`f_b=0.13`$ is the baryon fraction, and $`f_{\mathrm{hot}}<1`$ parameterizes the amount of hot gas that might arise from inefficient cooling and feedback. For low-mass systems ($`<30\mathrm{k}\mathrm{m}`$s-1), we solve the dwarf satellite problem (e.g. Klypin et al. 1999) following the suggestion of Bullock, Kravtsov & Weinberg (2000), and assume that the total amount of gas available to form stars is proportional to the mass in place at the epoch of reionization $`z_{re}=10`$. Remarkably, we find that the choices $`t_{}=15`$ Gyr and $`f_{\mathrm{hot}}=0.85`$ allow us to roughly reproduce the main characteristics of our Galaxy, including the stellar halo, disk, and bulge masses, as well as the observed star formation rates and gas mass fractions of isolated dwarf galaxies and Local Group irregulars (see Johnston & Bullock, in preparation).

Crucial in this agreement is the rather long star formation timescale we have assumed for low-mass systems, $`t_{}=15`$Gyr. This means that systems that are accreted early tend to have high gas mass fractions, and contribute more to the main galaxy than the to stellar halo. Indeed, the observed ongoing star formation in many dwarf galaxies effectively demands that the star formation timescales must be long, at least in an average sense (e.g. Grebel 2001). The star formation timescale is likely shorter in larger-mass systems like the Milky Way host, perhaps because the gas can reach higher densities and feedback effects are less important in stabilizing the star formation. Figure 1 shows the expected evolution of the disk, bulge, and stellar halo components for four example accretion histories (dotted, dashed, and thick-dashed from top to bottom in each panel). We have associated mass accreted in hot gas with disk material, cold gas associated with disrupted satellites as bulge material, and unbound stellar mass as stellar halo material. While the disk mass and bulge mass build-up presented in this figure is somewhat dependent on our assumptions about gas infall, the general feasibility of our picture for how galaxies like the Milky Way were assembled is clear. We go on to use this simple semi-analytic model to generate the input conditions for our more detailed numerical modeling.

### 2.2. The parent galaxy

The potential of the parent galaxy is represented by three analytic functions, corresponding to bulge, disk, and halo components. The first two are held at fixed mass and scale. For the halo, we use a spherical NFW (Navarro, Frenk & White, 1996) model that grows in mass and scale according to the smoothed accretion history of the parent galaxy, described above. We do so following the prescription developed by Wechsler et al (2002) to fit fully self-consistent N-body simulations.

### 2.3. The satellite galaxies

We find that over course of its history, each Milky-Way-type halo accretes $`100`$ luminous satellites. We simulate the orbital evolution of each satellite once accreted using 105 particles, initially distributed as an NFW profile with mass, scale and accretion time taken from the semi-analytically generated accretion history. The orbits of the satellites are chosen at random so that the distribution of final orbits matches that measured by Ghigna et al (1998) at the end point of a cosmological N-body simulation. The self-gravity of the satellite is calculated using a basis-function-expansion code (Hernquist & Ostriker 1992) and the influence of the parent is imposed from the analytic functions described above. In addition, all particles within two tidal radii of the satellite feel a drag force due to dynamical friction, calculated as suggested by Hashimoto, Funato & Makino’s (2003) simulations.

The total mass in luminous matter in the satellites is assigned based on the star formation algorithm described above (gas mass is ignored for now.) Based on its luminosity, we assign each satellite galaxy a King density profile, with a core radius set to match the distribution of scales observed in local dwarf galaxies (Dekel & Woo, 2003). We assign a luminosity to every particle based on its energy at the time of accretion. We do so using Eddington’s formula (Binney & Tremaine 1982) for the appropriate King profile embedded within an equilibrium NFW halo. The key element here is that the star particles tend to be much more tightly bound than than most of the dark matter in the infalling satellite halos.

## 3. Preliminary Results

As shown by the thick dashed lines in Figure 1, the total stellar halo mass at $`z=0`$ in our model matches well what is expected for our Galaxy, $`2\times 10^9M_{}`$, with the total value varying from galaxy-to-galaxy by a factor of $`2`$. The top panel of Figure 2 shows the resulting stellar halo profiles derived from our simulations for four example halos (four thick lines, see caption). The thin line, for reference, shows an $`r^3`$ power law. Rough agreement with the expected power-law is evident, although we do predict a steep fall-off beyond $`50`$kpc. Some understanding of this steepening can be gained from the lower panel of Figure 2, which shows our derived stellar halo profile along with the profile of all of the unbound dark matter mass in our simulations. We see that the stellar profile falls off more quickly at large radii because the star particles tended to sample the most tightly-bound cores of the infalling halos. It is those cores that tend to survive deep into the halo core before becoming disrupted.

Figures 3 and 4 show Aitoff sky projections for two of our stellar halo realizations. Projected in radial bins, each point represents a tracer star that samples 1 in 500 $`L_{}`$, and we have assumed a stellar mass-to-light ratio of 2 in the V-band. The outer halos, beyond $`50`$kpc, show pronounced substructure, while the inner halos are more smooth. We expect structures of this kind to be uncovered by deep halo surveys like the Deep Giant Star Survey, Spaghetti Survey, and the Sloan Digital Sky Survey (Majewski et al. 2000, Morrison et al. 2000, Yanny et al. 2000). We have colored each point based on the time it became unbound (scale bar at the top), and the substructure is even more evident in this representation. Specifically, more recent disruption events tend to be seen as coherent substructure. This is dramatically apparent in the great-circle stream in Figure 3 at $`r3060`$kpc, but smaller coherent streams can be seen as far in as $`15`$kpc, especially in Figure 4. If star formation is ongoing in satellites after they are accreted, then more chemically evolved stars will be associated with more recently unbound particles. This suggests that chemical probes may prove useful in substructure searches. While not shown, similar trends are seen even if we tag particles simply by the time when the satellite was accreted.

Note that the surviving satellites (boxes) are likely to have had significantly different chemical evolution histories than those that disrupted early and contributed to the central halo. The median accretion time for surviving satellites is $`5`$Gyr in the past, while most of the stellar halo was  in place $`>8`$Gyr ago. That the surviving systems are expected to represent a different population of dwarfs than those that get destroyed may help explain why most satellite dwarf abundance patterns do not look like that of the average inner stellar halo of the Milky Way (e.g., Shetrone et al. 2001). Interestingly, studies indicate that the outer halo may have abundances similar to those in some dwarfs (e.g., Fulbright 2002). This trend is roughly consistent with our expectations, as demonstrated in Figure 3 and especially in Figure 4.

#### Acknowledgments.

JSB is supported by NASA through Hubble Fellowship grant HF-01146.01-A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. KVJ was supported by NASA LTSA grant NAG5-9064, NSF CAREER award AST-0133617 and a travel award from Wesleyan University. We thank Francisco Prada, David Martínez-Delgado, and the rest of the organizers for an enjoyable and instructive meeting.

## References

Binney, J. & Tremaine, S. 1987, Galactic Dynamics (Princeton: Princeton University Press)

Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2000, ApJ, 528, 145

Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2001, ApJ, 548, 33

Dekel, A. & Woo, J. 2003, MNRAS, 344, 1131

Fulbright, J.P. 2002, AJ, 123, 404

Ghigna, S., Moore, B., Governato, F., Lake, G., Quinn, T., & Stadel, J. 1998, MNRAS, 300, 146

Grebel, E. K. 2001, ApSSS.277–231G

Hashinoto, Y. Funato, Y. & Makino, J. 2003, astro-ph/0208452

Helmi, A., White, S.D.M., & Springel, V. 2003, MNRAS, 339, 834

Hernquist, L. 1990, ApJ, 356, 359

Johnston, K. V., Sackett, P. D., & Bullock, J. S. 2001, ApJ, 557, 137

Klypin, A.A., Kravtsov, A.V., Valenzuela, O., & Prada, F. 1999, ApJ, 552, 82

Lacey, C. & Cole, S. 1993, MNRAS, 262, 627

Majewski, S. R., Ostheimer, J. C., Kunkel, W. E., & Patterson, R. J. 2000, AJ, 120, 2550

Morrison, H. L. et al. 2000, AJ, 119, 2254

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Searl, L. & Zinn, R. 1978, ApJ, 225, 357

Shetrone, M.D.,Côté, P., & Sargent, W.L.W. 2001, ApJ, 548, 592

Somerville, R. S. & Kolatt, T.S. 1999, MNRAS, 305, 1

Wechsler, R. H., Bullock, J. S., Primack, J. R., Kravtsov, A. V., & Dekel, A. 2002, ApJ, 568, 52

Yanny, B. et al. 2000, ApJ, 540, 825

Zentner, A.J. & Bullock, J. S. 2003, ApJ, 598, 49


text:
# Are Supershells Powered by Multiple Supernovae? Modeling the Radio Pulsar Population Produced by OB Associations

## 1 Introduction

Surveys for radio pulsars form the foundation for the considerable contribution which studies of neutron stars (NS) have made to astrophysics. Not only have such surveys identified a variety of unique or unusual pulsars, but the large samples of objects thus accumulated have provided vital information on the velocity distribution, luminosity function, beaming fraction and magnetic field evolution of rotation powered NS (e.g. Ostriker & Gunn 1969; Narayan & Ostriker 1990; Lorimer et al. 1993).

Surveys for pulsars have usually either targeted specific regions such as supernova remnants, globular clusters, unidentified gamma-ray sources or steep-spectrum radio sources (e.g. Biggs & Lyne 1996; Lorimer, Lyne & Camilo 1998), or have carried out unbiased, but usually shallower, surveys over large regions of the sky (e.g. Manchester et al. 1996, 2001). The former set of surveys are aimed at discovering particulary interesting pulsars but cover only a small fraction of the sky; the latter usually find large number of pulsars, and allow statistical studies of pulsar birth parameters to be performed. However, these studies generally require apriori assumptions regarding the pulsar birth rate and the initial spatial distribution of pulsars in the Galaxy.

In this paper, we perform the first study (to the best of our knowledge) of pulsars correlated to supershells in our Galaxy. Expanding giant H i supershells (see e.g. Tenorio-Tagle & Bodenheimer 1988 for a review) have been observed for several decades in 21 cm surveys of spiral galaxies. These nearly spherical structures have very low density in their interiors and high H i density at their boundaries, and expand at velocities of several tens of $`\mathrm{km}\mathrm{s}^1`$. Their radii are much larger than those of ordinary supernova remnants, often exceeding $`1`$ kpc; their ages are typically in the range of $`10^6`$ to a few $`\times 10^7`$ years. The Milky Way contains several tens of them (Heiles 1979, 1984; McClure-Griffiths et al. 2002), and in one case the estimated kinetic energy of expansion is as high as $`10^{54}`$ ergs (Heiles 1979). Similar supershells are also observed in other nearby galaxies (Rand & van der Hulst 1993; Kim et al. 1999).

Whereas it is clear that these H i supershells result from deposition of an enormous amount of energy into the interstellar medium, the energy source is still a subject of debate. Collisions with high-velocity clouds (Tenorio-Tagle 1981) could account for those cases where only one hemisphere is present, and where the required input energy is not too large. However, it is unclear how such collisions could produce the near-complete ringlike appearance observed in some cases (Rand & van der Hulst 1993). Loeb & Perna (1998) and Efremov et al. (1998) suggested that Gamma-Ray Bursts might be responsible for powering the largest supershells. Whereas it is likely that a fraction of the smaller Galactic shells is indeed powered by GRBs given estimates of their rates (e.g. Schmidt 1999; Perna, Sari & Frail 2003), however GRBs might not be able to power the largest supershells if they are beamed (and therefore their energies reduced compared to their isotropic values $`10^{53}10^{54}`$ ergs).

The most discussed interpretation in the literature is that supershells are a consequence of the collective action of stellar winds and supernova (SN) explosions originating from OB star associations (e.g. McCray & Kafatos 1987; Shull & Saken 1995). Several hundreds of stars exploding within a relatively short time are necessary to power the most energetic supershells.

The motivation for simulating and studying the population of pulsars associated with supershells is therefore twofold: first, it would provide an independent test of the multiple SN origin for them; second, if the pulsar population associated with a supershell is indeed identified, it would possibly provide an independent way to constrain pulsar birth parameters. This is because both the birth rate and the birth location of the pulsars would be independently constrained rather than assumed a priori.

The paper is organized as follows: §2.1 discusses the model for the superbubble growth and the corresponding pulsar birth rates; the distributions of the pulsar characteristics and the Montecarlo simulation of the observable population are discussed in §2.2. In §3, the results of our simulations are compared to the data for the supershells which lie in regions of sky covered by past, sensitive pulsar surveys, and predictions are made for some of the other supershells for which current data are not sufficient to make a comparison at this time. Our work is finally summarized in §4.

## 2 Simulating a population of pulsars from an OB association

### 2.1 OB associations and supershell growth

The evolution of a superbubble powered by a continuous energy source444 The approximation made here is that each SN, rather than outputing its energy suddenly at time $`t`$, spreads it output over an interval $`\mathrm{\Delta }t`$ until the next explosion. Mac Low & McCray (1988) showed that after five to ten SNe have exploded, the continuous energy injection model provides a very good approximation to the discrete one. has been calculated by Weaver et al. (1977; see also MacLow & McCray 1988; Shull & Saken 1995). Weaver et al. derived a similarity solution in terms of the equivalent mechanical luminosity of supernovae, $`L_{\mathrm{SN}}`$, ambient density $`\rho _0`$ and time $`t`$. This yields the time-evolution of the superbubble radius

$$R=\left(\frac{125}{154\pi }\right)^{1/5}L_{\mathrm{SN}}^{1/5}\rho _0^{1/5}t^{3/5}(267\mathrm{p}\mathrm{c})\left(\frac{L_{38}t_7^3}{n_0}\right)^{1/5},$$
(1)
and velocity

$$V(15.7\mathrm{kms}^1)L_{38}^{1/5}n_0^{1/5}t_7^{2/5}.$$
(2)
Here $`n_0`$ is the atomic number density (in cm-3), $`t_7=t/(10^7\mathrm{yr})`$ and $`L_{38}=L_{\mathrm{SN}}/(10^{38}\mathrm{ergs}\mathrm{s}^1)`$ is the equivalent of one SN of energy $`E_{\mathrm{SN}}=10^{51}`$ ergs occurring every $`\mathrm{\Delta }t_{\mathrm{SN}}=E_{\mathrm{SN}}/L_{38}=3.2\times 10^5`$ yr.555Note that the numerical value of the interval $`\mathrm{\Delta }t_{\mathrm{SN}}`$ has no special meaning but it is only resulting from the definition of $`L_{\mathrm{SN}}`$ in units of $`L_{38}`$. For a given superbubble, from the observed values of $`R`$, $`v`$ and $`n_0`$, we derive the corresponding age $`t_7`$ and the equivalent mechanical luminosity $`L_{\mathrm{SN}}`$. The number $`N_{\mathrm{SN}}`$ of stars that become supernovae (SN) over the time scale $`t`$ is given by $`N_{\mathrm{SN}}=L_{\mathrm{SN}}t/E_{\mathrm{SN}}`$. Although stellar winds from the OB association will produce a hot bubble before the first SN goes off, these winds are not important compared to SNe for the later dynamics of the supershell (McCray & Kafatos 1987; Mac Low & McCray 1988). Also, while the early times of the bubble evolution are sensitive to the details of the rate of energy injection (which is assumed constant in the Weaver et al. model), the late phase of the bubble is not666There might however be discrepancies if a considerable SN activity occurs near the shell rim (Oey 1996). (Shull & Saken 1995). We should point out that the factor $`N_{\mathrm{SN}}`$ estimated using the equations above is typically a factor of a few smaller than the equivalent number estimated using the energy $`E_\mathrm{E}=5.3\times 10^{43}n_0^{1.12}R^{3.12}V^{1.4}`$ ergs which would be required to produce the shell by a sudden explosion (Chevalier 1974). The latter, larger, number has been the one often used in the literature to estimate superbubble energies (Heiles 1979; Rhode et al. 1999; McClure-Griffith et al. 2002).

The initial mass function for massive stars can be written as (Garmany, Conti & Chiosi 1982)

$$f_{\mathrm{IMF}}(M_{})dN_{}/dM_{}M_{}^\beta ,$$
(3)
where $`\beta 2.6`$. We assume $`M_{\mathrm{min}}=3M_{}`$ and $`M_{\mathrm{max}}=140M_{}`$ for the minimum and maximum mass of the distribution, respectively.

The main-sequence lifetimes of massive stars are given approximately by (Stothers 1972; Chiosi, Nasi & Sreenivasan 1978)

$$t_{}\{\begin{array}{cc}3\times 10^7(M_{}/10M_{})^{1.6}\mathrm{yr}\hfill & \text{if }\mathrm{\hspace{0.33em}7}M_{}30M_{}\hfill \\ 9\times 10^6(M_{}/10M_{})^{0.5}\mathrm{yr}\hfill & \text{if }M_{}30M_{}\hfill \end{array}.$$
(4)
The least massive star that is expected to terminate as a Type II SN has initial mass $`M_{\mathrm{min},\mathrm{SN}}=7M_{}`$ (Trimble 1982). We take $`M_{\mathrm{max}}140M_{}`$ as the mass of the most massive star of the association. Our results are rather insensitive to the precise value of the minimum star mass, $`M_{\mathrm{min}}`$, in the distribution, as stars with $`M<M_{\mathrm{min},\mathrm{SN}}`$ would not contribute to the energetics of the shell. Similarly, our results are also rather insensitive to the precise value of the upper mass cutoff due to both the steep decline in $`dN_{}/dM_{}`$ at large $`M_{}`$, as well as to the fact that low metallicity stars with mass in the range $`40140M_{}`$ collapse directly to a black hole (BH) without any significant explosion energy777A fraction of these stars with a sufficient amount of angular momentum might explode as hypernovae and give rise to gamma-ray bursts. In those particular, but rare cases, an amount of energy would be released into the medium, but probably not significantly larger than that of a standard SN. (Heger et al. 2003). The maximum mass of the star that will leave behind a NS is rather uncertain. However, estimates of rates of heavy element production require stars as massive as 25 $`M_{}`$ to form supernovae (e.g., Hillebrandt 1982). Numerical simulations (e.g., Heger et al. 2003) find that for low-metallicity, non-rotating isolated stars, $`M=25M_{}`$ is the maximum initial star mass that would leave behind a NS star remnant, while stars in the range $`2540M_{}`$ form a BH by fallback while releasing a small amount of energy into the medium. Following these results, here we assume that only stars in the range $`725M_{}`$ leave behind a NS, while those with $`M>M_{\mathrm{SN},\mathrm{max}}=25M_{}`$ (but $`M<40M_{}`$) contribute to power the supershell, but do not leave behind a pulsar. For stars with mass below 7 $`M_{}`$ we assume that their energy contribution is negligible compared to that of the SNe. Note that, if stars had solar metallicity, they would be expected to form NS remnants for a much larger range of initial masses (Heger et al. 2003), and therefore there would be many more pulsars than what we predict here888As a matter of fact, there is a large spread in the metallicity of stars in the Galaxy (e.g. McWilliam 1997), from much lower than solar to about solar. The assumption of low metallicity for all the stars of the association is a conservative assumption for our purposes.. Finally, we consider a model of an OB association with coeval star formation; that is all stars are assumed to be formed at once with no age spread. Once again, note that this is a conservative assumption for our predictions, because if stars were formed with an age spread, the pulsars in the bubble would be generally younger than in the coeval case, and hence would be more likely to be detectable (as younger pulsars have typically larger luminosities and beaming fractions). Also, they would have had less time to travel far from their birthplaces, and therefore a larger fraction of them would be closer to the supershell at the present time.

### 2.2 Properties of the pulsar population

Modelling the intrinsic properties of the Galactic population of pulsars, especially with respect to the luminosity function and initial spin period, has been the subject of extensive investigation in the past few decades (e.g. Gunn & Ostriker 1970; Vivekanand & Narayan 1981; Phinney & Blandford 1981; Narayan & Ostriker 1990; Lyne et al. 1985; Narayan 1987; Stollman 1987; Emmering & Chevalier 1989; Lorimer et al. 1993; Johnston 1994; Cordes & Chernoff 1998). Despite all these efforts, no consensus has been reached on what the intrinsic and birth properties of the pulsars are. The most recent, comprehensive analysis, based on large-scale 0.4 GHz pulsar surveys, has been made by Arzoumanian, Cordes & Chernoff (2002; ACC in the following). Here, as a first step in estimating the expected number of pulsars for a supershell of given age and energy, we use their inferred parameters under the assumption that spin-down is only caused by dipole radiation losses. They find that pulsars are likely to be born with velocities distributed according to a two-component gaussian distribution, with $`40\%`$ of them having a characteristic speed (i.e. magnitude of the 3-D velocity) of 90 km s-1 and the remaining of 500 km s-1. The initial birth period distribution (taken as a gaussian in log), has a mean $`\mathrm{log}P_0(s)=2.3`$ with a spread greater than $`0.2`$999Here we adopt the value 0.3, after verifying that variations in the range 0.3-0.6 do not give statistically different results., while the initial magnetic field strength (also a log-gaussian), is found to have a mean $`\mathrm{log}B_0[G]=12.35`$ and a variance of $`0.4`$. These distributions of initial parameters have been derived under the assumption of no magnetic field decay. The spin evolution of the pulsars is then simply given by:

$$P(t)=\left[P_0^2+\left(\frac{16\pi ^2R^6B^2}{3Ic^3}\right)t\right]^{1/2},$$
(5)
where $`I=10^{45}`$ g cm2 is the moment of inertia of the star, and $`R=10`$ km its radius.

For simplicity and because the ACC radio luminosity model is not obviously scalable to frequencies other than 400 MHz, here we adopt the luminosity function proposed by Emmering & Chevalier (1989)

$$L_\nu =\gamma P^{1.61}\dot{P}^{0.5}.$$
(6)
Lorimer et al. (1993) found that the best fit to the data at $`\nu =400`$ MHz was obtained with $`\gamma =3.5`$ mJy kpc2, which is the value we adopt here. The flux density for a pulsar at a distance $`D`$ is related to its luminosity through the relation $`S_\nu =L_\nu D^2`$. We determine the coefficient of the luminosity function at other frequencies assuming a typical spectral index $`\alpha =1.8`$, where $`L_\nu \nu ^\alpha `$ (Maron et al. 2000). For a given $`P`$ and $`\dot{P}`$, we also allow for a spread in luminosity by considering a log-gaussian probability distribution with mean $`L`$ given by Equation (6) and spread $`\sigma _L`$. As pointed out by ACC, this is a way to parameterize the fact that the pulsar luminosity has an angle dependence with respect to the magnetic axis, and the particular value that is observed by a given observer depends on its orientation with respect to that axis. We find that, with a mean given by Equation (6), a spread $`0.8`$ (in the log) appears to give results consistent with the current data.

For the evolution of the beaming fraction $`f`$ with period we adopt the relation derived by Tauris & Manchester (1998)

$$f=0.09\left(\mathrm{log}\frac{P}{10}\right)^2+0.03,$$
(7)
where $`P`$ is in seconds.

Pulsars are radio active until they slow down to a point where the electric field near the polar cap becomes too small (e.g. Chen & Ruderman 1993). The transition for a given $`P`$ and $`\dot{P}`$ is probably not abrupt, and rather than being described by a “death line”, it is probably better described by a “death band”. Following ACC, we define the probability that objects with given $`\dot{P}/P^3`$ are radio-active as

$$\frac{1}{2}\left\{\mathrm{tanh}\left[\left(\dot{P}_{15}^3/P^310^{\mathrm{DL}}\right)/\sigma _{\mathrm{DL}}\right]+1\right\},$$
(8)
where $`\dot{P}_{15}\dot{P}/(10^{15}\mathrm{s}\mathrm{s}^1)`$. DL$`=0.5`$ and $`\sigma _{\mathrm{DL}}=1.4`$ represent the position and lateral extent of the death band, respectively.

For the pulsar velocity distribution at birth, we adopt the double-gaussian model of ACC described above. The orientation of the pulsar kick at birth is assumed to be random, and the Galactic rotational velocity (220 km s-1) is added vectorially to the initial kick velocity.

Finally, we compute the detector sensitivity using the standard detector equation (e.g. Dewey et al. 1984):

$$S_{\mathrm{lim}}=\frac{\sigma \beta (T_{\mathrm{sys}}+T_{\mathrm{sky}})}{G\sqrt{BN_p\tau _{\mathrm{obs}}}}\left(\frac{W}{PW}\right)^{1/2},$$
(9)
In the above equation, $`\sigma `$ is a loss factor, $`\beta `$ is the signal-to-noise ratio for a detection (assumed to be 8), $`T_{\mathrm{sys}}`$ is the system temperature, $`G`$ the telescope gain, $`B`$ the receiver bandwidth, $`N_p`$ the number of polarizations and $`\tau _{\mathrm{obs}}`$ the observation time. The specific values of these parameters for the various surveys we consider can be found in the papers referenced in §3 below. The observed pulse width is given by $`W=(W_{50}^2+\tau _{\mathrm{samp}}^2+\tau _{\mathrm{DM}}^2+\tau _{\mathrm{scatt}}^2)^{1/2}`$, $`W_{50}`$ is the intrinsic pulse width, assumed equal to 5% of the period, $`\tau _{\mathrm{samp}}`$ is the sampling interval, $`\tau _{\mathrm{DM}}`$ is the dispersion smearing across one frequency channel, and $`\tau _{\mathrm{scatt}}`$ is the broadening of the pulse due to interstellar scattering. The dispersion measure and scattering measure of each pulsar (needed to calculate the above quantities) are computed using the public code101010http://astrosun.tn.cornell.edu/~cordes/NE2001/ by Cordes & Lazio (2003), while the all-sky data are from the 408 MHz compilation of Haslam et al (1982), scaled to other frequencies assuming a spectral index of –2.6 for the brightness temperature of the diffuse Galactic synchrotron background. We neglect the effects of cut-offs in harmonic summing and filtering of the data (see e.g. Manchester et al. 2001), which has a small effect on the sensitivity for periods which are not very short and for large DMs (F. Crawford 2003, private communication).

## 3 Simulating a survey of pulsars in a supershell

The simulation of the pulsar population produced by the OB association that gave rise to a supershell is of the Montecarlo-type. First, the mass of a star is generated from the distribution in Equation (3); only stars whose lifetime is shorter than the present age of the supershell, $`t_{\mathrm{shell}}`$, are followed. The birthtime $`t_{\mathrm{birth}}`$ of a given pulsar is taken as the time that the progenitor star explodes. For each pulsar, the birth properties ($`B`$ field, initial velocity, initial period) are computed according to the ACC distribution described above. The period of each pulsar is then evolved for a time $`t_{\mathrm{puls}}=t_{\mathrm{shell}}t_{\mathrm{birth}}`$ according to Equation (5); if the pulsar is still above the death band, a luminosity is assigned (based on the probability distribution described in §3). To determine whether the pulsar is beamed towards our line of sight, a random number $`x`$ in the interval {0,1} is generated. If $`xf`$ (defined in Eq:(7), then the pulsar is considered beamed towards us. The birth place of the pulsars is assumed to coincide with the center of the supershell; however, because of the relative motion of the supershell with respect to the sun (due to Galactic differential rotation), the coordinates of the supershell on the sky are a function of time. We account for this effect for a pulsar of age $`t_{\mathrm{puls}}`$ by extrapolating backwards by a time $`t_{\mathrm{puls}}`$ to determine the coordinates of the supershell on the sky (and hence of the pulsar birth place) at that specific time in the past. The final position of the pulsar (i.e., at its current age) is determined by integrating its equations of motion in the Galactic potential (Kuijken & Gilmore 1989) for a time $`t_{\mathrm{puls}}`$. Finally, the flux of each pulsar is compared to the detector threshold after computing the dispersion and scattering measure for its position. The simulation is stopped when the number of supernovae reaches what is needed to power the supershell.

Galactic shells and supershells have been studied for several decades since the pioneering work of Heiles (1979, 1984). Pulsar surveys, on the other hand, have never been specifically directed at sampling supershell regions at high sensitivity. Therefore, for the purpose of our analysis in this paper, we use the results of existing surveys, while emphasizing the importance of future, deeper surveys of the supershell regions.

A Montecarlo realization of the 1.4-GHz pulsar population expected to be associated with the two Galactic supershells GSH 242–03+37 and GSH 088+02–103 (which are among the biggest in the Galaxy) is shown in Figure 1 and 2. In both cases, the left panel shows the expected luminosity versus age distribution of pulsars beamed towards us (independently of whether still inside the shell or not). The oldest pulsars have an age comparable to the age of the supershell. The right panel shows the current position on the sky of the pulsars with respect to that of the supershell. Note that some of the pulsars have moved very far from the supershell and therefore they do not appear in the right panels of the figures (which shows the region in the sky around the supershell). Our simulations show that there is a significant probability of finding an enhancement of pulsars in the innermost regions of the bubble. These pulsars are a combination of the youngest sources in the sample, those on the low-velocity tail of the distribution, and those whose velocity vector happens to be almost parallel to the line of sight. The region where GSH 088+02–103 lies has been searched in the Green Bank survey at 370 MHz (Sayer, Nice & Taylor 1997), while GSH 242–03+37 has been searched at 436 MHz by the Parkes Southern Pulsar survey (Manchester et al. 1996), and also had a few pointings at 1.4 GHz in the Parkes Multibeam Pulsar survey (Manchester et al. 2001; A. J. Faulkner, 2003, private communication). The predicted number of detected pulsars for the survey parameters corresponding to the region around each supershell are listed in Table 1. In Figures 1 and 2, as an example, we show the pulsars which would be visible with survey parameters corresponding to those of the current Parkes Multibeam Pulsar survey (filled symbols), as well as those which could not be detected by Parkes but would be visible to the planned Square Kilometer Array (SKA), which will have a sensitivity of about 1.4 $`\mu `$Jy for one minute of integration time (Kramer 2003).

For a quantitative comparison with current data, we selected the most energetic supershells among those at low Galactic latitude (Heiles 1979; McClure-Griffiths et al. 2002), where most of the pulsar surveys have been focused. Our results are summarized in Table 1. For each supershell, we ran 1000 Montecarlo realizations of the pulsar population associated with it, and determined the mean and variance of the distribution for the number of pulsars that are expected to be found within the supershell region on the sky at the current time. To be more specific, if $`l_0,b_0`$ are the Galactic coordinates of the supershell, and $`\mathrm{\Delta }l`$ and $`\mathrm{\Delta }b`$ its total extent (i.e. diameter) on the sky, then the pulsars associated with the supershells are those whose coordinates $`l,b`$ satisfy the condition $`|ll_0|<\mathrm{\Delta }l/2`$ and $`|bb_0|<\mathrm{\Delta }b/2`$, and are within a distance $`D\pm R`$, where $`D`$ is the distance to the shell and $`R`$ its mean radius. These are only a fraction of the ones produced within the association and potentially observable, as shown in Fig. 1 and 2. For each supershell, we computed the distribution of pulsars whose flux is above the threshold of the deepest survey which has covered that particular region of sky in which the shell lies, as well as the distribution of pulsars which would be detectable with SKA. The pulsars in the current catalogue111111http://www.atnf.csiro.au/research/pulsar/psrcat/ which are candidate associations are also reported. These are selected on the basis of their coordinates $`l,b`$ which have to satisfy the condition $`|ll_0|<\mathrm{\Delta }l/2`$ and $`|bb_0|<\mathrm{\Delta }b/2`$, as well as on their distance, which has to be within $`D\pm R`$. However, both the pulsars and the shell distances are known with some uncertainty. For pulsars, the distance as inferred from dispersion measure has some 25% uncertainty on a statistical average121212In some cases distances estimated from the DM may be off by factors of $`2`$ or greater; however, given that ours is a statistical analysis, we adopt the statistical average uncertainty for the model as a whole. (see e.g. Cordes & Lazio 2003), while the errors on the supershell distances are specifically indicated in the McClure-Griffiths et al. sample, and most of them are in the range 10-20 %. Unfortunately no distance errors are provided in the Heiles (1979) sample, and for those we assume a more conservative 25% uncertainty. We then consider a pulsar a possible association with the supershell if it falls within $`D\pm R`$ after accounting for the distance uncertainties in both the pulsar and the supershell.

As Table 1 shows, in a few cases, due to a combination of the large distances of the supershells and the comparatively poor sensitivity level at which the corresponding region was surveyed, no pulsars are expected to have been detected in association with the supershell, and therefore no meaningful constraints can be derived. For the two supershells which lie in the region covered by the Parkes Multibeam survey (GSH 285–02+86 and GSH 292–01+55; McClure-Griffiths et al. 2002), the predictions of the multiple SN scenario at the detection level of that survey are consistent with the number of candidate associations, although they do not rule out other interpretations, such as GRBs or collisions with high-velocity clouds, as in both these scenarios one would not expect any pulsar enhancement correlated with the supershell. Deeper surveys of those regions are needed to test the multiple SN scenario. However, an interesting constraint can already be derived for GSH 242–03+37 (Heiles 1979). Our results (see Table 1) show that there is a 95% probability that the supershell was not the result of multiple SNe. Deeper surveys are necessary to set much tighter constraints. Our current results, in fact, are rather dependent on the luminosity function (Lorimer et al. 1993) that we have adopted here; however, recent new discoveries of low-luminosity radio pulsars (Camilo et al. 2002a, 2002b) might indicate that young pulsars might be fainter than previously realized, and this could be an alternate explanation for the lack of pulsars in the supershell at the current sensitivity threshold. However, when the sensitivity is sufficiently high that even the faintest pulsars can be detected (as it will be with the planned SKA instrument), then the particular details of the luminosity function will be irrelevant for the proposed experiment. Also note that, while current pulsar surveys are able to probe mostly the pulsars within our Galaxy, the SKA survey will be able to detect a large fraction of pulsars also in the LMC, which has several giant supershells (Kim et al. 1999). The analysis that we are proposing here can therefore be made also for the LMC supershells in the SKA era.

Before concluding we should note that, besides depending on the total energetics131313We recall that the Weaver et al. (1977) solution that we have adopted yields an energy value that is a few times smaller than the Chevalier (1974) solution for sudden injection, often used in the literature for supershells. Adopting the latter model would entail us predicting many more supernova explosions, and hence detected radio pulsars, than calculated above. of the supershells (which determine the total number of SN explosions needed), our results also rely on the pulsar birth parameters, and in particular on their velocity distribution, magnetic field, and periods. Here we have adopted the ACC results for the case with braking index $`n=3`$, which corresponds to pure dipole losses. In reality, braking indices are often found to be different than 3 (e.g. Lyne et al. 1996), possibly due to a variety of causes (e.g. Menou et al. 2001 and references therein), and this affects the pulsar evolution. ACC determined the initial birth parameters also for the case $`n=4.5`$, and found that they were not significantly different. In particular, the velocity distributions for the two cases are consistent within the error bars. Given that both sets of initial parameters are determined so that they yield the same pulsar distribution today, different initial conditions should not affect much our predictions for the observability of pulsars in the bubbles today. At any event, future population studies resulting from the analysis of the Parkes Multibeam Survey will be very important for further constraining the underlying pulsar distributions (Lorimer 2003; Vranesenic et al. 2003) and, as a matter of fact, the study of the pulsar population associated with the largest supershells can be potentially used as an independent way to constrain pulsar birth parameters.

Finally, we note that, if the giant supershells are indeed powered by multiple SN explosions, then besides the several tens of NS visible as radio pulsars, there should be hundreds of other NSs (all the pulsars not beamed towards us and all the ones that have passed the death line), plus a smaller number of BHs, in the region delimited by the superbubble. The only hope to possibly detect such large concentrations of compact objects would be through accretion from the interstellar medium on a Bondi-type mode (see Blaes & Madau 1993 for NSs and Agol & Kamionkowski 2002 for BHs). However, current limits for detection of isolated, accreting NSs, may suggest that accretion occurs at a rate which is likely well below the Bondi value (Perna et al. 2003; Toropina et al. 2003), and therefore this large concentration of isolated NSs and BHs in supershells might have to await more sensitive X-ray instruments to be possibly detected in X-rays.

## 4 Summary

Pulsar searches so far have mostly focused either on specifc regions such as supernova remnants and globular clusters or on a large contiguous region of sky. In this paper we have, for the first time, simulated the pulsar population that should be associated with the giant supershells observed in our Galaxy and in nearby galaxies, if their interpretation as due to multiple SN explosions is correct. We have found that, for the largest supershells, several tens of pulsars should be present within the supershell region.

Failure to detect any pulsar enhancement in the largest superbubbles would put serious constraints on the multiple SN origin for them. Conversely, the discovery of the pulsar population associated with a superbubble would provide new insights into the distribution and birth properties of pulsars.

We thank Andy Faulkner and Froney Crawford for information on the Parkes Multibeam Pulsar Survey, and Mordecai-Mark MacLow for a very useful discussion on estimates of supershell energetics. We also thank the referee, Fernando Camilo, for a careful reading of our manuscript and helpful comments.


text:
# Correlated and Decorrelated Positional and Orientational Order in the Nucleosomal Core Particle Mesophases

## Abstract

We investigate the orientational order of transverse polarization vectors of long columns of nucleosomal core particles and their coupling to positional order in high density mesophases discovered recently. Inhomogeneous polar ordering of these columns precipitates crystalization of the 2D sections with different orientations of the transverse polarization vector on each column in the unit cell. We propose possible scenarios for going from the 2D hexagonal phase into distorted lamellar and related phases observed experimentally.

The problem of efficient genome compaction in the evolution of eucaryotic organisms was solved on the smallest scale via the formation of DNA - histone protein complexes, that make a fundamental unit of chromatin: the nucleosome . A nucleosome is composed of a 147 bp long DNA fragment, wrapped 1.75 times around the histone octamer core composed of four different histone proteins and a DNA linker that connects two consecutive nucleosomes . Digestion of the linker DNA between different nucleosomes along the eucaryotic genome gives rise to the nucleosome core particle (NCP). NCP has an approximately cylindrical shape of height and radius $``$ 55 Å, and a structural charge of $``$ -165. This complex is stable in aqueous solutions from mM to 750 mM monovalent salt ionic strength .

At high enough concentration the solution of NCPs ceases to be isotropic. Competing, as yet poorly understood , microscopic interactions between NCPs in concentrated solutions give rise to a slew of liquid crystalline phases of fascinating complexity . At low salt concentrations and sufficiently high osmotic pressure NCPs aggregate into columns of almost ideal cylindrical shape which contain dozens of individual particles. These columns then aggregate to form a lamellar phase . At higher salt concentrations this phase gives way to a 2D hexagonal columnar liquid crystalline phase leading at still higher salt concentrations to a 3D hexagonal crystal .

The high asymmetry of the NCP particle, with a pronounced dyadic axis marking the entry-exit point of the wrapped nucleosomal DNA, allows one to actually observe the local orientational order of the NCPs, described by its in-plane dyadic axis vector, conventionally called transverse polarization . In the lamellar phase the dyadic axes of individual NCPs are strongly correlated along the column leading to a non-zero average transverse polarization vector $`𝐏`$ of a column. Cryomicroscopy studies  give strong evidence that the dyadic axes, and consequently the vectors $`𝐏`$, of columns in the two layers of a lamella are in an antiparallel direction. In addition, the columns in the two layers are shifted to a small but constant distance in such a way that the antiparallel dyadic axes form a periodic polar vector field (Fig. 1). SAXS refine this picture  by showing that there are strong correlations between the NCP columns within a single lamella but different lamellae appear to be completely decorrelated (Fig. 1). In this sense the lamellar phase is decorrelated at a macroscopic scale and corresponds to simple smectic ordering. 2D hexagonal phase is characterized, in turn, by a positional order in the plane perpendicular to the columns, but dyadic axes of NCPs are orientationally disordered along the columns (Fig. 1). The corresponding average polar vector of the column is zero.

In what follows we will investigate the nature of positional and orientational order in phases that are induced by a condensation of the periodic antiparallel transverse polarization field modes. Basing our assumptions on experiments we postulate that the transverse polarization order in a lamella described with $`𝐏(𝝆)`$, where $`𝝆=(x,y)`$, is periodic and highly anisotropic. We will show that the transition between the 2D colloidal crystal with disordered dyadic axes and the smectic-like lamellar phase with the in-plane periodic polar order can be understood by introducing a common subphase substructure for both phases. The common subphase is a 2D colloidal crystal with lamellar structure and with the in-plane polar order of individual lamellae similar to that in the actual lamellar phase. However, the correlations between the positions of columns and the orientations of their dyadic axes from one lamella to another make this phase crystalline and different with respect to the smectic-like structure. The transformation from the hexagonal colloidal crystal to the smectic lamellar phase can take place as a sequence of two second-order transitions through the intermediate crystal lamellar phase, or as a direct first-order transition. Order parameters of both transformations can be easily expressed in terms of the periodic transverse polarization fields.

If the transverse polarization vector field condensation starts from a 2D hexagonal phase then in general one would have several different choices for the direction of the wave vector $`𝐤`$ with respect to the symmetry planes of the 2D hexagonal phase. In the hexagonal coordinate system the wave vector of the condensation into the crystal lamellar phase has to be described in Kovalev’s notation  by $`𝐤_{12}=\frac{1}{2}𝐛_1`$, where $`𝐛_1`$ is one of the reciprocal basis vectors. The irreducible star of the $`𝐤_{12}`$ vector has three arms given by $`𝐤_1=𝐤_{12}^{(1)}=\frac{1}{2}𝐛_1`$, $`𝐤_2=𝐤_{12}^{(2)}=\frac{1}{2}𝐛_2`$ and $`𝐤_3=𝐤_{12}^{(3)}=\frac{1}{2}\left(𝐛_1+𝐛_2\right)`$. The basis functions of a polar vector field with this periodicity are given by $`\psi _1=\left(2\widehat{e}_x+\widehat{e}_y\right)e^{i𝐤_1𝐫}`$, $`\psi _2=\left(\widehat{e}_x+2\widehat{e}_y\right)e^{i𝐤_2𝐫}`$ and $`\psi _3=\left(\widehat{e}_x+\widehat{e}_y\right)e^{i𝐤_3𝐫}`$. The physical realisation of the order parameter is given by the symmetry coordinates composed of the column’s polar vectors $`𝐏_i`$. Denoting by $`i`$ the index of the column, in such a way that its $`(x,y)`$ coordinates are $`i=1(0,0)`$, $`i=2(1,0)`$, $`i=3(0,1)`$ and $`i=4(1,1)`$ in the 2D hexagonal unit cell and the 2D polar vector components by $`(P_x,P_y)`$, one gets the following components of the order parameter

$`\eta _1`$ $`=`$ $`P_{x1}P_{x2}+P_{x3}P_{x4}+\frac{1}{2}\left(P_{y1}P_{y2}+P_{y3}P_{y4}\right)`$
$`\eta _2`$ $`=`$ $`\frac{1}{2}(P_{x1}+P_{x2}P_{x3}P_{x4})+P_{y1}+P_{y2}P_{y3}P_{y4}`$
$`\eta _3`$ $`=`$ $`\frac{1}{2}\left(P_{x1}+P_{x2}+P_{x3}P_{x4}+P_{y1}P_{y2}P_{y3}+P_{y4}\right).`$ (1)
The Landau free energy of this system is now given via the orientation probability density function $`\rho (x,y,z)=\rho _0+_{i=1}^3\eta _i\psi _i(x,y,z)`$, in the form

$$=f\left(\eta _i\psi _i(x,y,z)\right)𝑑x𝑑y𝑑z=F(\eta _i).$$
(2)
For the sake of clarity, we have omitted in Eq. 2 the terms stemming from the large-scale elastic deformations of this order. Consequently, the free energy of the system can be considered as a simple function of the order parameter components.

Let us also stress an evident duality between the polar vector carried by a NCP column and the polar vector of its displacement in the solvent which has important consequencies for the structure of its ordered phases. Since the displacements $`𝐮_i`$ of NCP columns and the polar vectors $`𝐏_i`$ of the dyadic axes span the same irreducible representation of the symmetry group of the 2D hexagonal phase the components of the order parameter can be expressed in the same form in terms of the columns’ displacements $`(u_{xi},u_{yi})`$. Minimization of the free energy dependent on both $`𝐏_i`$ and $`𝐮_i`$ shows that they are proportional. This fact permits one to obtain also the positional order and thus also the final structures of the phases.

To establish the exaustive list of the phases induced by the considered condensation mechanism we use the fact that $`F(\eta _i)`$ has to be an invariant function of the irreducible representation spanned by the periodic transverse polarization vector field. Since the basis invariants of this represenatation are $`I_1=\eta _1^2+\eta _2^2+\eta _3^2`$, $`I_2=\eta _1^2\eta _2^2+\eta _2^2\eta _3^2+\eta _3^2\eta _1^2`$ and $`I_3=\eta _1^2\eta _2^2\eta _3^2`$, the free energy Eq. 2 has to have the form $`=F(I_1,I_2,I_3)`$. The minimization of this free energy w.r.t $`\eta _i`$ gives the equation of state. There are six different solutions of this equation of state which we list below in terms of $`\eta _i`$, the unit cell multiplication in the ordered state with respect to the 2D hexagonal phase i.e. $`V_{ord}/V_{hex}`$ and their space groups.

Let us show that the first solution, i.e. $`\eta _10`$, $`\eta _{2,3}=0`$ describes the crystal lamellar state. From Eq. 1 it follows that in this case $`𝐏_1=𝐏_3`$ and $`𝐏_2=𝐏_4`$, or in other words, that the orientations of dyadic axes carried by columns 1 and 3 or columns 2 and 4 are parellel and have the same degree of correlation along the column. Thus, due to the duality, the displacements of the columns $`1`$ and $`3`$ as well as $`2`$ and $`4`$ are also identical. This means that the periodic transverse polarization field induces a lamellar structure, where a single lamella is defined by the columns $`1,3,1,3\mathrm{}`$ for the top part and $`2,4,2,4\mathrm{}`$ for the bottom part of a single lamella. The order parameter in this crystalline lamellar phase induced by the periodic transverse polarization vector ordering is given by

$$\eta _1=(2P_{x1}+P_{y1})(2P_{x2}+P_{y2}).$$
(3)
The simplest phenomenological Landau free energy expansion of an isolated transition from a hexagonal 2D crystal $`\eta _i=(0,0,0)`$ to a lamellar 2D crystal with $`\eta _i=(\eta ,0,0)`$, with basis invariants $`I_1=\eta ^2`$, $`I_2=\eta ^4`$ and $`I_3=\eta ^6`$ can be written as

$$=a_1I_1+b_1I_2+\mathrm{}.$$
(4)
This free energy describes a second order transition between a hexagonal and a lamellar crystalline phases, driven trivially by the parameter $`a_1`$, that can depend on various solution conditions, e.g. the salt concentration. The same periodic transverse polarization vector field can induce also more complicated transitions. Assuming now a general form of the order parameter we can start with the following free energy ansatz

$$_1(\eta _i)=a_1I_1+a_2I_1^2+b_1I_2+\mathrm{}..$$
(5)
The phase diagram in this case is presented in Fig. 2. One can observe that in addition to the lamellar phase, stable for $`b_1>0`$, the diagram shows the $`(\eta ,\eta ,0)`$ phase, to be discussed below, stable for $`b_1<0`$.

The growth of the complexity of the phase diagram as more non-linear terms are included into the Landau free energy is clearly apparent. E.g. if one starts from the free energy

$$_2(\eta _i)=_1(\eta _i)+a_3I_1^3+c_1I_3+d_{12}I_1I_2+\mathrm{}.,$$
(6)
this leads to the form of the phase diagram presented on Fig. 2, where an additional phase, this time $`(\eta _1,\eta _2,0)`$, intercalates between the $`(\eta ,\eta ,0)`$ and the $`(\eta ,0,0)`$ phases.

The structure of other phases listed in Table 1, induced by the periodic ordering of the transverse polarization vector field obtained in highly non-linear free energy models can be quite complicated. We have already analysed the lamellar phase, beside this we also have:
1. $`(\eta ,\eta ,\eta )`$. The corresponding vector field is given in this phase by the vectors of the different columns in the form :

$`𝐏_1=(P_x,2P_y)`$ $`𝐏_2=(0,0);`$
$`𝐏_3=(P_x,P_y)`$ $`𝐏_4=(2P_x,P_y)`$ (7)
The vectors of orientation (and, consequently, the displacements) of the columns are expressed in terms of one polar vector $`(P_x,P_y)`$. Note that, because of the specific interaction between the columns in this phase, column 2 is frustrated : its average polar vector is zero. Consequently, its displacement with respect to the position in the 2D hexagonal phase is zero and the dyadic axes of individual NCPs are disordered along this column. All three components of the order parameter are equal to the following combination of $`𝐏_i`$ vectors : $`\eta =P_{x1}+2P_{y1}+P_{x3}P_{y3}2P_{x4}P_{y4}`$. Analyzing the displacement field in this phase, see Fig. 2, it is obvious that the columns of a hexagonal phase are now displaced to make an inverted hexagonal phase, observed also in experiments , whose vertices are given by the columns $`3^{}`$, $`2`$, $`1^{}`$, $`2`$, $`4^{}`$ and $`2`$.
2. $`(\eta ,\eta ,0)`$. The corresponding vector displacement fields in this phase are given by

$`𝐏_1=\frac{3}{2}(P_x,P_y)`$ $`𝐏_2=\frac{1}{2}(P_x,P_y)`$
$`𝐏_3=\frac{1}{2}(P_x,P_y)`$ $`𝐏_4=\frac{3}{2}(P_x,P_y)`$ (8)
Both non-zero components of the order parameter in this phase are equal to $`\eta =\frac{3}{2}(P_{x1}+P_{y1})+\frac{1}{2}(P_{x2}+P_{y2})+\frac{1}{2}(P_{x3}P_{y3})\frac{3}{2}(P_{x4}+P_{y4})`$. The columns of a hexagonal phase undergoing the orientational and displacement fields, see Fig. 2, form now an orthorhombic phase with a large unit cell. It can be understood as a succession of two series of columns undergoing 1D displacement waves : in the first series the wave is transeverse and in the second one it is longitudinal. Note that the point symmetry of this structure shows only slight deviation from the hexagonal symmetry.
3. Along the same line the vector fields and the structures of the low-symmetry phases $`(\eta _1,\eta _2,0)`$, $`(\eta _1,\eta _1,\eta _2)`$ and $`(\eta _1,\eta _2,\eta _3)`$ can be obtained.

We have shown how a 2D ordering of a transverse vector polarization field can induce a transition into a lamellar crystalline phase. The experimentally observed lamellar phase is however not of this type and is basically equivalent to a smectic phase. The difference is subtle and referrs to the fact that in SAXS in addition to the smectic peaks only reflections with two non-zero indices from the NCPs correlated along the column and within a single lamella are present. No reflection with three non-zero indices which can evidence the correlation between different lamellae is observed. The lamellae are thus decorrelated and the 3D long range order of the transverse polarization vector is destroyed. This means that the macroscopic symmetry of this phase is higher than that of the crystal lamellar phase with $`\eta _i=(\eta ,0,0)`$ and equivalent to the symmetry of a usual lamellar smectic phase ($`G=D_\mathrm{}hT_d^z`$). One can speculate that this decorrelation is related to the difference in correlation lengths of microscopic interactions which insure the stability of NCP mesophases on the one hand and their colloidal character on the other hand. Indeed, the displacements of columns in a 2D colloidal crystal are much greater than typical values of atomic shifts in a close-packed hexagonal crystal. Resulting distance $`d_L`$ between the lamellae in our case can become greater than the correlation length of the interaction between the NCPs but still smaller then the correlation length of the smectic density wave $`\xi _{\rho (z)}`$.

Since the symmetry of the crystal lamellar phase is lower than that of the actual lamellar smectic phase, this decorrelation transition must be accompanied by vanishing of a physical quantity responsible for this structural change. In other words the order parameter of the lamellar correlation transition can be defined. Group-subgroup relationship between these two lamellar phases permits one to show that the physical quantity vanishing through decorrelation transition is the macroscopic transverse polarization vector wave with the wave vector parallel to the smectic layer and the polarization directed perpendicular to the single lamella (1D transverse vector wave). Condensation of this wave leads directly from the smectic-like decorrelated structure to a correlated crystal lamellar phase of the $`\eta _i=(\eta ,0,0)`$ type, with periodic distribution of transverse polarization in columns in each layer. The order parameter of the transition is given by the first harmonic of the wave i.e. $`u^z(k)\widehat{e}_zexp(iky)+u^z(k)\widehat{e}_zexp(iky)`$, where $`u^z(k)`$ is the amplitude of this harmonic and $`k`$ is the wave vector perpendicular to the direction of the columns within a lamella, with the magnitude $`k=\frac{2\pi }{a}`$, where $`a`$ is the distance betweem neighboring columns within a single lamella.

These polarization waves exist localy for each lamella, but the distribution of the origins of these waves in each lamella is random for the whole sample and consequently the average macroscopic value of the wave is zero in the decorrelated phase (Fig. 1). Similar scenario but in a different context was first proposed by P. Bak . The order-disorder character of the transition and the wave nature of the order parameter were clarified. It was also shown that the wave vector of the macroscopic order parameter for the transitions of this type is equal to the wave vector of the local waves (wave vector of an atomic chain in  or of the polar vector wave in an individual lamella in our case). The transition from the decorrelated smectic phase into the correlated $`\eta _i=(\eta ,0,0)`$ type phase is described trivially via the Lifshitz free energy  $`=\alpha _1|u^z(k)|^2+\alpha _2|u^z(k)|^4+\mathrm{}`$ and is thus of second order.

In this correlated-decorrelated transition scenario the crystalline lamellar phase is a common subphase of both the hexagonal crystalline phase as well as the smectic lamellar phase. The hexagonal crystal $``$ crystal lamellar phase transition is driven by condensation of the 2D periodic transverse polarization field, whereas the smectic lamellar $``$ crystal lamellar phase transition is driven by the condensation of the 1D polar vector wave and is characterised by the correlation between the transverse polarization vector distributions in each individual lamella.

Let us finally discuss the influence of chirality on the state of NCP columns. Each NCP particle being chiral their aggregation results in a chiral column with the rotational symmetry of a cylinder without neither mirror planes nor inversion, and the periodicity approximately equal to the individual NCP height. According to the Curie principle ordering of the polar vector $`𝐏_i`$ along a chiral column leads to a simultaneous ordering of the axial vector $`𝐭_i`$, since both $`𝐏_i`$ and $`𝐭_i`$ span the same representation of the column symmetry group. In a chiral periodic aggregate the axial vector has an evident physical interpretation. Non-zero average value of $`𝐭_i`$ corresponds to the correlated tilt of individual NCPs with respect to the column axis (Fig. 1). In the 2D hexagonal phase the NCP columns carry zero average polar vector, then in a column the tilt of the particles is also zero : the axes of individual particles (usually called S-axes) are parallel to the column axis. On the contrary, in the lamellar phase the ordering of the polar dyadic axes induces a tilt of the NCPs S-axes. This simple consequence of the NCP chirality can explain peculiarities of the SAXS spectrum observed in .

One of the authors (R.P.) would like to thank F. Livolant and her collaborators for their hospitality while being a Chercheur Associé fellow of C.N.R.S at Orsay in 2002.