text:
# New radiative transfer models for obscuring tori in active galaxies

## 1 Introduction

Unified models for active galaxies postulate the existence of a circumnuclear dust torus. This torus obscures the central regions from view at certain solid angles, at wavelengths where the dust is optically thick. Various authors have studied the behaviour of a dusty torus in the presence of a strong radiation field. Models have been constructed for the absorption of the optical to FUV emission by the dust, and the re-emission of this energy in the infrared (Pier & Krolik, 1992; Granato & Danese, 1994; Efstathiou & Rowan-Robinson, 1995). These models have been used to interpret IRAS observations of active galaxies. Their main result is that the spectrum from a dusty disk is much broader than a single grey body. However, observed SEDs are often still broader than the models, and show different behaviour of the 10 $`\mu `$m silicate feature.

As discussed in previous work (van Bemmel & Barthel, 2001), it is very likely that there are two dust components in active galaxies. One component is a large scale dust distribution, heated by stars. The other is the obscuring torus. In narrow-line radio galaxies and quasars it is difficult to disentangle the two components, since their broad-band spectra are not well known, due to their distance and the sensitivity of previous infrared instruments. Previous model studies of obscuring tori were only concerned with another subset of the active galaxies, namely Seyfert galaxies, which happen to be nearby and therefore bright and well observed infrared sources. Moreover, these studies focused on reproducing the general width of the SED. Furthermore, the emphasis was on fitting a few selected objects, rather than a complete sample. In previous studies, the model geometry has been fixed, and the dependence of the resulting SED on the geometry and other disk parameters is not well understood. Finally, the models assume standard Galactic dust, which might not be appropriate for luminous active galaxies. Several studies (e.g. Laor & Draine, 1993; Maiolino et al., 2001a, b) have shown that the lack of a 10 $`\mu `$m silicate emission feature and shallow 2200 Å absorption can point to dust that is dominated by relatively large grains, while Galactic dust is dominated by small grains.

More recent models treat the torus as a clumpy medium, which produces a broader infrared spectrum, more consistent with the observations (Nenkova et al., 2002). However, these authors do not vary the geometry of the torus, nor the distribution of the mass. Although it is likely that the torus does not have a smooth density distribution, it is not yet understood how other parameters influence the width of the resulting broad-band spectrum. Moreover, they assume that all infrared emission observed in active galaxies arises in the torus.

Here we will re-examine the infrared SEDs of active galaxies from a theoretical point of view. We employ a new and very versatile radiative transfer code radical (Dullemond & Turolla, 2000), and examine a wide parameter space. Specifically, we study the dependence of the emission on observing angle, optical depth, disk geometry and size, dust mass distribution and dust grain size distribution. The initial goal will be to constrain the parameters of the torus, by comparing the colours to observations, and by fitting the infrared broad-band spectrum of Cygnus A. As a first order approximation we will assume the torus has a smooth distribution of dust. Eventually, these new models can be used to disentangle the two dust components in active galaxies, which will allow us to study star-formation in active galaxies, and the parameters of the obscuring torus in more detail than before.

## 2 Model setup

In our disk model a geometry for the dust, an input spectrum and dust opacities need to be specified. The numerical method is described below in Sec. 3. We have constructed a grid of models, built by varying different parameters. This results in a set of over 150 model SEDs, in which all variables are cross-correlated. Below we describe the different parameters of the models (see also Table 1).

### 2.1 Disk geometry and density distribution

The density of the disk in angular direction is a Gaussian, with a width of $`h(r)`$. The mass density for any point in the disk is then given by:

$$\rho (\mathrm{r},\theta )=\frac{\mathrm{\Sigma }(\mathrm{r})}{h(\mathrm{r})\sqrt{\pi }}\mathrm{exp}\left[\left(\frac{z}{h(\mathrm{r})}\right)^2\right]$$

with $`z=r\mathrm{tan}(\theta )`$ and $`\mathrm{\Sigma }(r)`$ the surface density. Due to the setup of the model, the surface where $`\tau _{20}`$ equals unity does not coincide with the scale height of the disk. For higher opacities this surface lies at higher $`z`$. Therefore the actual opening angle of the torus is smaller for higher opacities, although the width of the Gaussian is kept constant.

The geometry of the disk is thus described by four parameters: the scale height $`h(r)`$, the surface density $`\mathrm{\Sigma }(r)`$, the inner radius $`r_0`$ and the outer radius $`r_{\mathrm{out}}`$. Each of these parameters is varied. For the scale height we assume:

$$h(r)=h_0\left(\frac{r}{r_0}\right)^\alpha $$

and the surface density is defined as:

$$\mathrm{\Sigma }(r)=\mathrm{\Sigma }_0\left(\frac{r}{r_0}\right)^\beta $$

with $`r_0`$ the inner radius of the disk.

For $`h_0`$ we assume it equals the inner radius, i.e. $`h(r_0)=h_0=r_0`$, so the disk is thicker if it is further away from the emission source. This implies that the opening angles are comparable for similar values of $`\tau _{20}`$. The resulting opening angles are consistent with the size of the broad-line region that has to be shielded in a more edge-on orientation, and with the observed opening angles of ionization cones in active galaxies.

The exponents $`\alpha `$ and $`\beta `$ are variable model parameters. They determine the geometry and radial density distribution. We consider a pillbox disk ($`\alpha =0`$), a conical disk ($`\alpha =1`$) or a flared disk ($`\alpha =9/7`$), see also Fig. 1. The flaring exponent of 9/7 is chosen following models for young stellar objects (Chiang & Goldreich, 1997), since there are no better constraints for disks around active galactic nuclei. $`\beta `$ determines the surface density distribution, i.e. the mass concentration with radius. We constructed models with a flat mass distribution ($`\beta =0`$), which implies that the mass is concentrated in the outer regions for the conical and flaring disk. In addition, we have models where the mass is concentrated in the inner regions, either similar amounts in logarithmic intervals ($`\beta =2`$), or less centrally concentrated ($`\beta =1`$).

The radial optical depth of the disk at the equator is defined at 20 $`\mu `$m ($`\tau _{20}`$), and four values are used: 5, 10, 20 and 50. A value below 1 at 20 $`\mu `$m is inconsistent with observed anisotropies in the mid-infrared (van Bemmel et al., 1998, 2000; Heckman et al., 1992). The upper value is limited by the numerical possibilities of radical. $`\mathrm{\Sigma }_0`$ is chosen such that the model has the desired radial equatorial optical depth at 20 $`\mu `$m.

The inner radius ($`r_0`$) is either 3 or 30 pc, the outer radius ($`r_{\mathrm{out}}`$) 30 or 300 pc, which results in three different sizes for the disk. These will be referred to as large disk (L $`=3300`$ pc), inner disk (I $`=330`$ pc) and outer disk (O $`=30300`$ pc).

### 2.2 Dust characteristics

Using the extinction efficiency tables from Laor & Draine (1993) we have constructed an input opacity table. The mass absorption coefficient $`\kappa _\lambda `$ is obtained by summing the values of $`\kappa _\lambda (a)`$ for all separate grain sizes and dust species, using a normalized size distribution $`S_i(a)`$ for each dust species. The normalized size distribution is given by:

$$S_i(a)=\frac{a^qm_{d,i}X_i}{_aS_i(a)}$$

where $`a`$ is the grain radius, $`m_{d,i}`$ the mass of a dust grain of dust species $`i`$, $`X_i`$ the fraction of dust species $`i`$ and $`q`$ a dimensionless exponent that can vary between $`0.5`$ and $`3.5`$. In Fig. 3 the effect of varying $`q`$ on the resulting opacity is shown. We use only graphite and silicate dust in a 50/50% ratio ($`X_i=0.5`$). The grain sizes range from 0.001 to 10 $`\mu `$m.

For $`q=3.5`$ the Standard Galactic Dust or MRN dust is obtained, as described by Mathis et al. (1977). We have made models with $`q=3.5`$ and $`q=2.0`$ to test the effect of the grain size distribution on a conical disk ($`\alpha =1`$) with a constant density ($`\beta =0`$). No silicate 10 $`\mu `$m emission feature has yet been detected in any active galaxy, constituting a major problem for radiative transfer models of obscuring tori so far. The lack of 10 $`\mu `$m emission might indicate that dust in active galaxies is dominated by larger grains than Galactic dust (Laor & Draine, 1993). Since other studies point in the same direction (Maiolino et al., 2001b), we have adopted $`q=2.0`$ as standard.

### 2.3 Input spectrum

The input spectrum is defined from 1011 to 1017 Hz i.e. from 3 mm to 3 nm (400 eV), and divided in 112 equally large logarithmic spacings. This number was chosen to optimize both calculation speed and frequency resolution. The upper limit for the frequency lies at the point where the dust opacity drops precipitously and the exact value is of little consequence. The lower limit lies far enough in the millimeter to ensure complete coverage of the infrared output spectrum.

For the shape of the input spectrum a power-law is assumed:

$$L(\nu )=L_0\nu ^j$$

This allows for two more variables: the spectral slope $`j`$ and the total luminosity. The results do not depend strongly on $`j`$ when it is varied between 0 and 4; we will adopt $`j=2`$. The luminosity is scaled to be 1012 L, integrated over the highest energies (100–400 eV), which is a typical intrinsic X-ray luminosity for radio-loud active galaxies. In Fig. 2 we show the spectrum that we use as input for radical.

### 2.4 Grid setup

The radial grid is set up using $`r_0`$ and $`r_{\mathrm{out}}`$. A fixed number of radial points is defined between the inner and outer radius. The grid is refined at the inner edge to prevent flux overshoots in the innermost grid cells. For the small disks, we have used 40 radial gridpoints; for the large disk 60. The angular grid does not need to extend from pole to pole. The dust density in the wings of the Gaussian is negligible at a few scale heights, and the disk can be mirrored in the equator. The grid is therefore defined from the equator to a maximum angle that extends well beyond the height where the opacity reaches unity. The code assumes that the remaining space has the same density as the grid cells on the edge of the defined angular raster. We use 20 gridpoints in angular direction for all models. The resulting grid is depicted in Fig. 1.

## 3 The radiative transfer code radical

In order to compute the spectrum emitted by the irradiated torus, a 2-D radiative transfer calculation is performed. This is done using radical111For more information, see also http://www.mpa-garching.mpg.de/PUBLICATIONS/DATA/radtrans/radical/, which is an all-purpose 2-D radiative transfer code for axisymmetric circumstellar media (Dullemond & Turolla, 2000). In this section we describe how we use this code to solve the transfer problem of a torus, partially obscuring the active nucleus.

The dust in the torus is assumed to be in equilibrium with the radiation field. Each dust grain has a fine-tuned temperature such that it emits as much energy per second as it absorbs. In the case of an optically thin torus this temperature is easily computed, since the flux of the nucleus is known and simply goes as $`1/r^2`$. In such a case, no radiative transfer is needed, since the temperatures of the grains can be computed. A simple volume integral over the emission is then sufficient to obtain the emitted spectrum.

Observations of radio-loud active galaxies suggest that the opacity of the obscuring disk exceeds unity at least up to mid-infrared wavelengths ($``$40 $`\mu `$m). The calculation of the emission changes significantly in optically thick models. The nuclear radiation is now stopped when it is intercepted by the torus at the inner edge, and slowly diffuses into the dust. The inner edge is therefore expected to be very hot, while the interior of the torus is much cooler. The spectral energy distribution of such a torus with a large range in temperatures will be broader than expected from an optically thin torus. Solving the temperature of the dust grains is much more difficult than in the optically thin case. This temperature depends not only on the incident nuclear flux, but also on the infrared radiation of the torus itself. A grain located many optical depths deep into the torus will see little of the direct nuclear radiation, but instead will see the infrared emission of its neighboring grains. The radiative equilibrium of this grain depends therefore on the emission of the grains elsewhere in the torus. Every element of the torus is radiatively connected to every other part.

### 3.1 The transfer algorithm

At present, the best way to solve this problem is to perform a 2-D radiative transfer calculation. The absorption and re-emission of radiative energy will be modeled as it seeps through the torus. The radiative transfer along each ray through the torus is given by:

$$\frac{dI_\nu }{ds}=\rho \kappa _\nu (B_\nu (T)I_\nu )$$
(1)
where $`B_\nu (T)`$ is the Planck function at the local temperature, $`I_\nu `$ the intensity, and $`s`$ the path length along the ray. The solution of the radiative transfer proceeds as follows. We start with an initial guess for the dust temperature distribution. Using a radiative transfer integrator based on the algorithm of Extended Short Characteristics (see Dullemond & Turolla, 2000, for details), the radiation field at every location in the disk is computed as a function of frequency $`\nu `$ and direction $`\mu ,\varphi `$ of the radiation.

From this radiation field the mean intensity $`J_\nu `$ is computed:

$$J_\nu \frac{1}{2}I_\nu 𝑑\mathrm{\Omega }$$
(2)
where $`d\mathrm{\Omega }`$ is the solid angle. The dust temperature at every point can then be found by solving the equation:

$$_0^{\mathrm{}}B_\nu (T)𝑑\nu =_0^{\mathrm{}}J_\nu 𝑑\nu $$
(3)
Once this (new) temperature structure is found, it is compared to the old one. If the relative difference is larger than 0.001, we go back to the integration of the transfer equation (Eq. 1), and keep iterating until the relative successive difference drops below the convergence criterion of 0.001. We tested the convergence by starting the model with a cold (1 K) and a hot (1500 K) disk and comparing the resulting spectra. The results are indistinguishable, which indicates that the models are fully converged.

### 3.2 Computing the spectrum

The spectrum is obtained from the density and temperature distribution using a ray-tracing algorithm. For each discrete frequency an image of the torus is produced. These images are then integrated to obtain the spectrum $`F_\nu `$. To make sure that the pixels of these images accurately sample all relevant scales in the system, we distribute the pixels of the image not in a rectangle (as in normal digital images), but in concentric circles around the central source. Each circle contains 64 pixels, and the circles themselves have radii such than for every radial gridpoint on the computational domain there are two circles in the image. In this way there will be an equal number of pixel-circles per logarithmic interval in radius, and this guarantees that both the very inner parts of the disk (few pc) as well as the far outer parts (tens to hundreds of pc) will be sampled properly.

The spectra can be computed at any requested inclination. We select inclinations ($`\theta `$) of 5, 30, 50, 60, 70 and 90 degrees. For small $`\theta `$ the inclination cannot be smaller than 0.1 radian, due to limitations of the code, so its lowest value is 5.7 degrees. All 2-D/3-D effects are included. For edge-on inclination ($`\theta =90^{}`$) one will see only the outer edge of the disk, while at more face-on inclination (say $`\theta =5^{}`$) one will also see the central source and the warm inner regions of the disk. The angle at which the inner regions are shielded from view depends on the opening angle of the torus and the optical depth in our models.

### 3.3 Flux conservation and resolution tests

In the ideal situation, the torus will emit exactly the same amount of luminosity as it has absorbed from the central source. This requirement of flux conservation is a powerful check on the results of the code, and this check is in fact always performed by radical at the end of each run. Moreover, radical is built in such a way that radiation in optically thin areas of the system does not loose or gain flux through numerical diffusion. Unfortunately, in areas of large optical depth, flux conservation is not intrinsically guaranteed by the code, and occasionally the flux can be slightly over-predicted. In most cases we have energy conservation to within about 10%, but sometimes (larger disk and/or high opacities) this error can go up to 20%. Refinement of the grid and fine-tuning of the boundaries of the computational domain usually help to reduce this error, but there are feasibility limits to the number of grid-points that can be used. In this paper all models stay within this 20% tolerance, and most of the models in fact stay within $``$10%.

The frequency resolution has to be high enough to prevent flux losses or overshoots. In Fig. 4 the effect of a lower frequency resolution is shown: with half the number of frequency points the flux deviates less than a factor of two. A higher frequency resolution does not significantly change the results.

Changing the scale height at the inner radius also affects the output spectrum. As shown in Fig. 4, decreasing $`h_0`$ mainly results in a lower output flux in the far-infrared. Since the disk is thinner, there is less cool dust in the disk. The decrease of flux is proportional to the decrease in $`h_0`$. It only marginally affects the shape of the spectrum, therefore it can be used to tune the results to observations for individual objects. When comparing colours, the assumed value of $`h_0`$ is irrelevant.

### 3.4 Improvement on previous methods

The large parameter space probed by our models is unique and allows for a full comparison of observations with these models. The dependence of the output spectrum on disk parameters has not been studied in this much detail before. For simplicity we have chosen to deal with a smooth dust distribution, before assessing the effects of a clumpy medium. One of the main problems has always been to reproduce the broad infrared SED observed in active galaxies: the effects of the geometry of the torus, its density, size, or opacity are not clear. We compare the models to available broad band data in order to assess the width of the resulting spectrum. Instead of fitting a large sample of objects, we will compare the colours of our models to the observed colours. Only for Cygnus A, a powerful nearby radio galaxy, will we study the actual fit in detail.

A second important issue has been the occurence of 10 $`\mu `$m emission at small inclinations in most literature models, which is not observed. Although our frequency resolution does not allow an in depth study of the behaviour of the 10 $`\mu `$m feature, we will discuss the effects of altering the grain size distribution on our models. The powerful radio-loud active galaxies lack proper observations of the 10 $`\mu `$m region, and since we mainly aim to understand these objects, no direct comparison to infrared data will be presented in this paper. In a consecutive paper, high resolution models will be discussed, allowing for a direct comparison with new observations from the SIRTF satellite for a range of active galaxy types.

## 4 Results

In order to understand the effect of the different variables, we have defined a reference model, to which all other models are compared. The relevant parameters are given in Table 1. The reference model is chosen to be a conical disk with constant surface density, a size from 3 to 30 pc and optical depth of 10. This will allow a direct comparison to literature models, which generally have similar characteristics. However, we will use $`q=2.0`$, which does introduce a different behaviour when considering the 10 $`\mu `$m silicate feature, but which is a more realistic value when considering active galaxies. The optical depth of 10 is chosen to account for observed differences at 40 $`\mu `$m restframe wavelength in quasars and radio galaxies (Heckman et al., 1992; van Bemmel et al., 1998, 2000). The size of the disk corresponds to sizes found in Seyfert galaxies, although it may be significantly larger in powerful radio galaxies.

### 4.1 The reference model

The model spectrum is shown in Fig. 5. The spectrum peaks between 10 and 20 $`\mu `$m, depending on inclination. For the different inclinations there are clear flux differences with inclinations (anisotropies) up to 100 $`\mu `$m. The knee at 110 $`\mu `$m is a dust opacity feature. There is a slight indication of 10 $`\mu `$m silicate absorption at larger angles.

The temperature and density distribution for the reference model are depicted in Fig. 6. The temperature in the disk is lower than in its optically thin surroundings, which is to be expected. In the inner regions the temperature drops steeply with radius, flattening at larger radii. Within the disk, the temperature in angular direction is more or less constant, with a steep increase at the edges to reach the temperature of the surrounding medium.

Below is described what the effect of varying the different parameters has on the model results, which is shown in Figs. 7, 8 and 9. These figures are a concise representation of the effects of changing the parameters. A complete set of hundreds of figures has been used to study the effect of the variables on anisotropy and 10 $`\mu `$m absorption in the models.

The depth of the 10 $`\mu `$m absorption shown in Fig. 9 is determined using a similar method as Pier & Krolik (1992): a power-law is fitted through the 6.5 and 13.9 $`\mu `$m points. However, since the model curves have a grey-body shape, fitting a power-law will mimick silicate emission where none is present. This is circumvented by adding to the power-law fits the residual from the model with 5 inclination, where no silicate emission or absorption occurs, resulting in a continuum that approaches a grey-body shape in the region 5–15 $`\mu `$m. Subsequently, the difference between the continuum and the 9.7 $`\mu `$m gridpoint is measured. The depth $`\tau _{9.7\mu m}`$ is positive for aborption.

* $`\alpha `$: changing the physical shape of the disk to a pillbox ($`\alpha =0`$) or flaring disk ($`\alpha =9/7`$). The flaring disk and conical disk have similar spectra, slightly varying with inclination. The pillbox disk is less luminous at all inclinations, but this is more obvious at smaller angles. It is hotter than the reference model, with the peak at shorter wavelengths for smaller inclinations. The wavelength of the emission peak does not shift significantly with inclination, but there is an indication of increased anisotropy. The flared models have a higher degree of isotropy than the reference model, but show a similar shift in the peak wavelength with inclination. They also show a hint of 10 $`\mu `$m absorption at small inclinations, but this effect depends strongly on the optical depth and the size of the disk.
* $`\beta `$: changing the surface density from a constant to a more centrally concentrated configuration ($`\beta =2`$). For fixed $`\alpha ,`$ the models differ slightly at large inclinations, and marginally at small inclinations. The more centrally concentrated, the more the spectrum moves to shorter wavelengths. At large $`\theta `$, there is a hint of 10 $`\mu `$m absorption, which deepens when $`\beta `$ changes from 0 to –2. The depth does not depend on the value of $`\tau _{20}`$. The 10 $`\mu `$m absorption is absent in the 30–300 pc disk.
* $`\tau _{20}`$: changing the optical depth at 20 $`\mu `$m from 5 to 50. The higher opacity models have broader infrared spectra and shorter peak wavelengths at small angles. For larger inclination the high opacity models are not much broader, but clearly have more far-infrared emission. The peak wavelength hardly changes at larger angles. An important effect is the increase of anisotropy with increasing $`\tau _{20}`$, but the wavelength at which the dust becomes optically thin does not vary significantly. There is a hint of 10 $`\mu `$m absorption in the low opacity models at large inclination.
* size: the outer disk emits at much longer wavelengths than the inner disk, and shows a stronger anisotropy at longer wavelengths. The large disk has a much broader spectrum at small angles, indicating a larger temperature range. At large inclinations the outer and the large disk are almost identical. The large disk has clear 10 $`\mu `$m absorption. A hint of this is seen in the inner and outer disk.
* $`q`$: using Galactic dust instead of a large-grain dominated mixture. For larger inclinations the 10 $`\mu `$m absorption becomes deeper. There is a hint of 10 $`\mu `$m emission, which strongly depends on the optical depth and disk size. For smaller optical depth, the emission is clearly present. For larger optical depth, the absorption is only seen in intermediate inclinations. A direct comparison of $`\tau _{9.7\mu m}`$ shows minor silicate emission at inclinations around 30, but up to a factor of 5 increase of silicate absorption at angles larger than 40. In the outer disk, the absorption is deeper at larger inclinations, but there is a clear emission feature for small angles. The anisotropy is less pronounced than in the reference model, although it is visible up to 100 $`\mu `$m as well. The peak wavelength is unchanged.

It is illuminating to examine the spatial distribution of the emission in the disk. In Fig. 10 we show images of the reference model, at 30 and 200 $`\mu `$m, and with inclination of 30 and 70. The images display the entire angular grid, causing them to have a circular cut-off, independent of the inclination of the disk. At 30 $`\mu `$m the disk is optically thick and only the inner regions radiate. At 70 the inner regions are not seen directly, but some emission is seen from the opposite inner side of the torus. At 200 $`\mu `$m the dust is optically thin and the entire disk will radiate. At 30 the outermost regions are not as bright as the rest, since they are somewhat cooler. At 70 the inner regions are now seen direct, since the dust is transparent. They are the brightest part of the disk, since this is where the column density is highest.

The effect of varying the parameters is shown in Fig. 11. The top row shows again the reference model. The second row shows the pillbox models ($`\alpha =0`$), which are clearly less bright. This difference is most pronounced at 200 $`\mu `$m. They are of course flatter at 70. The third row displays the $`\beta =2`$ models, i.e. models with equal mass per logarithmic interval. Here the inner regions dominate all the images, since they represent the densest parts of the disk. This is especially clear in the 200 $`\mu `$m images, where the disk is optically thin. The final row contains the images for a disk with $`\tau _{20}=50`$. At 30 $`\mu `$m the emission is more centrally concentrated and the contrast with the outer regions is larger than for the reference model. However, for the 200 $`\mu `$m images the opacity is irrelevant, and the images are comparable.

### 4.2 Comparison with models in literature

The resulting spectrum from our reference model can be compared to literature models with comparable disk and dust parameters. The main difference between our model and the literature models is the absence of significant 10 $`\mu `$m emission at small inclinations, and only weak 10 $`\mu `$m absorption at larger inclinations. However, since all literature models use Standard Galactic Dust, a comparison at 10 $`\mu `$m is unrealistic.

A qualitative comparison of our reference model with literature models with comparable parameters reveils that the general shape of our output spectrum is best matched by the Pier & Krolik (1992) models (hereafter PK92). The models by Granato & Danese (1994) (hereafter GD94) also match the shape quite well, but the Efstathiou & Rowan-Robinson (1995) (hereafter ER95) models are very different, with a much flatter SED in the mid-infrared and often a peak in the near-infrared. Our reference model becomes optically thin longward of 100 $`\mu `$m. Only the PK92 models match this, other models all become optically thin between 60 and 80 $`\mu `$m. In all literature models changing the inclination causes a continous decrease in flux at a fixed wavelength. The GD94 models show a very discrete change in flux between 45 and 75 degrees inlination. As a consequence, there is no obvious difference between their models with 0 and 45 degrees inclination, and neither between the models with inclinations of 75 and 90 degrees, except for the very thick and compact geometries.

For a more quatitative comparison we have calculated an anisotropy value, defined by taking the flux ratio for face-on over edge-on oriented disks at 20 $`\mu `$m. For the literature models that have parameters closest to our reference model, we find this ratio to be $`5`$. This is consistent with the anisotropy value in our reference model.

Clumpy torus models, as described by Nenkova et al. (2002) are in general flatter than our reference model, which is especially prominent at small inclinations. Their anisotropy value is also $`5`$, as in our reference model, and they do not become optically thin shortward of 100 $`\mu `$m. Depending on the mean free path of the clumps, the models have a more pronounced peak around 60 $`\mu `$m, and become more luminous.

A more detailed comparison is not feasible, since the models in literature use different setups and calculation methods. All observed differences between our results and literature models can be ascribed to this.

## 5 Understanding the behaviour

The output of a dusty disk depends mainly on two properties: the temperature and the mass of the dust. If the dust is further away from the emission source, it will be cooler, and therefore emit at longer wavelengths. At shorter wavelengths the dust is self-shielding, since the dust is optically thick for its own radiation. The inner regions will become very hot, but these can only get rid of the energy at wavelengths where the dust is optically thin. As a consequence, the inner regions will be the brightest at the shorter wavelengths, and only in the far-infrared will the entire disk start to radiate. This effect is clearly seen in the images (Fig. 10).

In a pillbox model, the $`\tau =1`$ surface is ’squeezed’ more to the equatorial plane and, hence, compared to the other models less radiation of the central source is intercepted. Pillbox models are therefore less luminous. The difference between conical and flaring disks is less pronounced. Nevertheless, for a flaring disk, some dust at large radii will have a direct line of sight to the central source and be warmer. This effect is best observed at large inclinations (see Fig. 8).

When we vary $`\beta `$ from 0 to –2, more mass will be concentrated in the inner regions. This dust will be hotter, and therefore brighter in the near- and mid-infrared. Since the dust is more centrally concentrated for $`\beta =2`$, there is less mass in the outer regions, and therefore less cool dust, as seen in Fig. 11. As a consequence the far-infrared emission will decrease in flux, and the entire spectrum shifts to shorter wavelengths.

Increasing the optical depth increases the amount of absorbed radiation. For higher values of $`\tau _{20}`$ it will be more difficult to transport absorbed energy into the torus. As a consequence the inner parts get hotter, but the outer parts will remain cooler. This results in a wider temperature distribution, which translates into a broader energy distribution for small inclinations (Fig. 7). For large inclinations, the increased amount of cold dust is most important. Since the inner regions are shielded from direct view, the energy distribution shifts to longer wavelengths.

The distance to the central source is critical in the models. At large inclinations, we cannot view the inner edge of the torus direct. The hot dust is shielded from view by the cool outer regions. For this reason, the large disk and outer disk have comparable spectra. Overall the inner disk is much warmer and brighter at the shortest wavelengths. At small inclination, the whole disk is visible and both the inner and large disk, which have hot dust close to the central source, are much brighter at shorter wavelengths. The large disk and the outer disk show very similar behaviour at the longest wavelengths where the dust is optically thin. The large disk produces less near- and mid-infrared emission than the inner disk, even though its inner radius is similar, since the temperature gradient is steeper. This is also related to the model setup: since the disk is larger, a smaller value for $`\mathrm{\Sigma }_0`$ is required to obtain the requested equatorial optical depth, and hence the total opacity is lower.

Assuming $`q=2`$ instead of using Galactic dust mainly affects the 10 $`\mu `$m silicate feature. Although our models clearly lack significant silicate emission, they still display clear absorption, albeit not as deep as the Galactic dust models. The detailed behaviour of the silicate feature will be the purpose of future high frequency resolution models.

## 6 Comparing the models to observations

### 6.1 Colour comparison

Using the known infrared colours of samples of active galaxies, we can attempt to constrain the parameters of the obscuring medium. We will use the IRAS bands (12, 25, 60 and 100 $`\mu `$m) which are also available as ISO bands. We have selected all 3C radio-loud active galaxies from the literature (Golombek et al., 1988; Neugebauer et al., 1986) with a least three solid detections in three of the IRAS bands, and with small redshifts ($`z<0.5`$). It should be mentioned, however, that IRAS detected only a few radio-loud active galaxies, which all lie close to the detection limits of the instrument. In addition we selected a sample of starburst, ultra-luminous and Seyfert galaxies with four IRAS detections from literature (Sekiguchi, 1987; MacKenty, 1989).

For all radio-loud objects we assured that the infrared emission is of thermal nature. From Heckman et al. (1992) we have constructed the colours for samples of quasars and radio galaxies. All data are overplotted on the models, while varying different parameters in Fig. 12, 13 and 14. The $`\alpha `$ symbol on the axis, annotated with two wavelengths, represents the spectral index or colour, defined as $`F\nu ^\alpha `$, while the $`\alpha `$ symbol in the plots, without annotation, refers to the geometry of the disk. We will refer to the colours as hot ($`\alpha _{12}^{25}`$), warm ($`\alpha _{25}^{60}`$), and cold ($`\alpha _{60}^{100}`$). For an increase in spectral index the flux at the shorter wavelength increases with respect to the flux at the longer wavelength. Hence an increase in the value for the colour implies more flux at shorter wavelength.

The objects within a given class cluster together. The quasars have a narrow range of hot and cold colours, both lie between 0 and $`1`$. The warm colour ranges from 0.3 to $`2`$. The narrow-line radio galaxies have cooler hot colours than the quasars, ranging from $`0.5`$ to $`2`$. The warm colours show approximately the same values, from 0 to $`1.5`$, but the cold colour has a large range from 1 to almost $`2`$. For the broad-line radio galaxies the hot colour is comparable to that of the quasars, but the warm colour is clearly warmer, having typical values between 1 and $`0.5`$. The cold colour is only determined for one object, the rest lack 100 $`\mu `$m detections. The starburst and ultra-luminous galaxies lie clearly in a different region than the active galaxies in the hot vs. warm plot, while their cold colours are comparable to the cold colours of the radio-loud active galaxies. The Seyfert 1 colour distribution is indistinguishable from the colours for the radio-loud active galaxies, but the Seyfert 2 galaxies overlap both the active galaxies and the starbursts in their hot and warm colours.

There are no objects which have a positive hot colour in the entire sample, and only a few which have a positive warm colour. Nevertheless, the colours of the three radio-loud classes do not separate out too much. From an infrared point of view, these three classes are a somewhat arbitrary classification of a continuum distribution of types.

#### 6.1.1 Models with $`\beta =0`$ and $`\tau _{20}=10`$

In Fig. 12, we have fixed $`\beta `$ and $`\tau _{20}`$, and varied the disk size (3 different lines in plot), $`\alpha `$ (the different plots) and the orientation (along the lines, edge-on to face-on is from left to right). The overplotted data are for objects which are presumed to represent a continuum of different orientations, from quasars to narrow-line radio galaxies. The filled symbols are for ensembles of quasars and radio galaxies.
Changing $`\alpha `$ from $`\alpha =0`$ (pillbox) to $`\alpha =1`$ (conical) has a profound effect on the colours of the disk. A further increase to $`\alpha =2`$ (flaring disk) has little effect. These effects are in line with the SEDs shown in Figs. 7 and 8.

#### 6.1.2 Models with $`\alpha =1`$ and $`\tau _{20}=10`$

In Fig. 13, the value for $`\alpha `$ is fixed, as is $`\tau _{20}=10`$. Decreasing $`\beta `$ to $`2`$ the models move up, except the outer disk, indicating more positive warm and cold spectral indices. The hot colour only increases with decreasing $`\beta `$ in the large disk.

#### 6.1.3 Models with $`\alpha =1`$ and $`\beta =0`$

So far we have fixed the value of $`\tau _{20}`$, but this can of course also vary. In Fig. 14 we have fixed $`\alpha =1`$ and $`\beta =0`$ and varied $`\tau _{20}`$ in the different plots. Along the lines we still vary the inclination.
Increasing $`\tau _{20}`$ will widen the range of the hot and warm colours of the models with changing the inclination. The warm colour also becomes cooler with increasing opacity. The cold colours have approximately the same range when varying $`\tau _{20}`$, but become cooler for higher values.

#### 6.1.4 Comparison with observations

As judged from the plots, for all models the cold colours appear inconsistent with the observations. This implies that either the observed 60 $`\mu `$m flux or the 100 $`\mu `$m flux, or both are not fitted by the models. The hot colours in the models span the same range as the observed colours for active galaxies, but in the hot vs warm colour plot, the warm colours in the models do not cover the entire range of observed values. This indicates that already at 60 $`\mu `$m the spectrum starts to deviate. Hence, both the 60 and 100 $`\mu `$m emission in the active galaxies are dominated by a different component.

At shorter wavelengths, the difference in hot and warm colours between starbursts and active galaxies indicates that the 12 and 25 $`\mu `$m emission in active galaxies arises from a component unrelated to star-formation, most likely the active nucleus. Longward of 60 $`\mu `$m the colours of active galaxies are indistinguishable from the colours of non-active galaxies, which implies that the emission process is the same in both groups. This favours star-formation as the dominant process of far-infrared emission in active galaxies. A similar conclusion is reached by studies of infrared colours of Seyfert galaxies longward of 100 $`\mu `$m by Spinoglio et al. (2002), who find that this emission is independent of the presence of an active nucleus. Furthermore, Prieto et al. (2001) conclude that coronal line strength is strongly correlated with near- and mid-infrared emission, but not at all with far-infrared, indicative of a different heating mechanism for the far-infrared emission.

In the hot vs warm colour plot the narrow-line and broad-line radio galaxies are well fitted by the models. For the other colour plots only a few narrow-line radio galaxies are still within the range of model colours. The outlying quasar in the hot vs warm colours is 3C 48, known to have a large star-formation rate, and possible all its infrared emission arises in star-forming regions (Canalizo & Stockton, 2000). Its warm colour is comparable to the warm colour of objects with known large star-formation rates, but its hot colour is clearly warmer, indicative of very hot dust. This indicates that a significant contribution from an active nucleus will cause objects to have a more positive warm colour, while their cold colours are indistinguishable from the colours of star-forming objects.

The inner disk does not fit the observed colours well, except in the hot vs warm colour plot. The large and outer disk do better for $`\alpha =1`$ and $`\alpha =2`$ disks and large inclinations, but also fail in the hot vs cold colours. The pillbox disks do not fit any of the observations. For decreasing $`\beta `$ all models move away from the observations. Increasing $`\tau _{20}`$ moves the models closer to the observed colours, but the discrepancy in cold colours is not explained.

We note with interest that the ultra-luminous galaxies are well fitted by the high opacity models in the cold vs. hot and cold vs. warm colour plots. Also Cygnus A lies closer to the models than the remainder of the radio galaxy population. For the ultra-luminous galaxies it is known that star-formation dominates the far-infrared output, and also Cygnus A is has strong evidence for significant star-formation in the host galaxy. Hence, star-formation can mimic the SED of an obscuring torus and its relative contribution must be assessed by other means, e.g. PAH emission, optical absorption lines, to fully disentangle the two infrared components

#### 6.1.5 Conclusions from colour-colour plots

Most objects have negative cold spectral indices which are not explained by the models. This could very well be due to a contribution from other components, such as additional cold dust from star-formation regions. Observational evidence for a second dust component heated by stars in active galaxies is found in several studies (Prieto et al., 2001; Spinoglio et al., 2002; van Bemmel & Barthel, 2001). Therefore, we conclude that the torus emission only dominates the SED to about 25 $`\mu `$m, while the observed 60 and 100 $`\mu `$m flux arise in a second infrared component.

With respect to the unification theory, there is no trend in the data for radio galaxies and Seyfert 2 to be more consistent with large inclination models, or quasars and Seyfert 1 with small inclinations. In the hot vs warm plots all radio galaxies lie closer to the models than the quasars, while the objects with known large star-formation rates occupy a clearly distinct region, having cooler hot and warm colours. The broad-line radio galaxy colours are consistent with the lack of large scale cold dust, as discussed in van Bemmel & Barthel (2001). If the torus is the only infrared component in broad-line radio galaxies, they should also be consistent with the models in their cold colours, but we lack detections at 100 $`\mu `$m to confirm this. As emphasized above, the differences between models and observations are indicative of a second component in the far-infrared.

### 6.2 Comparison with Cygnus A

For individual objects, the model parameters can be constrained by a direct fit to the observed infrared spectrum. One of the best objects to test our models is the narrow-line radio galaxy Cygnus A, since it is the best observed radio-loud active galaxy in the infrared, and has a rather well determined central continuum. IRAS and ISO have observed Cygnus A in the infrared and additional (sub-)millimeter data are available to determine the non-thermal contribution. From optical observations the size of the obscuring torus has been constrained to lie between a few parsecs inner radius, up to maximum of a few hundred parsecs outer radius (Imanishi & Ueno, 2000). Also, we know that the angle of orientation of the radio jet is $`50`$ degrees (Carilli & Barthel, 1996), allowing us to vary $`\theta `$ around this value. For the models we fix $`\alpha =1`$ and $`\beta =0`$, since these are the best values from the colour-colour fitting. We use $`q=2`$ for the dust grain size distribution.

We have constructed models normalized to the observed hard X-ray luminosity and spectral slope, and taken into account size limitations of the torus in Cygnus A. The inner and outer radius are varied with one extra value: $`r_0=3`$, 10 and 30 pc, $`r_{out}=30`$, 100 and 300 pc. Models are constructed for the same four values of $`\tau _{20}`$ and spectra are obtained for varying inclinations.

In Fig. 15 we show several fits to the spectrum of Cygnus A, with different variables. A small scaling factor is needed to fit the models, which is either due to a different intrinsic luminosity, or different value of $`h_0`$, which are both kept fixed.

#### 6.2.1 Tuning the fit: size of the disk

The size of the disk has the largest effect on the width of the SED. Good fits are obtained for all optical depths for a disk with $`r=10100`$ pc. The reference model does not fit the SED, as expected from the colour-colour analysis. It overproduces near-infrared emission, and cannot account for all the far-infrared emission. The large and outer disk can fit the far-infrared output of Cygnus A, but both underpredict the observed mid-infrared flux. The inner radius of the disk can be carefully tuned to the shape and intensity of the mid-infrared observations. For an inner radius of 10 pc the mid-infrared is well fitted by the model.

With an outer radius of 30 pc all four mid-infrared points are well fitted by the different opacity models. The model with an outer radius of 100 pc fits all the data, although the fit in the mid-infrared is somewhat worse. The longest wavelength far-infrared point is only fitted by the 10–100 pc disk with $`\tau _{20}=5`$. Increasing the outer radius to 300 pc improves the fit for other opacities at the longest wavelengths, but then the fit in the mid-infrared deteriorates. The models with an outer radius of 100 pc or more peak at 60 $`\mu `$m, as observed, and the slope of the far-infrared and sub-millimeter tail is reasonably consistent with the observations.

We have no problem reproducing the width of the observed SED with a single large torus model, as did previous models. However, from optical imaging and spectropolarimetry we know that Cygnus A has large scale dust and star-formation (Jackson et al., 1998; Fosbury et al., 1999), which is expected to contribute to the far-infrared emission. From the colour-colour analysis we also conclude that the far-infrared emission beyond 60 $`\mu `$m is generally not produced by a torus. Therefore, the 10–30 pc disk provides the best model for the torus in Cygnus A.

#### 6.2.2 Tuning inclination and opacity

Small variations in inclination have little effect on the resulting SED and models with inclinations between 5 and 60 degrees fit equally well. Models with higher inclinations are ruled out.

The opacity cannot be fully constrained from the models. For the 10–100 pc disk the $`\tau _{20}=10`$ and 20 models fit best, which is also the case in the 10–30 pc disk. However, in the latter the opacity is not constrained by the far-infrared, which leaves only four datapoints to optimize the fit. The effect of varying $`\alpha `$ and $`\beta `$ has been described above, and is similar in the Cygnus A models. Changing them will not significantly improve the fit.

#### 6.2.3 Conclusions from the Cygnus A fit

From the observations of Cygnus A we can constrain the size and inclination of the obscuring torus. The geometry and density distribution are fixed, varying these parameters would not significantly change the output spectrum. Assuming significant star-formation in the host galaxy produces the far-infrared emission, the best fit for the torus is obtained for a disk with size of 10–30 pc, inclination of 50 degrees and optical depth of 10–20. A single large disk fits all infrared points, but does not allow a contribution from the observed large scale dust, and is therefore not fully consistent with optical observations of young stars and large scale dust in Cygnus A.

## 7 Conclusions

Radiative transfer models have been constructed to study the emission of a dusty torus around an active galactic nucleus. We have varied optical depth, disk geometry, disk size, density distribution and dust grain size distribution, to study the effect on the resulting SED. In most models, the dust becomes optically thin around 100 $`\mu `$m, which implies that the emission is isotropic at longer wavelengths. The models are compared to the observed infrared colours of a sample of radio-loud active galaxies, and furthermore fitted to observations of Cygnus A. The problem with previous models was often that the resulting SED was not broad enough. We have identified this as a matter of size and/or opacity: the larger the disk, or the higher the opacity, the broader the temperature distribution, and, as a consequence, the broader the spectrum. In the fit to Cygnus A we have no problem to reproduce the entire width of the observed infrared SED.

From the comparison to observed colours, we conclude that the tori in active galaxies are relatively large. The inner disk model does not fit any of the observations: the inner radius more likely lies around 10 pc. The outer radius is not well constrained. Pillbox models are rejected on basis of the colour-colour analysis, both the conical and flaring disks are possible. The surface density seems constant with radius; for more centrally concentrated models the colours are not consistent with observations. Varying the optical depth can sometimes improve the level of agreement, but its effect is diverse. An optimum value cannot be derived.

The objects that are best fitted by the models are the radio galaxies. In the hot vs warm colours all broad-line radio galaxies fall entirely in the range of the model colours. This is consistent with observational evidence that in other objects a second dust component contributes, which is not related to a central obscuring torus. This second component alters the shape of the spectrum, and hence the colours, which explains the inconsistency of the models with the observed cold colours. We conclude that the torus dominates only at 12 and 25 $`\mu `$m. As a consequence, selecting active galaxies is best done using their 12 and/or 25 $`\mu `$m power. However, at these wavelengths the quasars are expected to be brighter, for which bias one has to correct.

Hes et al. (1995) find a relation between the 60 $`\mu `$m luminosity and radio power of active galaxies: the more powerful active galaxies are more luminous at 60 $`\mu `$m. They interpret this as evidence that the 60 $`\mu `$m emission is dominated by reprocessed radiation from the obscuring torus. However, the above analysis suggests that the 60 $`\mu `$m emission in active galaxies is dominated by a second component, and hence its power is not indicative of the power of the active nucleus. If star-formation dominates the 60 $`\mu `$m emission in active galaxies, the correlation would imply that the more powerful active galaxies have higher star-formation rates. A better understanding of the process dominating the 60 $`\mu `$m emission can be obtained by comparing individual fits to objects with and without star-formation, such as 3C 48 and the broad-line radio galaxies.

For Cygnus A, the inner radius of the disk is constrained to lie around 10 pc. This value is found to be most consistent with the observations in the colour-colour plots. The outer radius can be anywhere within 30–300 pc, but the best fit is obtained with a value of 30–100 pc. All values are within the observational limits set by optical studies. Although a single torus model fits the entire infrared SED, there is strong observational evidence for a second dust component in Cygnus A. In addition, the quality of the fit in the mid-infrared improves somewhat for a smaller disk. An outer radius of 30 pc is preferred to allow for the observed star-formation to contribute to the far-infrared emission in Cygnus A.

The consequences for unification models are ambiguous. On one hand, we can fit the entire SED of Cygnus A with a single torus model, while we know that other processes must contribute to its infrared emission. This is in itself evidence in favour of the presence of an obscuring torus. On the other hand, the cold colours of most objects are inconsistent with the torus models, indicating that a second component must be present, and a good fit to the entire SED is not possible. We know that this second component must also be present in Cygnus A, and therefore, the relative contribution of star-formation must be assessed with other means to thoroughly disentangle the infrared SED of active galaxies. Only with a proper decomposition of the SED can the actual emission from an obscuring torus be studied in detail.

## 8 Future work

Infrared observations provide a powerful tool to study the emission from an obscuring torus in any type of active galaxy, especially at near- and mid-infrared wavelengths, where the emission from the torus is observed directly. The models presented in this paper are a next step in increasing our understanding of such tori. However, improvement of the models can be expected on several fronts.

As mentioned, the 10 $`\mu `$m absorption feature is sensitive to the parameters of the model. Unfortunately, no data presently exist for radio-loud active galaxies, due to their faintness at mid-infrared data and lack of capable instruments. SIRTF is expected to have the sensitivity to observe a large sample of radio-loud active galaxies, over a large range of redshifts. By constructing models with higher frequency resolution around 10 $`\mu `$m, predictions can be made for future observations of 10 $`\mu `$m emission or absorption in radio-loud active galaxies. In addition the spectral index of the mid-infrared spectrum is sensitive to the parameters, especially $`\alpha _{12}^{25}`$, and to some extent $`\alpha _{25}^{60}`$. Since this is the wavelength range where we are reasonably sure we observe emission from the torus, high resolution spectra are needed. The SIRTF satellite will be able to obtain spectra up to 90 $`\mu `$m for active galaxies over a large range of redshifts, which will provide a good test for the models.

The models presented here are not optimal. Improved models could lead to a better understanding of the physics involved in the obscuring torus around an active nucleus. They could be made self-consistent with a scale height equal to the local radiative pressure scale height. The input spectrum could be adapted to better reflect the emission spectrum from an accretion disk, which often shows a strong bump in the soft X-rays. In addition, there is evidence that star-formation can occur in the torus, which will increase the scale height of the disk. We plan to take this into account as well. Finally, we have not treated indirect heating by the hot gas around the accretion disk. This gas absorbs the X-rays and re-emits the energy in the UV. This UV serves as an additional heating source for the dust. However, improvement of the models is only relevant when better observations are available to constrain the models, and when we have a better understanding of the contribution of star-formation to the infrared SED.

###### Acknowledgements.

The authors wish to thank Xander Tielens for initiating the collaboration and for advice and critical comments during the writing of the paper. Thanks to Peter Barthel for advice, comments and support, and thanks to Julian Krolik for critical discussions. Finally, we acknowledge expert reading by the referee, Željko Ivezić, leading to significant improvement of this paper and many new ideas for the future. CPD acknowledges support from the European Commission under TMR grant ERBFMRX-CT98-0195 (‘Accretion onto black holes, compact objects and protostars’).


format: text

text:
# Grain Boundary Scars and Spherical Crystallography

## Abstract

We describe experimental investigations of the structure of two-dimensional spherical crystals. The crystals, formed by beads self-assembled on water droplets in oil, serve as model systems for exploring very general theories about the minimum energy configurations of particles with arbitrary repulsive interactions on curved surfaces. Above a critical system size we find that crystals develop distinctive high-angle grain boundaries, or scars, not found in planar crystals. The number of excess defects in a scar is shown to grow linearly with the dimensionless system size. The observed slope is expected to be universal, independent of the microscopic potential.

Spherical particles on a flat surface pack most efficiently in a simple lattice of triangles, similar to billiard balls at the start of a game. Such six-fold coordinated triangular lattices BM:1977  cannot, however, be wrapped on the curved surface of a sphere; instead, there must be extra defects in coordination number. Soccer balls and C60 fullerenes KHOCS:1985 ; J:2000  provide familiar realizations of this fact – they have 12 pentagonal panels and 20 hexagonal panels. The necessary packing defects can be characterized by their topological or disclination charge, $`q`$, which is the departure of their coordination number $`c`$ from the preferred flat space value of 6 ($`q=6c`$); a classic theorem of Euler Euler:1750 ; HP:1996  shows that the total disclination charge of any triangulation of the sphere must be 12 Sphere:12 . A total disclination charge of 12 can be achieved in many ways, however, which makes the determination of the minimum energy configuration of repulsive particles, essential for crystallography on a sphere, an extremely difficult problem. This was recognized nearly 100 years ago by J.J. Thomson JJT:1904 , who attempted, unsuccessfully, to explain the periodic table in terms of rigid electron shells. Similar problems recur in fields as diverse as multi-electron bubbles in superfluid helium Leiderer:1995 , virus morphology CasparKlug:1962 ; MarzecDay:1993 ; Viper:2001 , protein s-layers Sleytr:2001 ; PMS:1991  and coding theory Sloane:1984 ; ConwaySloane:1998 . Indeed, both the classic Thomson problem, which deals with particles interacting through the Coulomb potential, and its generalization to other interaction potentials remain largely unsolved after almost 100 years Smale:1998 ; Altschuler:1994/7 ; EH:1997 .

The spatial curvature encountered in curved geometries adds a fundamentally new ingredient to crystallography, not found in the study of order in spatially flat systems. To date, however, studies of the Thomson and related problems have been limited to theory and computer simulation. As the number of particles on the sphere grows, isolated charge 1 defects are predicted to induce too much strain; this can be relieved by introducing additional dislocations, consisting of pairs of tightly bound 5-7 defects Nelson:2002  which still satisfy Euler’s theorem since their net disclination charge is zero. Dislocations, which are point-like topological defects in two dimensions, disrupt the translational order of the crystalline phase but are less disruptive of orientational order Nelson:2002 . While they play an essential role in crystallography on a spherical surface, the configuration and orientation of these excess defects remains undetermined, and can only be fully understood through a combination of theory and experiment. Experimental realizations that probe the subtle structures have, however, been sorely lacking.

We present an experimental realization of the generalized Thomson problem that allows us to explore the lowest energy configuration of the dense packing of repulsive particles on a spherical surface. We create two dimensional packings of colloidal particles on the surface of spherical water droplets, and view the structures with optical microscopy. Above a critical system size, the thermally equilibrated colloidal crystals display distinctive high-angle grain boundaries, which we label “scars”. These grain boundaries are found to end within the crystal, which is not observed to occur on flat surfaces because the energy penalty is too high.

Our experimental system is based on the self-assembly of one micron diameter cross-linked polystyrene beads adsorbed on the surface of spherical water droplets (of radius $`R`$), themselves suspended in a density-matched toluene/chlorobenzene mixture DHNMBW:2002 . The particles are imaged with phase contrast using an inverted microscope. The curvature of the spherical water droplet limits the imaged surface area to between $`5\%`$ and $`20\%`$ of the full surface area of the sphere, depending on the size of the droplet. After determining the center of mass of each bead, the lattice geometry is analyzed by Delaunay triangulation algorithms Delaunay:1934  appropriate to spherical surfaces.

We analyze the lattice configurations of a collection of 40 droplets. A typical small spherical droplet with system size, $`\mathrm{R}/\mathrm{a}=4.2`$, where $`a`$ is the mean particle spacing, is shown in Fig. 1A. The associated Delaunay triangulation is shown in Fig. 1B. The only defect is one isolated charge $`+1`$ disclination. Extrapolation to the entire surface of the sphere is statistically consistent with the required $`12`$ total disclinations.

Qualitatively different results are observed for larger droplet sizes as defect configurations with excess dislocations appear. Although some of these excess dislocations are isolated, most occur in the form of distinctive ($`5757\mathrm{}5`$) chains, each of net charge $`+1`$, as shown in Fig. 1D. These chains form high-angle ($`30^{}`$) grain boundaries, or scars, which terminate freely within the crystal. Such a feature is energetically prohibitive in equilibrium crystals in flat space. Thus, although grain boundaries are a common feature of 2D and 3D crystalline materials, arising from a mismatch of crystallographic orientations across a boundary, they usually terminate at the boundary of the sample in flat space because of the excessive strain energy associated with isolated terminal disclinations. Termination within the crystal is a feature unique to curved space. Thus, our results provide important guidance to determine the configuration of excess defects on a sphere.

To ascertain that the colloidal particles are equilibrated, we use particle tracking routines and subsequent automated triangulation to measure the mobility of the observed screening dislocations. The diffusion of thermally excited colloid particles on the surface of the water droplets results in local rearrangements of the crystal structure. Thermal fluctuations create and destroy dislocations once every few seconds, on average, indicating that the defect arrays reach equilibrium much faster than the observation time of 10-60 min. The equilibration time can also be estimated by the time required for a dislocation to diffuse across typical defect structures. The measured diffusion constant allows us to calculate equilibration times that range from a few seconds to hundreds of seconds, depending on the size of the crystal. Since the spherical crystals exist for 10 to 60 minutes, the system has sufficient time to reach equilibrium. Thus, our observations reflect the equilibrium ground state, as opposed to a history-dependent nonequilibrium effect, which is the case for crystals in flat space.

To quantify the behavior of the scars, we determine the number of excess dislocations per chain for each droplet, with the convention that dislocations are counted as being part of the same array if they are within three lattice spacings, and plot the results as a function of $`\mathrm{R}/\mathrm{a}`$ in Fig. 2. Scars only appear for droplets with $`\mathrm{R}/\mathrm{a}5`$. These results provide a critical confirmation of a theoretical prediction that $`\mathrm{R}/\mathrm{a}`$ must exceed a threshold value $`(\mathrm{R}/\mathrm{a})_c5`$, corresponding to $`M360`$ particles, for excess defects to proliferate in the ground state of a spherical crystal BNT:2000 . The precise value of $`(\mathrm{R}/\mathrm{a})_c`$ depends on details of the microscopic potential but its origin is easily understood by considering just one of the 12 charge $`+1`$ disclinations required by the topology of the sphere. In flat space such a topological defect has an associated energy that grows quadratically with the size of the system, since it is created by excising a $`60^{}`$ wedge of material and gluing the boundaries together CL:1995 . The elastic strain energy associated with this defect grows as the area. In the case of the sphere the radius plays the role of the system size. As the radius increases, isolated disclinations become much more energetically costly. This elastic strain energy may be reduced by the formation of linear dislocation arrays, i.e. grain boundaries. The energy needed to create these additional dislocation arrays is proportional to a dislocation core energy $`E_c`$ and scales linearly with the system size. Such screening is inevitable in flat space (the plane) if one forces an extra disclination into the defect-free ground state. Unlike the situation in flat space, grain boundaries on the sphere can freely terminate BNT:2000 ; DM:1997 ; PGM:1999 ; BCNT:2002 ; Toomre , and our experimental results confirm these theoretical expectations.

One systematic approach to determining the ground state of a collection of $`M`$ particles distributed on the sphere and interacting via an arbitrary repulsive potential BNT:2000 ; BCNT:2002  treats the disclination defects themselves as the fundamental degrees of freedom, with the 6-coordinated particles forming a continuum elastic background. The agreement between the predicted and observed values of $`(\mathrm{R}/\mathrm{a})_c`$ supports the validity of this theoretical approach. The original particle pair potential is replaced by a long range defect pair potential given by $`\chi (\beta )R^2(1+_0^{\frac{1\mathrm{cos}\beta }{2}}𝑑z\frac{\mathrm{ln}z}{1z})`$, for a pair of defects separated by an angular distance $`\beta `$. The potential is attractive for opposite charged defects and repulsive for like-charged defects. The underlying microscopic potential enters only in determining the proportionality constant (equivalent to an elastic Young modulus) and $`E_c`$. Many predictions of this model are therefore universal in the sense that they are insensitive to the exact microscopic potential. This enables us to make definite predictions even though the colloidal potential is not precisely known. It also means that our model system serves as a prototype for any analogous system with repulsive interactions and spherical geometry. To further test the validity of this approach, we show a typical ground state for large M in Fig. 3. The system size here is $`\mathrm{R}/\mathrm{a}=12`$, similar to the droplet in Fig. 1D. The results are remarkably similar to the experimentally observed configuration in Fig. 1D; the only difference is a result of thermal fluctuations, which break the two defect scars in the experiment. This agreement between theory and experiment also provides convincing evidence that these scars are essential components of the equilibrium crystal structure on a sphere.

The theory predicts that an isolated charge +1 disclination on a sphere is screened by a string of dislocations of length $`\mathrm{cos}^1(5/6)R0.59R`$ BNT:2000 . We can use the variable linear density of dislocations to compute the total number of excess dislocations $`N`$ in a scar. We find that $`N`$ grows for large $`(\mathrm{R}/\mathrm{a})`$ as $`\frac{\pi }{3}\left[\sqrt{11}5\mathrm{cos}^1(5/6)\right]\frac{\mathrm{R}}{\mathrm{a}}0.41(\frac{\mathrm{R}}{\mathrm{a}})`$, independently of the microscopic potential. This prediction is universal, and is in remarkable agreement with the experiment, as shown by the solid line in Fig. 2.

We expect these scars to be widespread in nature. They should occur, and hence may be exploited, in sufficiently large viral protein capsids, giant spherical fullerenes, spherical bacterial surface layers (s-layers) and the siliceous skeletons of spherical radiolaria (aulosphaera)ORA:1983 , provided that the spherical geometry is not too distorted. Terminating strings of heptagons and pentagons might serve as sites for chemical reactions or even initiation points for bacterial cell division Sleytr:2001  and will surely influence the mechanical properties of spherical crystalline shells.

Acknowledgements

This work was supported by the Department of Energy, the National Science Foundation, the Harvard Materials Research Science and Engineering Center (DMR-0213805), the Fonds der Chemischen Industrie and the Emmy Noether Programme of the DFG.

Correspondence and requests for materials should be addressed to MJB (email: bowick@physics.syr.edu) or ARB (email:abausch@ph.tum.de).


text:
# Temperature dependence of the charge carrier mobility in gated quasi-one-dimensional systems

## I Introduction

The purpose of this paper is to understand how electron-electron interactions affect the ac conductivity in the stochastic regime of transport in low dimensional materials. Transport in one-dimensional or quasi one-dimensional molecular wires formed with self-assembled discotic molecules1a ; 1b ; 1c ; 1d ; 2a ; 2b ; 3 , nanotubes4 ; 5  at high temperatures or with disorder, and MBE gate engineered confined charged channels, constitutes currently a very hot and technically important topic6 . Other important systems include 2- dimensional layers such as those encountered in smectic liquid crystals7 ; 8 , gate charge injected or doped polymers9 , gate addressed FET devices based on organic materials 10 , thin film transistors, and same such systems when they exhibit a metal insulator transition and superconductors.

One-carrier conduction phenomena are now reasonably well understood, and many workers are in the process of designing structures in order to exploit the subtle phenomena associated with one-carrier resonant tunnelling through large molecules 12a ; 12b ; 12c . The molecules have energy levels structure which can be put in and out of resonance to the energy of the incoming waves. This can be done using gates and nano-gates. The changes in these energy levels caused by contamination or interaction or external fields can then be recorded as a new current-voltage response and used in various device applications. At the same time it is also well known that small low-dimensional structures are more susceptible to electron-electron interactions. As long as one is dealing with a few steady state carriers which can escape before others arrive, mean-field approaches are in order. Simple minded mean-field approaches will however for obvious reasons not work in low-dimensional systems in which carriers interfere seriously with each others pathways. Interactions are known, for example, to play a role in carbon nanotubes at low temperatures, and we believe that many ac response phenomena in conjugated polymers and molecular wires have in the past been wrongly attributed to pure disorder effects. In reality one is often dealing with a combination of both electron-electron interactions and disorder and it is not easy to differentiate the two experimentally unless one has the right theoretical tools. Electronic interactions also come into play with high charging, such as occur for example in electrode or gate charged FET organic interfaces10  and thin film transistors. The same applies to highly doped materials. Charging energy is particularly important in high-band gap, low dielectric constant and low-dimensional molecular materials such as chemically or injection doped smectic and discotic liquid crystals1a ; 1b ; 1c ; 1d ; 7 .

Charges can be injected into organic materials chemically or by high fields and by light and one can measure the frequency dependent conductivity. The frequency dependent conductivity is in principle one of the easiest ways of seeing the effect of interactions. The long range coulomb forces produce, naively speaking, disorder-like potentials which confine and scatter the carriers. They generate in principle a frequency dependence on the motion which looks similar to the one produced by true disorder. The problem in practice is, however, that one always has some disorder and some electrode polarisation effects, and then it is difficult to disentangle these processes from each other.

## II The computational model

Using many particle Monte-Carlo simulations, we study the effect of charge-charge interactions on the ac conductivity of finite and electrode addressed one-dimensional chains following the work in reference 13, . A schematic representation of the system is shown in figure 1. We consider that a molecular columnar wire, which does not interact with its neighboring columns, is charged via a gated electrode and is attached to two metal electrodes, which serve the purpose of source and drain. The gate electrode lies at a distance of $`d`$=1 nm from the column, and its distance from the edge electrodes equals $`L_G`$=1.4 nm. We have also considered the more realistic situation that the gate electrode is $`d`$=30 nm away from the column. A number $`N`$ of charges (negative or positive, in our case negative) is placed along this linear gate, which is parallel to the organic column. Charges (which can be either holes or electrons, and in our case we consider them to be holes) are created in the column and move along its length under the influence of an externally applied field, which can be either DC or AC. These holes can enter or leave the column through the metal electrodes, but asymptotically they establish a steady-state distribution in the column.

In practice, we consider a column of 300 sites, with a lattice constant $`a=0.35`$ nm, which corresponds to a $`0.1`$ micron device. The number $`N`$ of charges in the gate is fixed by the assumed value of the gate field near the plate right at the start and remains constant throughout the simulation. We assume that the negative voltage at the gate causes an equal number of charges (holes) to appear in random positions in the column.

The Coulomb energy $`V_n(t)`$ of a hole located on site with index $`n`$ at distance $`r_n=na`$ from the source electrode, is the sum of repulsive forces due to the presence of other holes and the attractive force due to the negatively charged gate:

$$V_n(t)=V_{n,\mathrm{rep}}(t)+V_{n,\mathrm{att}}=\frac{e^2}{4\pi ϵ_0ϵ_r}\underset{m=1}{\overset{L}{}}\frac{H_m(t)}{|r_nr_m|}+V_{n,\mathrm{att}},$$
(1)
where $`H_m(t)`$ is the hole density of site m at time $`t`$, and assumes the values of 1 or 0, depending on whether this site is occupied by a hole at time $`t`$.

For the calculation of the Coulomb attraction energy on the holes we use the jellium model. We consider that the number of charges in the linear gate are $`N`$, the length of the column is $`L`$, the gate electrode starts at a distance $`L_G`$ from the source and drain electrodes and its distance from the column is $`d`$. Under these conditions, the attraction energy for a hole located on site $`r_n`$ is equal to

$$V_{n,\mathrm{att}}=\{\begin{array}{c}\frac{e^2}{4\pi ϵ_0e_r}\frac{N}{L2L_G}\mathrm{ln}\left(\frac{LL_Gr_n+\sqrt{(LL_Gr_n)^2+d^2}}{L_Gr_n+\sqrt{(L_Gr_n)^2+d^2}}\right),r_n<L_G\hfill \\ \frac{e^2}{4\pi ϵ_0e_r}\frac{N}{L2L_G}\mathrm{ln}\left(\frac{\left(r_nL_G+\sqrt{(L_Gr_n)^2+d^2}\right)\left(LL_Gr_n+\sqrt{(LL_Gr_n)^2+d^2}\right)}{d^2}\right),L_Gr_nLL_G\hfill \\ \frac{e^2}{4\pi ϵ_0e_r}\frac{N}{L2L_G}\mathrm{ln}\left(\frac{r_nL_G+\sqrt{(r_nL_G)^2+d^2}}{r_nL+L_G+\sqrt{(r_nL+L_G)^2+d^2}}\right),r_n>LL_G\hfill \end{array}.$$
(2)
During each step, the transition probabilities for all the charges in the column are computed, so that a charge in site $`n`$ has a hopping transition rate $`W_{n,n\pm 1}`$ towards its neighbors

$$W_{n,n\pm 1}=\{\begin{array}{cc}\nu _0\mathrm{exp}\left(\frac{\mathrm{\Delta }E_{n,n\pm 1}}{kT}\right)& \mathrm{for}\mathrm{\Delta }E_{n,n\pm 1}>0\\ \nu _0& \mathrm{for}\mathrm{\Delta }E_{n,n\pm 1}0\end{array},$$
(3)
The energy difference between the two sites is computed by the formula

$$\mathrm{\Delta }E_{n,n\pm 1}=V_{n,n\pm 1}V_neFa,$$
(4)
where the symbol $`V_{n,n\pm 1}`$ represents the Coulombic energy on site $`n\pm 1`$ when we displace the hole located on site $`n`$ to the site $`n\pm 1`$. The jump frequency $`\nu _0`$ is fixed to 1012 Hz. The external electric field $`F`$ can either be constant $`F_0=`$20 kV/cm, or oscillating with a period $`\omega `$ so that $`F(t)=F_0\mathrm{cos}(\omega t)`$. When a neighboring site is occupied the transition rate towards that site is 0, i.e. we use excluded volume rules.

For the first and last site of the column the quantity $`\mathrm{\Delta }E`$ is modified, since now there is only one neighbour to move to, and injection/absorption to the electrodes is possible. There exists an energy barrier between the valence band of the organic material and the Fermi level of the metal electrode, whose height we denote by $`E_b`$. The energy difference for a hole injection from the electrode to the first column site is

$$\mathrm{\Delta }E_{n=1}^{\mathrm{inj}}=E_b+V_1eFa,$$
(5)
and for hole injection from the other electrode to the opposite edge is

$$\mathrm{\Delta }E_{n=L}^{\mathrm{inj}}=E_b+V_N+eFa.$$
(6)
In the same way, the energy difference for absorption of a hole at $`n=1`$ to the adjacent electrode is

$$\mathrm{\Delta }E_{n=1}^{\mathrm{abs}}=E_bV_1+eFa,$$
(7)
while the same quantity for a hole leaving the column at the other side is given by

$$\mathrm{\Delta }E_{n=L}^{\mathrm{abs}}=E_bV_NeFa.$$
(8)
The case of $`E_b=\mathrm{}`$ corresponds to blocking barriers, where the number of holes in the system is fixed for the entire duration of the simulation and no hole can enter or leave the column.

Once the transition probabilities $`W`$ have been computed for all holes in the system, the time $`\delta t`$ required for each hop is drawn by an exponential distribution,

$$\delta t=\frac{\mathrm{ln}R}{W_{n,n1}+W_{n,n+1}},$$
(9)
with $`R`$ being a uniformly distributed random variable, assuming values between 0 and 1. The shortest time $`\delta t_{\mathrm{min}}`$ of this ensemble (including the time needed for injection/absorption) is chosen and this jump takes place. The total time advances by the same amount $`\delta t_{\mathrm{min}}`$, and the direction of the jump is decided with a probability

$$P_{n,n\pm 1}=\frac{W_{n,n\pm 1}}{W_{n,n1}+W_{n,n+1}}.$$
(10)
When a hole moves from a site to another, the work done by the electric field during the time interval $`\delta t_{\mathrm{min}}`$ is calculated by $`\delta A=\pm eaF(t)`$. The total work done over a time interval is added to the total work $`\mathrm{\Delta }A`$ for that interval.

The above procedure is repeated many times with the holes moving around, leaving and entering the column. The time that a given realization of the system lasts depends on the convergence of the results, which in turn depends on the specific parameters of the problem. Thus, typical times simulated are in the range of 0.01-0.1 ms. In order to eliminate any initial distribution bias, we repeat the whole algorithm for a number of different realizations (typically of the order of 10-100), with different random initial configurations. During all these realizations we monitor the average of the holes mobility $`\mu `$ in the system with time, until this average is well converged. In almost all cases, the error bars for the mobility estimation were less than 10%.

The number of holes $`N_h`$ in the column is a number which varies with time and converges asymptotically to an average value, depending on the exact simulation conditions. When we use blocking electrodes, of course, $`N_h`$ is always equal to $`N`$, the number of negative charges in the gate.

The average hole mobility $`\mu `$ is calculated via the work $`\mathrm{\Delta }A`$ done on all holes during an oscillation period $`\mathrm{\Delta }t=2\pi /\omega `$, and is given by13

$$\mu =\frac{2\mathrm{\Delta }A}{eF_0^2\mathrm{\Delta }tN_h}.$$
(11)
## III Results

We have simulated the gated system described above for both DC and AC externally applied field. When we considered an AC field, we also studied the following two cases: a) include a finite barrier height between the metal electrodes and the column, and b) assume an infinite barrier height, which corresponds to blocking electrodes.

### III.1 DC Mobility

In figure 2 we present the average mobility for the gated one-dimensional column with 300 sites in a DC field $`F`$=20 kV/cm, as a function of the number of charges $`N`$ in the gate. The injection energy barrier height is equal to $`E_b`$=0.21 eV. Neglecting all carrier-carrier interactions, including those from the gate charge, the mobility would scale as

$$\mu =\mu _0\mathrm{exp}\left(\frac{E_b}{kT}\right).$$
(12)
Different curves correspond to runs at different temperatures. Figure 2 shows that the temperature dependence of the mobility is strongly dependent on the value of $`N`$, the charge density. The increase in mobility with temperature is due to two factors. a) The carriers find it easier to come into the chain, and as the temperature increases because of the barrier factor (i.e. as in equation 12) b) the charges can also hop more easily inside the chains against the self-field produced by the other carriers, which means that in equation 12 the factor $`\mu _0`$ should be replaced by $`\mu _0(T,N)`$. The calculated low temperature mobility is very small at $`T`$=50 K (not shown). In this limit the charges are practically immobile and cannot overcome the charge-charge interactions which keep them vibrating around their “equilibrium” positions.

An extremely interesting phenomenon is observed when we look at the temperature dependence of the calculated mobility at different densities. Figure 2 shows that there are densities at which the “activation energy” is practically zero. This means that for these particular values of $`N`$, a charge distribution is established in the given field geometry, where the injection barrier is compensated by the internal space charge and the gate field. One can also look at it as the variation of the calculated mobility with the number of charges $`N`$. When the temperature is constant, the average mobility decreases with increasing $`N`$. This is expected because the carriers have less space available to move in. This observation is actually true for room temperatures and above, but figure 2 shows that this rule is clearly violated at lower temperatures, and for high concentrations. For low temperatures, the mobility decreases with charge concentration, but when the concentration exceeds a certain value (around $`N`$=30, corresponding to a concentration of 0.1) the mobility increases with increasing concentration until $`N`$=100. The mobility decreases again, but exhibits another peak at $`N`$=150. For concentrations higher than 0.5 the average mobility drops extremely rapidly, and the charges are practically immobile. In other words, in the mobility versus $`N`$ plot, we can observe three maxima (located roughly at $`N`$=5, $`N`$=100, and $`N`$=150) and three minima (for $`N`$=30, $`N`$=120, and $`N`$=180, which is the highest concentration we used). One can understand this behavior if one notes that the combination of the intrinsic (work function) barrier, the charging energy coming from other carriers, and the charging energy from the gate within a particular gate geometry, together produce a density dependent pinning barrier which modulates the absorption and injection efficiency. One can see that when absorption/injection efficiency is high, the mobility can be high, and this is caused by the forward drift of all carriers. As soon as the density exceeds the critical number $`N_c`$ the injection is reduced and the incoming carrier has to wait until one has exited . To calculate $`N_c`$ we note that the Coulomb barrier at entry on column is defined by

$$E_{\mathrm{CB}}=\underset{n=1}{\overset{N}{}}\frac{e^2}{4\pi ϵϵ_0na(N)}E_{\mathrm{gate}},$$
(13)
where $`a(N)`$ is the average distance between two carriers at total number $`N`$. With $`ϵ3`$, $`E_{\mathrm{CB}}`$ exceeds $`kT`$ when the number of carriers $`N`$ on the chain is larger than $`N_c4050`$ (the value increases slightly with the temperature).

When the Coulomb repulsion energy at the entry site starts to get much bigger than $`kT`$ the carriers block each other along the column. They can only move sequentially, “in a phase” forced on them by the internal fields. The total energy barrier seen by the carrier on the first site in the column is given by

$$E_{\mathrm{tot}}=E_b+\underset{j}{}\frac{e^2}{4\pi ϵϵ_0r_{ij}}E_{\mathrm{gate}},$$
(14)
where $`E_{\mathrm{gate}}`$ is the “attractive” energy of the gate counterions and depends on the distance and gate geometry, and the sum runs over the position of the carriers. When we increase the distance to the gate counter charge distribution from 1 nm to a more realistic value for conventional gate devices of $`d`$ = 30 nm (see figure 3) 10 ; 14a ; 14b , the mobility versus density behaves quite differently as shown in figure 3. The internal charge density as a function of gate voltage is also shown in figure 3. Note that we now have fewer charges for the same gate charge because the coulomb repulsion is relatively stronger and prevents charges from coming into the column. Note also that the mobility versus $`N`$ oscillations observed in figure 2 no longer appear in the same place.

For higher temperatures, the picture is much easier to understand because the charges have enough energy to overcome the local attractive energy barriers and there is no pinning effect. When $`T`$=450 K and $`N`$=5, for example, it is much easier for holes to enter the system, so that there are always 8 to 9 charges in the column which undergo a forward drift motion, as opposed to 1 or 2 charges in the column for $`T`$=100 K.

For the injection barrier $`E_b`$ value used, the average number of holes in the column during the simulation is close to the value expected from the gate charge, except of course when we use $`d`$=30 nm. Here, the Coulomb repulsion with weak screening keeps the charges out of the column (see figure 3). For $`d`$=1 nm, the gate to charge density matching is almost one to one for all temperatures, and at the lower end of concentrations. At higher concentrations, we observe a reduction in the average number of holes in the column of the order of 10-15%, e.g. for $`N`$=100 we have around 90 holes. This is not surprising and is an effect of the “not metallically screened coulomb repulsion” in the chain 14a ; 14b .

### III.2 AC Mobility

#### III.2.1 Finite electrode barrier height

For the relatively low frequencies of 106 Hz (figure 4a), and for intermediate densities, the mobility increases almost one order of magnitude as we increase the temperature (the injection barrier is limiting the mobility). But, at the low concentration end, when carrier interactions are small the temperature dependence is weak. At very low frequencies, one should of course expect similar results to the ones observed in the limit of a DC field. So, the temperature independence at low concentrations implies that at 106 Hz one is already looking at bulk motion. The barrier no longer affects the “single-particle” mobility. It is important to note that apparently all the temperature dependence in figure 4a is already due to many-body effects. When the frequency of the electric field is increased further, as shown in figure 4b, to the frequency value of 109 Hz for example, the change with temperature is now weaker even for higher densities. On this time scale, the interactions are getting less important as one might expect. Indeed, for very high frequencies, e.g. 1010 Hz as shown in figure 4c, the mobility actually decreases with increasing temperature, and this behavior is monotonic with $`N`$. The change in mobility is however less than half an order of magnitude. This behavior is unusual for hopping systems15 . We believe that this is because the effective disorder is increasing with temperature. The disorder increases because the number of carriers in the column is increasing and this reduces the free volume, as shown in figure 2.

For a given temperature, the calculated frequency dependence is shown in figure 5. Here we have plotted 3 curves: a) $`T`$=193 K, b) $`T`$=343 K, and c) $`T`$=493 K. Consider first the low temperature curve. Here the average mobility at small concentrations is practically constant for all frequencies. This means that the injection barrier has very little effect on the frequency behavior at this temperature. Let us try to understand why. In the spirit of a simple resistance sum16  for a one particle random walker and rigorous in the low frequency limit, we have for a column of 300 sites a transport admittance $`Y`$:

$$Y(\omega )=\frac{e^2a^2}{kT}\left(\frac{1}{i\omega +W_{\mathrm{inj}}}+\frac{299}{i\omega +W_{\mathrm{bulk}}}\right)^1.$$
(15)
The injection rate $`W_{\mathrm{inj}}`$ is activated as in equation 12, and the bulk rate $`W_{\mathrm{bulk}}`$ is a constant for a single particle. From equation 15 one would expect the frequency dependence to come in when the frequency is of the order of the injection rate, that is roughly 108 Hz, and to saturate at the bulk value of 1012. This, however, does not happen in figure 5. At $`N=3`$ the mobility stays roughly at saturation value all the way. The reason for this behavior is as shown in figure 2. The particles in the simulation are moving in a potential made up of the gate field, the electron-electron interaction and the barrier, i.e. as given by equation 14. The barrier at a given concentration may indeed be much lower than the injection barrier, and that is why there is no, or only a weak frequency dependence for this particular concentration. Looking at figure 5 we see that at higher concentrations the mobility becomes strongly frequency dependent (see the corresponding dc limit in figure 2), and rapidly increases with frequency at around 109 Hz. This is now an interaction effect as in equation 14, and it looks as if in the spirit of one body theory, interactions are causing an effective disorder. The interactions appear to have little effect on time scales shorter than $`10^9`$ s. As we raise the temperature further to 343 K, the frequency dependence is now weak for all $`N`$. Above 343 K, there is in effect no strong inhomogeneity in the hop time distribution. This means that the variation in hopping activation energies is comparable to $`kT`$ . Finally, in the extreme high $`T`$ limit, when $`T`$=493 K the mobility actually decreases with increasing frequency. This is most unusual for hopping transport and is in this model due to the extra charging that high temperature causes, as shown in figure 2. At high densities, a real time analysis shows that the charges move like a collective “line of billiard balls”. A carrier can only come in when one has escaped, and the mobility decreases with $`N`$ at high temperatures. In figure 6, we present the steady state number of holes in the column as a function of temperature for different gate charges. This number is an average over all different frequencies used. It is evident that there is a considerable net charging effect in the system. The charging is more prominent for smaller $`N`$. For example, when $`N`$=3 and $`T`$=593 K the average number of holes in the system is around 7, while for $`T`$=193 K it is 2 (the absolute value of these numbers depends of course on the value of the interface energy barrier $`E_b`$, but the basic trends remain the same). This charging effect is responsible for giving the lower value of the mobility observed at higher temperatures (see for example figure 2). Normally, in hopping transport we expect the mobility to go up with temperature. Here however, the steady state density of holes, which also influences the magnitude of the hole-hole scattering, is temperature-dependent. There are competing influences: whereas temperature makes hopping easier, the number of holes also goes up, and this increases therefore the hole-hole interactions, or what is also the degree of disorder in the system.

#### III.2.2 Blocking electrodes

We also studied the case where there is an infinitely high energy barrier among the electrodes and the molecular wire. This corresponds to blocking electrodes, which means that the system is “closed” in the sense that no carrier can (subsequent to charging) enter or leave the column. In figure 7 we can see that at high frequencies the average mobility is constant and almost the same for different concentrations. This is because in this limit we are looking at bulk hopping. As we lower the frequency, the mobility drops steadily to very small values. The carrier displacement is at low frequencies limited by the boundary and the decrease in mobility is expected in a confined system. The effect is also concentration dependent. Going up in frequency, we observe that the higher the hole densities, the later the saturation frequency. For the same frequency, the difference in the mobility for the low and high concentrations can be up to 2 orders of magnitude. This is true for concentrations in the range from 0.01 to 0.12 charges per site.

The corresponding temperature/frequency dependence is shown in figure 8. The curves seem to lie practically all on the same points except in the limit of high frequencies. Here, the average mobility is temperature dependent and the interesting point is that the mobility decreases with temperature. This is a many-body effect and says that temperature is enhancing the coulombic disorder by allowing configurations which provide a wider distribution of activation energies.

In order to understand this behavior more accurately, we monitored the motion of the charges in the column as a function of time. The main conclusion is that the charges follow the periodic pattern imposed by the external alternating field. When the field frequency is high, the charges do not have enough time to come closer to each other or to the electrodes, so they fluctuate rapidly around their initial positions following the field, and this gives rise to the given calculated mobility value shown in figures 7 and 8. As the temperature increases, the thermal energy of the charges allows them to move more easily against the direction of the fields. This increases the electron-electron disorder energy, and in this way temperature now reduces the overall mobility! This picture explains the decrease of the mobility with temperature at a given value of (high) frequency, as shown in figure 8.

For a given high frequency range, this picture remains basically the same. Suppose that 3 charges are located in the column at some equilibrium positions. For a frequency 1012 Hz, for example, they will make very short moves towards one or the other side. When the frequency is switched to 1010 Hz, the charges will again follow the field, but despite the larger distances traveled in any one direction, the charge closest to the electrode will still not have enough time to seriously be blocked by the confining electrodes. Similarly, the other charges will not have enough time to feel the changes in each others positions, and they will continue to oscillate, essentially undisturbed, thus giving rise to the same mobility as for higher frequencies. This explains the constant value of the $`\mu (N)`$ mobility versus frequency at high frequencies.

When the frequency is small so that the displacement can reach the electrode boundary, the charges are pushed towards one side of the column. They come very close to each other and thus the mobility decreases with $`N`$ (only the one furthest from the edge can move easily, but when the charge tries to escape, the external field pushes it back). As the frequency gets lower and lower, the effect becomes more pronounced because now the biasing field remains at the same magnitude and in the same direction for a longer time $`1/\omega `$. This is the reason why the mobility falls rapidly and monotonically with lower frequencies. In fact, when using blocking electrodes, the DC mobility is rigorously 0 as it should be, with the (+) charges ‘piling’ up near the (-) electrode.

## IV Conclusions

We have continued the study of charge transport in interacting quasi-one- dimensional gated systems with emphasis on the temperature dependence of the mobility. We have seen that the temperature changes the mobility in the presence of a DC field but that the effect is very much dependent on charge concentration. The temperature dependence is partially due to the assumed injection barrier and partially due to the fact that carriers can “push their way in” like a line of “billiard balls”. Both factors work together to produce the final result shown in figure 2. Temperature can cause the net steady state charge to vary and this then also affects the charge-charge interaction efficiency (figure 2).

The many body interactions do indeed cause a frequency dependence of the mobility. However, using the same parameters as for doped discotic HAT618 , we find that the frequency dependence comes in at too high a frequency. It seems that the experimental data in references 1a, ; 1b, ; 1c, ; 1d, ; 18,  can only be understood if we also include disorder.

When we have blocking electrodes, the interactions influence the density and the temperature dependent ac mobility. The charges increase the confinement and the mobility decreases with concentration and with temperature in the way one would expect for finite ordered segments. Finally, we now need to extend our theoretical tools to handle 2-dimensional systems, so that we can model smectic liquid crystals, mixes of smectics, FET’s, thin-film transistors (TFT), and also the transport through large molecules such as nanotubes and nanotubes filled with fullerenes (pea pods)17 . So far the work has focused on quantum transport. There are and will be, we believe, many new and exciting experimental results coming up in this field in the stochastic and diffusion limit.

###### Acknowledgements.

One of the authors, B. Movaghar, would like to thank NEDO Int. Joint Research Grant for financial support.