text:
# A First Mass Production of Gas Electron Multipliers

## Abstract

We report on the manufacture of a first batch of approximately 2,000 Gas Electron Multipliers (GEMs) using 3M’s fully automated roll to roll flexible circuit production line. This process allows low-cost, reproducible fabrication of a high volume of GEMs of dimensions up to 30$`\times `$30 cm2. First tests indicate that the resulting GEMs have optimal properties as radiation detectors. Production techniques and preliminary measurements of GEM performance are described. This now demonstrated industrial capability should help further establish the prominence of micropattern gas detectors in accelerator based and non-accelerator particle physics, imaging and photodetection.

A number of new radiation detector designs, collectively referred to as Micropattern Gas Detectors (MPGDs)  have recently emerged in response to the extraordinary demands of next-generation High Energy Physics (HEP) experiments, namely the ability to respond to a high counting rate and integrated particle flux, superior radiation resistance and fine spatial resolution. Common to these designs is the presence of a large voltage drop (several hundred volts) across microstructures immersed in a suitable gas mixture. Electrons originating from ionization of the gas in a conversion volume drift to the region of the microstructures where the intense electric field allows gas amplification to occur. Due to the confined amplification regions, slow positive ions are removed immediately from the amplification volume, increasing rate capability by several orders of magnitude compared to wire based gas detectors.

The attractive features of these detectors have lead to a growing number of applications in many fields. For instance, MICROMEGAS chambers  can be found nowadays in medical digital X-ray imaging equipment , where a high sensitivity in low intensity radiation fields results in a diminished dose to the patient, while profiting from an enhanced image contrast. Similarly, photocathode-coated MPGDs promise to surpass photomultiplier tubes in light detection efficiency, reduced cost and speed . Other emerging applications are industrial imaging  and X-ray astronomy . Reviews of these can be found in .

Recently it has been proposed to extend MPGD use to the field of non-accelerator particle physics , where uses would be numerous in view of their simplicity, the possibility to easily construct MPGDs out of radioclean materials and their very low energy threshold. It is in the context of the development of a new type of neutrino detector  that we attempted to manufacture MPGDs in large numbers and with near-perfect reproducibility, using an industrial approach. This effort may nevertheless have repercussions in satisfying the large demand for MPGDs in accelerator based physics. For this first attempt we chose a popular design, the Gas Electron Multiplier (GEM)  given its particular simplicity of design. A GEM consists of a $`50\mu m`$-thick polyimide (KAPTON$`^{\text{TM}}`$) film copper clad on both sides, perforated with a regular matrix of small holes (diameter few tens of $`\mu `$m) produced by photolithography. When a voltage difference is applied between the two sides of the GEM, a large electric field is produced in the holes. Electrons that enter the holes undergo gas amplification. A remarkable advantage of GEMs is the possibility of building multi-stage amplification layers , where electrons are transferred from one GEM to the next, undergoing successive avalanches and yielding very large charge gains. The resulting high-efficiency for single electron detection looks particularly attractive to us, when considering the small energy depositions expected from low-energy neutrino recoils .

We report here on preliminary tests and observations made on a first batch of $`2,000`$ GEMS (Fig. 1) produced using 3M’s high volume, wide web, roll to roll, adhesiveless flexible circuit (FLEX) circuit making process (Fig. 2). At the time of this writing every indication points at their having a satisfactory performance as radiation detectors. Our preliminary results are encouraging but testing is in an early stage. For example we have not yet studied the resistance to radiation (ageing) of the GEM. We however feel that the widespread demand for GEMs by numerous research groups justifies the early release of our findings.

Flexible circuits are utilized in a variety of applications such as inkjet printer cartridges, hard disk drives, liquid crystal display modules, and IC packages among others . These applications have a variety of needs that are met with various FLEX circuit constructions. These can be grouped into two categories: 3-layer and adhesiveless FLEX circuits . A 3-layer FLEX is appropriately named since it is constructed from a copper foil, a polymeric film, and an adhesive to bond foil to film. In many applications such as Hard Disk Drives and GEMs, the presence of the adhesive would create outgassing and ionic problems. An adhesiveless FLEX circuit (also referred to as a 2-layer FLEX) requires alternative means for securely bonding the copper to the polymer. The two primary methods for fabricating an adhesiveless FLEX circuit are a) direct metallization of the polymeric film and b) casting of liquid polyimide onto the Cu foil. After the substrate has been created, the copper and polyimide materials must be patterned to form the desired geometry for the application. The copper pattern can be formed by using either an additive or subtractive circuitization process. The process flows for each are illustrated in Fig. 3.

The additive process consists of applying a photo resist imaged with the desired copper pattern to a $`50\mu m`$-thick polyimide film, directly metallized on both sides. The copper is then electroplated to the appropriate thickness onto the exposed flash layer. This plating technique can allow for a wide range of copper thicknesses ranging from 4 to 36 $`\mu `$m. As shown in Fig. 4, this additive circuitization process can achieve very fine copper features , down to $`20\mu m`$ trace and $`20\mu m`$ space on 1-metal layer FLEX circuits ($`30\mu m`$ traces and $`30\mu m`$ spaces on 2-metal layer FLEX circuits). As can be seen in the trace cross section in the figure, the sidewalls on the additive copper are nearly vertical.

GEM foils were also manufactured using the subtractive process flow outlined in Fig. 3. The subtractive structure is believed to be similar to the construction described by Bouclier et al. . The side walls of the copper openings have a somewhat shallower slope than in the additive process. GEMs produced with these two different circuitization processes showed some performance differences described below.

One lot from each fabrication method has been produced so far, each containing $``$30 identical panels of 33 GEM elements (Fig. 1). Production of a much larger number of panels per lot, up to a few hundred, is possible. In both cases the chosen design was the so-called “standard GEM” , i.e., $`80\mu `$m holes in an hexagonal pattern with $`140\mu `$m pitch and a biconical transversal hole cross section. In other words, the innermost part of the holes exhibits a reduced opening of $`55\mu `$m, a characteristic also found in most GEMs built elsewhere. The use of Dupont E-film KAPTON$`^{\text{TM}}`$ as the substrate does not allow to reduce this opening any further. In present lots the copper thickness was fixed at $`12\mu `$m to insure the success of these first trials. In successive attempts this will be further reduced, a feature of interest for tracking devices where multiple scattering in the detector must be minimized.

The surface quality of both lots has been studied via SEM (Fig. 5) and AFM (Fig. 6). Slight copper opening irregularities are observable mostly on one side of these first additive GEMs. This may lead to gain inhomogeneities across the GEM surface: therefore we have concentrated at first on the characterization of the subtractive lot. The additive fabrication process has proven to be challenging: small polyimide ribs stemming from the interior of the holes were initially observed to envelope the edges of copper openings. Additional treatment of the lot removed these but resulted in a slightly diminished copper to polyimide attachment111It must be noted that the method used to test copper to polyimide attachment is probably too stringent, consisting of firmly attaching adhesive tape to the GEM surface and swiftly peeling it off. Only some additive 3M GEMs are seen not to pass the test. With any luck a GEM should not have to withstand such abuse during normal operation. and the previously mentioned irregularities. The smoother copper surface quality obtained with this method (Fig. 6) is nevertheless a redeeming quality that justifies further exploration: it may be of importance in applications where total inhibition of field effect electron emission is sought, as is the case in  and other efforts concerned with single-electron detection . The subtractive surface quality exhibits apparently innocuous copper microcrystallite growths (Fig. 5): while we have not observed any sparking nor unexpected behavior from their presence, an attempt will be made to remove them in new lots. In order to inhibit their growth it should suffice to reduce the concentration of dissolved copper in etchant baths.

Figs. 7-10 incorporate the extent of our preliminary characterization of subtractive 3M GEMs. Fig. 7 shows typical leakage currents measured in a number of randomly selected GEMs. They consistently display values comparable to previously produced GEMs. However, we have not yet observed any need to “cure” or “burn” 3M GEMs against shorts able to suddenly raise this current into the $`\mu `$A range, as is sometimes necessary with other GEMs. The good behavior of these leakage currents over periods of several hours probably comes from the homogeneous surface treatment that the fully-automated roll to roll process guarantees, together with the use of high-purity polyimide, free of any fillers. Each part of each GEM foil receives an identical treatment in every fabrication step, something hard to achieve in manual production runs, especially over large surface areas. For these same reasons we expect a good gain uniformity over large GEM surfaces. An optimal energy resolution in the presence of an uncollimated 55Fe source, as evidenced in Fig. 8, points in this direction. As a matter of fact, first tests of gain uniformity (Fig. 9) yield values already comparable to other MPGDs . Finally, Fig. 10 displays the gas gain measured using the GEMs as an isolated detector, i.e., without a charge collection backpanel anode. We observe no deviation from the expected behavior, nor any anomaly in the onset of discharges (at about VGEM = 450V in Ar + 10% DME and 600V in Ar:CO2).

While the R&D on these GEMs has barely started, all observations are presently very encouraging. First trials with a liquid crystal polymer (LCP) substrate show near-cylindrical hole walls, which can be of interest in applications where excessive dielectric charge-up via ion deposition is a concern (this can lead to a diminished gain uniformity across the surface). Other advantages of LCP compared to KAPTON$`^{\text{TM}}`$  are a much smaller maximum water absorption (0.02% vs. 2%, which may result in lower outgas, of relevance in HEP applications where extreme gas purity is required), better dielectric properties and a higher chemical and heat resistance. The last may result in GEMs more compatible with soldering and operation in commonly used etching detector gases such as CF4, and possibly more resistant to sparking. We expect to be able to report on LCP-GEMs soon.

Hopefully the methods presented here will enable the production of large area MPGD’s. These will be required in large next-generation time-projection chambers, the leading candidate for the tracking system at the next linear collider  and also a possible contender in future underground physics experiments . Proposals for hadron-blind GEM-based detectors in heavy-ion physics programs may similarly benefit .

Tested GEM samples can be obtained from collar@uchicago.edu. JIC and PB would like to thank Q. Guo for his assistance in performing SEM and AFM measurements and T. Witten for helpful discussions. JM and IPJS thank Kirk Arndt and Tom Smith of the Department for Physics at Purdue University for technical support.


format: text

text:
# Core radius evolution of star clusters

## 1 Introduction

Star clusters are vital laboratories for testing theories of stellar dynamics and stellar evolution (e.g. Meylan & Heggie, 1997). The globular clusters orbiting the Milky Way constitute a sample of about $`150`$ evolved stellar systems with ages $`10^{10}`$ years. These provide useful constraints on the long term evolution of clusters but information about the paths which led to these end states can only be inferred via comparisons with the end points of numerical cluster simulations. In contrast, the system of rich star clusters in the Large Magellanic Cloud (LMC) includes clusters of comparable mass to Galactic globular clusters but which span the entire age range from less than $`5\times 10^6`$ years (e.g. R 136) to about $`10^{10}`$ years (e.g. NGC 2019). This allows us to observe massive star clusters at various epochs in their evolution and makes it possible to identify evolutionary trends, under the assumption that the external processes affecting cluster evolution have not changed significantly over the lifetime of the LMC. With this aim in mind, HST project GO7307 (e.g. Beaulieu et al., 1999) obtained deep images of four pairs of LMC clusters. Within each selected pair the clusters are approximately coeval, have similar masses and metallicities and lie at comparable radii from the LMC centre but differ significantly in other properties such as core radius. The pairs range in age from about $`10^7`$ years up to $`10^{10}`$ years. The youngest clusters (NGC 1805 and NGC 1818) are old enough to have expelled all their residual gas and so are pure $`N`$-body systems.

Fig. 1 shows a plot of the observed core radii of star clusters in the LMC as a function of their ages, with the clusters in the GO7307 HST sample identified by name. There is a clear trend that the spread in core radius is an increasing function of cluster age. This trend was first identified by Elson et al. (1989) from ground-based images and was recently confirmed by Mackey & Gilmore (2003) using HST archive data. Mackey & Gilmore (2003) calculated core radii for $`53`$ LMC clusters in a uniform manner with the aim of obtaining an unbiased data set. They discuss in some detail the uncertainties in the data and conclude that selection effects are not responsible for the observed trend. They suggest, as Elson et al. (1989) had earlier argued, that the increased spread in core radius is an evolutionary phenomenon: all clusters formed with relatively small core radii ($`23`$ pc) and subsequently some clusters experienced core expansion while others did not.

Elson et al. (1989) discuss the possibility that inter-cluster variations in the slope of the initial mass function (IMF) could have produced the observed spread in core radius. They point out that mass loss from the increased numbers of massive stars found in clusters with IMFs flatter than the standard Salpeter IMF could lead to core expansion of the observed magnitude on a timescale of about $`10^{8.5}`$ years. As Elson (1991) notes, however, the IMF slope required to explain the largest core radii would lead to disruption of the clusters on a timescale of a few $`\times 10^7`$ years. More directly, based on a detailed observational analysis of the clusters in the GO7307 HST sample, de Grijs et al. (2002b) have shown that the slopes of their luminosity functions are similar to each other within the observational uncertainties, strongly suggesting that their IMFs were very similar down to stellar masses of about $`0.81.0`$M. The absence of inter-cluster population variations implies that stellar evolutionary mass loss cannot be invoked to explain the observations.

As a star cluster evolves, encounters lead to a redistribution of energy among the stars. Heavier stars and binaries acquire lower velocities, on average, and sink towards the cluster core. This process of mass segregation could enhance the heating effect of a population of primordial binaries because the binaries would tend to fall to the centre where encounters occur more frequently. Observations of young clusters (e.g. Hillenbrand, 1997) and numerical simulations of star formation (e.g. Bonnell et al., 2001) show that so-called primordial mass segregation can also result from the cluster formation process. de Grijs et al. (2002b) have shown that there is strong evidence for mass segregation in all the young and intermediate-age clusters in the GO7307 sample. Further, the amount of mass segregation appears to be similar in all clusters and cannot therefore account for differences in the cluster core radii.

The aim of this paper is to investigate, by means of $`N`$-body simulations, whether dynamical processes could give rise to the observed core radius–age relation. In this context, the pair of clusters NGC 1868 and NGC 1831 are of particular interest. These clusters have very similar ages and metallicities and their masses differ by at most a factor of two (see Table 1 of de Grijs et al., 2002b). However, as Fig. 1 shows, the core radius of NGC 1831 is almost three times larger than that of NGC 1868. Given that the cluster ages are only around $`500`$Myr any process which has affected the core radii must be capable of producing significant evolution even on this relatively short time scale. In this paper, we consider the dynamical effects of (1) variations in the external tidal field experienced by the clusters and (2) differences in the initial numbers of dynamically active binary stars in the clusters. Given that the crossing times at the half-mass radii of these clusters are about $`10^7`$ yrs and the times to orbit the LMC are about the same order of magnitude, it is reasonable to suppose that either of these processes could have begun to produce observable effects on cluster structure on the relevant time scale.

During the past decade, numerical simulations of star clusters have become increasingly realistic owing to dramatic increases in computational speeds and significant improvements in programming algorithms (e.g. Aarseth, 1999). The $`N`$-body code used in this paper incorporates the stellar evolution models of Hurley et al. (2000, 2002). This makes it possible to analyse the output from the simulations in an identical manner to that used for observational data, thus allowing meaningful confrontation of models with observations. This is particularly important when discussing cluster properties such as core radius because the relation between the theoretical and observational definitions of this quantity depend on the conversion between the mass function and luminosity function of the cluster stars. For a real cluster this conversion is non-trivial because it depends on, among other quantities, the cluster metallicity and is further complicated by the possible presence of mass segregation (see e.g. de Grijs et al., 2002a). In this paper, all comparisons between models and observations are made in terms of directly observable properties.

In the interests of clarity, we include here a note on the definition of the core radius. As Casertano & Hut (1985) discuss, theoreticians, numerical simulators and observers assign different meanings to the term “core radius”, reflecting the quantity of greatest relevance to their purpose. In theoretical work, the core radius $`r_\mathrm{c}`$ refers to the natural length scale of the system under study, for example the scale radius of the widely used King (1966) profile. The operational definition of the core radius commonly used in discussions of $`N`$-body simulations is that of the density radius $`r_\mathrm{d}`$. This is the density-weighted average of the distance of each star from the density centre (Aarseth, 2001), where the local density at each star is estimated from the mass within a sphere containing the star’s six nearest neighbours and the density centre is defined as the density-weighted mean position of the stars. Observationally, since the available information is either in the form of cluster surface brightness profiles or luminosities of individual stars instead of mass densities and stellar masses, the core radius $`R_\mathrm{c}`$ corresponds to the radius at which the projected surface brightness falls to half its central value. In the case of a cluster of equal-mass stars in equilibrium, there is generally a well-defined relation between $`r_\mathrm{c}`$ and $`R_\mathrm{c}`$. For example, for a Plummer (1911) model it is straightforward to show that $`R_\mathrm{c}=0.6435r_\mathrm{c}`$. It is important to note, however, that during the evolution of the cluster this relation is not static and that there is no general relation between $`r_\mathrm{d}`$ and either $`r_\mathrm{c}`$ or $`R_\mathrm{c}`$. In this paper, we use the term “core radius” exclusively to refer to the observational quantity $`R_\mathrm{c}`$. In Section 3.4 we will compare the evolution of $`r_\mathrm{d}`$ and $`R_\mathrm{c}`$, in order to emphasize the care which must be taken in comparing observations with simulations.

In Section 2 we describe briefly the $`N`$-body code used in the simulations, as well as the pipeline used to calculate observational properties of the simulated clusters. Section 3 considers the effects of time variation of an external tidal field on the evolution of the cluster core radius while Section 4.1 deals with clusters moving in a steady tidal field but with different initial binary fractions. In Section 4.2 we investigate the combined effect of both mechanisms by studying binary-rich clusters in a time-varying tidal field and the results of these simulations are discussed in the context of the LMC clusters in Section 4.3. Finally in Section 5 we summarise our conclusions.

## 2 $`N`$-body simulations

$`N`$-body simulations are a very valuable tool in the study of star clusters because they allow many different processes to be studied with a minimum of simplifying assumptions. The simulations discussed below were performed using the NBODY4 code (e.g. Aarseth, 1999) running on the HARP-3 and GRAPE-6 special purpose computers (Makino et al., 1997) at the Institute of Astronomy, Cambridge and on a GRAPE-6 board at the American Museum of Natural History. The NBODY4 code uses the 4th order Hermite integration scheme to follow the orbits of the cluster stars and utilises state-of-the-art regularisation schemes (Mikkola & Aarseth, 1993, 1998) for the accurate and efficient treatment of the internal motion of binary and higher-order subsystems. It also incorporates a full set of stellar evolution models, in the form of parameterized functions, which are used to follow the evolution of single stars and binaries (Hurley et al., 2001). Of particular importance to the present study is the treatment of mass loss. Any gas which is lost from an evolving star and is not accreted on to a binary companion is assumed to leave the cluster instantaneously, thereby reducing the depth of the cluster potential well. A realistic treatment of mass loss is therefore essential to allow the study of the evolution of such cluster properties as the core radius.

Within the $`N`$-body code the equations of motion are integrated in scaled units such that $`G=1`$ and the total mass and virial radius of the cluster are also set to unity (Heggie & Mathieu, 1986). For a star cluster in virial equilibrium the initial energy in these units is $`1/4`$ and the crossing time is $`2\sqrt{2}`$. Given the total mass of the cluster in solar masses and an appropriately chosen length scale (which determines the conversion from $`N`$-body length units to physical units) it is straightforward to obtain the conversion factors for time and velocity from $`N`$-body units to Myr and $`\mathrm{km}\mathrm{s}^1`$, respectively.

In order to facilitate direct comparison of the results of our simulations with the observed data, we determine the core radius for each model cluster using the same processing pipeline as that used by Mackey & Gilmore (2003). In summary, the stellar data from the model clusters (in the form of position and luminosity information for each member star) are binned radially to obtain a projected surface brightness profile. The three-parameter model (Elson et al., 1987)

$$\mu (r)=\mu _0\left(1+\frac{r^2}{a^2}\right)^{\frac{\gamma }{2}}$$
(1)
is then fitted to the cluster profile. Here, $`\mu _0`$ is the central surface brightness, $`\gamma `$ determines the fall-off of the surface brightness profile at large radii and $`R_\mathrm{c}`$ is related to the length scale $`a`$ via $`R_\mathrm{c}=a(2^{2/\gamma }1)^{1/2}`$. Full details of the fitting procedure used and discussion of the associated uncertainties are presented in Mackey & Gilmore (2003). In order to reduce the numerical noise due to the small number of stars in our simulated clusters we observe each cluster along three perpendicular lines of sight and obtain the median core radius from these three estimates. This helps to reduce the impact of individual bright stars which might otherwise skew the core radius estimation. Further, each $`N`$-body run is performed four times using different random realisations of the initial cluster model and the results are averaged in order to ensure the robustness of any observed trends.

Initially, all the model clusters are generated according to a Plummer (1911) density law

$$\rho (r)=\frac{\rho _0r_\mathrm{P}^5}{(r_\mathrm{P}^2+r^2)^{5/2}},$$
(2)
following the prescription of Aarseth et al. (1974). Here $`r_\mathrm{P}`$ is the Plummer length scale and $`\rho _0`$ is the central mass density. Observations of the young clusters in the LMC show that they do not appear to be tidally limited (Mackey & Gilmore, 2003). In this case, the Plummer profile is a reasonable model to use for the initial density profile given that the fall-off in its outer regions lies between that of a tidally truncated King model and the observed profiles of the young LMC clusters.

The initial distribution of stellar masses is drawn from the IMF of Kroupa et al. (1993) which is a three-part power law fit to stellar data in the solar neighbourhood. We take the lower and upper mass limits to be $`0.1`$M and $`50`$M, respectively. We assume solar metallicity for the stars in our clusters. This is an appropriate choice for the younger clusters – for example, Johnson et al. (2001) have shown that NGC 1805 and NGC 1818, whose ages are about $`10`$ and $`25`$Myr respectively, are well matched by isochrones of approximately solar metallicity. Since the pairs of clusters in the GO7307 HST sample were chosen to have similar metallicity, we know a priori that the observed range of cluster core radii cannot be due to metallicity variations. We have therefore chosen to perform all our simulations at solar metallicity, and postpone discussion of the effects of metallicity variations on star cluster evolution to future work.

Table 1 summarises the simulations which we have performed. The clusters in these runs initially contain $`4500`$ or $`5000`$ stars and have a mean mass of about $`2300`$ M. Clearly, the LMC clusters are more massive and contain greater numbers of stars than our model clusters and therefore direct comparison of the results of our simulations with the observed data must be approached with caution. However, the effects of increasing $`N`$ in our simulations can be estimated by considering the various time scales in the problems we study. Detailed discussion of these scaling issues is presented in Sections 3.4 and 4.3 where we argue that none of our conclusions is strongly dependent on the small size of our simulated clusters.

## 3 Time-varying tidal fields

### 3.1 Model Initialisation

We first consider the effects of the external tidal field of the LMC on the evolution of the cluster core radii. Theoretical arguments and numerical simulations strongly suggest that tides have played a significant role in the evolution of star clusters and cluster systems in galaxies such as the Milky Way (e.g. Gnedin & Ostriker, 1997; Vesperini & Heggie, 1997; Vesperini, 1998; Baumgardt, 1998; Fall & Zhang, 2001) . In addition, observations of tidal tails around at least twenty Milky Way globular clusters (e.g. Leon et al., 2000; Odenkirchen et al., 2001) provide compelling evidence of the ongoing disruption of star clusters by the Milky Way. The gravitational field of the LMC (mass $`10^{10}`$M) is much weaker than that of a galaxy like the Milky Way (mass $`10^{12}`$M). However, as Mackey & Gilmore (2003) show (their Fig. 15), there is a deficit of low-mass clusters older than $`10^{10}`$ years. If we assume that the cluster population originally contained a spectrum of cluster masses similar to that of the present-day young clusters, their absence at late times strongly suggests that tides have been at work on the cluster system in exactly the manner expected from the models of Gnedin & Ostriker (1997). Given that the GO7307 HST sample of 8 clusters was selected on the basis that they lie at similar distances from the LMC centre but display a wide range of core radii, the steady tidal field experienced by a cluster on a roughly circular orbit cannot be the origin of the observed trends. Further, given the intrinsic weakness of the LMC field at a radius of $`4`$kpc, where our clusters lie, one might not expect the tidal field to play a significant role in cluster evolution. We therefore consider clusters on non-circular orbits. Such clusters experience tidal shocks during their orbits and we investigate whether these impart sufficient energy to the cluster to produce differences in the evolution of $`R_\mathrm{c}`$Elson (2000) discusses the possible role of disk shocking in the internal evolution of the intermediate-age LMC clusters, but questions whether any of the clusters is on an appropriate disk-crossing orbit for this to be an important effect. In this section we consider instead the tidal shocking produced by a point mass galaxy.

In order to investigate the role of the LMC tidal field in the evolution of its cluster system, it is necessary to incorporate an appropriate external gravitational field into the simulations. We model the LMC potential as that of a point mass of $`9\times 10^9`$M. This is a rather drastic oversimplification of the actual situation. We note, however, that this model exaggerates the effect of the tidal field, as any extended distribution of mass will have a more slowly varying gravitational potential. Further, in the radius range $`28`$ kpc, in which our model clusters orbit, the gradient of the potential is within a factor of two of that of the detailed LMC mass model of van der Marel et al. (2002). Thus our simulations constitute a plausible upper limit to the possible effects of tides.

In the $`N`$-body simulations, the equations of motion are integrated in an accelerating but non-rotating reference frame centred on the centre of mass of the cluster and whose axes are aligned with fixed spatial directions. This choice of frame simplifies the equations of motion for a cluster on a non-circular orbit despite having less symmetry than the more usual rotating frame used for steady tidal fields whose axes always point towards the attracting centre (Heggie, 2001). Implementation of the time varying tidal field requires some care because mass loss from evolving stars or the escape of stars from the cluster causes the centre of mass to acquire a non-zero velocity relative to the frame of integration. It is therefore necessary to re-centre the frame of reference to coincide with the rest frame of the centre of mass following each such event. All internal and external forces must then be reinitialized to ensure consistency with the re-centred frame during subsequent integrations.

We have performed simulations for clusters on elliptical orbits with perigalactic radii $`R_\mathrm{p}2`$kpc and apogalactic radii $`R_\mathrm{a}8`$kpc, probing the region in which most of our sample clusters lie but largely avoiding the complications of the inner regions of the LMC which may contain a stellar bar (e.g. Zhao & Evans, 2000). We also modelled clusters on circular orbits at a radius of $`6`$kpc which corresponds to the time-averaged radius of the elliptical orbits. This allows us to isolate the effect of the time variation of the tidal field because both sets of clusters experience the same mean tidal field.

We calculate the initial perigalactic tidal radius $`r_{\mathrm{t},\mathrm{p}}`$ in the approximation that the cluster is on a circular orbit at a radius $`R_\mathrm{p}`$. In this case, $`r_{\mathrm{t},\mathrm{p}}`$ is determined from the cluster mass $`M`$, the mass of the point mass galaxy $`M_\mathrm{G}`$ and $`R_\mathrm{p}`$ via (e.g. King, 1962)

$$r_{\mathrm{t},\mathrm{p}}=R_\mathrm{p}\left(\frac{M}{3GM_\mathrm{G}}\right)^{1/3}.$$
(3)
The apocentric tidal radius $`r_{\mathrm{t},\mathrm{a}}`$ is determined in a similar way from $`R_\mathrm{a}`$. The unit of length, which determines the conversion between $`N`$-body length units and physical units, is then chosen such that the limiting radius of all clusters is initially equal to $`r_{\mathrm{t},\mathrm{p}}`$. However, at the start of the simulations the clusters are placed at the apocentre of their orbits about the LMC which means that they do not fill their tidal radii. If, instead, the clusters initially extended to $`r_{\mathrm{t},\mathrm{a}}`$ and were not truncated at $`r_{\mathrm{t},\mathrm{p}}`$ then the removal of stars outside $`r_{\mathrm{t},\mathrm{p}}`$ during the first few perigalactic passages could mask any heating effects at radii interior to $`r_{\mathrm{t},\mathrm{p}}`$ due to the time variation of the tidal field. In contrast to a King model, the Plummer density law (2) formally extends to infinity and so the outer radius of a cluster is not well-defined. The length unit we use yields a pericentric tidal radius close to the (random) position of the outermost star in our initial clusters. The internal stellar velocities are scaled to ensure virial equilibrium at the start of the simulation.

The tidal radius of a cluster on an elliptical orbit does not have a unique value. It varies continuously both as the cluster moves along its orbit and when the cluster loses mass through stellar evolution or the escape of member stars. However, for clusters whose internal relaxation times are longer than their orbital periods, tidal stripping tends to truncate the star cluster at $`r_{\mathrm{t},\mathrm{p}}`$ (Oh & Lin, 1992). These two time scales are of similar magnitude for our simulated clusters so we expect two-body relaxation processes to lead to some expansion of the cluster beyond $`r_{\mathrm{t},\mathrm{p}}`$. We take account of this by defining escaping stars to be those stars outside twice the apogalactic tidal radius $`r_{\mathrm{t},\mathrm{a}}`$. Stars whose orbits carry them outside this radius are generally lost to the cluster (e.g. Heggie, 2001). Stars are thus free to populate the region between $`r_{\mathrm{t},\mathrm{p}}`$ and $`2r_{\mathrm{t},\mathrm{a}}`$.

The choice of length unit determines the physical density of the cluster and hence the physical rate at which two-body relaxation processes affect the evolution of the cluster. Our aim is to investigate the influence of a time-varying external tidal field so it is important to ensure that all internal evolutionary processes proceed in an identical manner for each of our simulated clusters. We therefore use the same length unit for the clusters on circular orbits as we chose for those on elliptical orbits. This means that the clusters on circular orbits also do not initially fill their tidal radii and so necessarily allows expansion of these clusters. However, since the circular orbits are placed at the time-averaged radius of the elliptical orbits this expansion should on average be the same for both sets of clusters. Any differences in the rates of expansion are due to the time variation in the tidal field.

The effect of steady tidal fields on star cluster evolution has been the subject of papers by a number of authors (e.g. Giersz & Heggie, 1997; Aarseth & Heggie, 1998; Baumgardt, 2001; Portegies Zwart et al., 2002). However, until recently the problem of clusters moving in time-varying tidal fields had received relatively little attention from $`N`$-body simulations. Baumgardt (1998) uses the results of some simulations of single-mass star clusters orbiting a point mass to determine the appropriate choice of escape radius for simulations of mass loss. More recently, Baumgardt & Makino (2003) presented an extensive survey of the dynamical evolution of star clusters in external tidal fields. Our investigation is complementary to theirs because we are looking at the evolution of the observable global structure of the cluster while Baumgardt & Makino (2003) focus their attention mostly on the evolution of the stellar mass function within the cluster. The results of all our simulations are presented in terms of quantities which can be compared directly with those calculated from the observed data, thereby enabling us to discuss how different dynamical processes impact on observations.

We follow the evolution of the clusters for about $`5t_{\mathrm{rh}}`$ (where $`t_{\mathrm{rh}}`$ is the half-mass relaxation time) This corresponds to $`T_{\mathrm{max}}1`$Gyr in most cases. By this time, the clusters on elliptical orbits have lost more than $`50`$ per cent of their stars making estimates of properties such as the core radius increasingly noisy. At regular intervals during the evolution, the data from the clusters are recorded and all relevant properties of the clusters are calculated.

### 3.2 Results

The solid curve in Fig. 2 presents the evolution of the core radius for clusters on circular orbits about the LMC at a distance of about $`6`$kpc from the centre. The broken curve shows the evolution of $`R_\mathrm{c}`$ for clusters on elliptical orbits with the same time-averaged galactocentric radius. It is apparent from this figure that the core radius evolution in both cases is indistinguishable, especially in view of the noise associated with the output of small-$`N`$ simulations. During the first $`10^8`$yrs, the core radius is seen to decrease slightly from its initial value. This is because the calculation of $`R_\mathrm{c}`$ is dominated by the light of the most massive stars which rapidly become more centrally concentrated due to the effects of mass segregation. As Fig. 3 shows, after $`10.5`$Myr the half-mass radius of the more massive stars is already $`30`$ per cent smaller than that of the less massive stars.

After its initial contraction phase, the core radius begins to increase. This occurs because the clusters do not fill the tidal radius corresponding to their galactocentric positions at the start of the simulations. A combination of mass loss from the rapid evolution of massive stars and two-body relaxation can therefore drive the overall expansion of the cluster seen in Fig. 3. The top panel of Fig. 4 shows that the cluster’s half-mass radius $`r_\mathrm{h}`$ is increasing with time (see fourth panel of Fig. 5 for long-term $`r_\mathrm{h}`$ behaviour). During this period its $`r_\mathrm{d}`$ (not shown) is roughly constant at about $`0.3`$pc.

The expansion of $`r_\mathrm{h}`$ is driven by a combination of two processes which are illustrated in the lower panels of Fig. 4. The middle panel shows the evolution of the virial ratio $`Q`$. Initially the cluster is in equilibrium with $`Q=0.5`$. After about $`10`$Myr, the cluster is temporarily displaced from equilibrium ($`Q>0.5`$) by mass loss from the most massive stars. This is accompanied by the increase in the half-mass radius of the cluster shown in the top panel of Fig. 4. It is interesting to note that the $`R_\mathrm{c}`$ of the cluster does not change as a result of this mass loss. Stellar evolution cannot account for the sustained expansion of the cluster, however. The cluster loses about $`5`$ per cent of its mass as a result of the above mass loss events. Using the virial theorem it is straightforward to show that the gravitational radius $`R^{}`$ of a stellar system with initial radius $`R`$ and mass $`M`$ which loses an amount of mass $`m`$ is given by (Hills, 1980)

$$R^{}=\frac{1}{2}\left(\frac{Mm}{\frac{1}{2}Mm}\right)R.$$
(4)
A $`5`$ per cent mass loss should lead to an increase of less than $`10`$ per cent in the size of the cluster which is insufficient to account for the cluster expansion seen in Figs. 2 and 3.

The bottom panel of Fig. 4 shows the evolution of the velocity anisotropy parameter $`\beta =1v_\mathrm{t}^2/2v_\mathrm{r}^2`$ (Binney & Tremaine, 1987) inside the $`10`$ per cent Lagrangian radius of the cluster (the radius containing $`10`$ per cent of the total mass of the cluster) and in the outer parts (between the $`80`$ per cent and $`90`$ per cent Lagrangian radii). Initially, the velocity distribution throughout the cluster is isotropic ($`\beta =0`$). While the central parts of the cluster retain this isotropy, in the outer parts the distribution becomes increasingly radially anisotropic ($`\beta >0`$). This indicates that the continued growth of $`r_\mathrm{h}`$ is being fueled by encounters in the inner parts of the cluster which produce a population of stars on extended radial orbits. The velocity distribution in the outer parts of isolated clusters simulated by Baumgardt et al. (2002) exhibited similar behaviour during post core-collapse evolution; those authors similarly concluded that the growth of the outer Lagrangian radii at late times was due to two-body processes in the inner regions.

Given the similarity between the core radius evolution of clusters on circular and elliptical orbits, we now proceed to investigate the extent to which the clusters on elliptical orbits have been affected by the variation in the tidal field. Fig. 5 presents the evolution of certain key internal properties of two typical clusters, one on a circular orbit and one on an elliptical orbit. The top panel shows the galactocentric radius $`R_\mathrm{g}`$ along the elliptical orbit to emphasize the clear correlation between perigalactic passages and changes in the internal properties of the cluster on the elliptical orbit. The second panel shows that, as expected, during a perigalactic passage the rms velocity $`v_{\mathrm{rms}}`$ of the stars in the cluster increases sharply. This velocity kick ejects the most energetic stars from the cluster. As a result of this mass loss, $`v_{\mathrm{rms}}`$ settles to a value slightly below that of the cluster on the circular orbit once it returns to virial equilibrium.

The third and fourth panels of Fig. 5 show how $`t_{\mathrm{rh}}`$ and $`r_\mathrm{h}`$ evolve during the cluster orbit. The most noticeable effect is that the relaxation time for the cluster on the elliptical orbit is significantly reduced after the second perigalactic passage, so that the cluster evolves more rapidly following this encounter. This is due both to the decrease in the total number of stars in the cluster and to the increase in the mean stellar mass $`\overline{m}`$ (shown in the bottom panel of Fig. 5). The increase in $`\overline{m}`$ occurs because, on average, low-mass stars acquire larger velocities during pericentre passages and are hence more likely to escape from the cluster. The value of $`r_\mathrm{h}`$ is seen to increase during the perigalactic passage but rapidly settles back to a value slightly below the circular-orbit value. Thus while the tidal shocks experienced by the clusters on elliptical orbits do not affect the evolution of $`R_\mathrm{c}`$, they are nonetheless strong enough to modify other important properties of the cluster. As Fig. 5 shows, the impact of the shocks on the cluster properties increases with time – this is due to the decreasing number of bound stars in the cluster.

### 3.3 Strength of tidal shocks

Before drawing our final conclusions based on these simulations, we now investigate whether the lack of impact on $`R_\mathrm{c}`$ is due to the weakness of the shocks we have simulated or whether the value of $`R_\mathrm{c}`$ is intrinsically robust to tidal shocking. The strength of a tidal shock is generally presented in terms of the dimensionless parameter $`\beta _{\mathrm{sh}}`$ given by (e.g. Spitzer, 1987)

$$\beta _{\mathrm{sh}}=\frac{t_{\mathrm{rh}}}{t_{\mathrm{sh}}},$$
(5)
where $`t_{\mathrm{sh}}`$ is an estimate of the destruction time of the cluster owing to energy input from the gravitational shocks. The time scale $`t_{\mathrm{sh}}`$ is given by

$$t_{\mathrm{sh}}=\frac{|E_\mathrm{h}|}{dE_\mathrm{h}/dt},$$
(6)
where $`E_\mathrm{h}`$ is the energy of the cluster evaluated at the half-mass radius and $`dE_\mathrm{h}/dt`$ is the rate at which stars at the half-mass radius acquire energy due to the shock. In a simulation containing a finite number of stars, evaluation of this expression requires some care. Averaging the energy changes of stars near the half-mass radius does not produce robust results because individual stars may gain and lose energy during the perigalactic passage through other means than interaction with the external tidal field. Instead we calculate the total energy binding the cluster as a function of time and estimate the energy input to the cluster during perigalactic passages. Fig. 6 compares the time evolution of the binding energy $`E_\mathrm{b}`$ of a cluster on a circular orbit with that of a cluster on an elliptical orbit. For a cluster in an external potential $`E_\mathrm{b}`$ is defined as (e.g. Fukushige & Heggie, 1995)

$$E_\mathrm{b}=\mathrm{T}+V+\mathrm{E}_\mathrm{T}$$
(7)
where $`T`$ and $`V`$ are the total kinetic and potential energy of the cluster, respectively, and $`E_\mathrm{T}`$ is the energy in the tidal field. The binding energy must be negative for the cluster to remain bound. The effects of the external field are visible as dramatic changes in $`E_\mathrm{b}`$ at about $`180`$ and $`520`$Myr. Note that the sudden, temporary change in $`E_\mathrm{b}`$ for the circular run at about $`700`$Myr is due to a collision between two binaries that results in the the break-up of one of the binaries and the subsequent ejection of its components. The ratio of the energy produced by this event to the binding energy of the cluster emphasises the relatively fragile nature of the clusters we are simulating.

During perigalactic passages, some stars acquire velocities above the escape velocity of the cluster and this leads to an enhanced mass-loss rate (see Fig. 7). However, the finite time required for such stars to escape means that $`E_\mathrm{b}`$ can temporarily become positive because it is calculated by summation over all stars within twice the tidal radius of the cluster. After the cluster has passed pericentre, stars continue to leave across the tidal boundary and $`E_\mathrm{b}`$ returns to a negative value indicating that the remaining cluster is still gravitationally bound. Following each shock, the cluster undergoes a revirialisation process during which the energy acquired from the tidal interaction is redistributed among the cluster stars. This process results in further loss of mass from the cluster because, on average, the revirialisation process gives velocity kicks to stars which are three times larger than those arising from the initial shock (see e.g. Gnedin & Ostriker, 1999).

We take the energy input by the shock to be the difference between the energy of the cluster after this revirialisation process is complete and the energy of the cluster on the circular orbit at this time (about $`T=540`$Myr). The time between the shocks is $`340`$Myr, which leads to $`t_{\mathrm{sh}}6800`$Myr. At this time $`t_{\mathrm{rh}}`$ for the cluster is about $`170`$Myr, so the shock strength is about $`0.02`$. Comparing this value with Fig. 13 of Gnedin et al. (1999) and noting that the initial concentrations $`c=\mathrm{log}_{10}(r_\mathrm{t}/r_\mathrm{d})`$ of our clusters are about $`1.1`$ (where we assume that $`r_\mathrm{t}`$ is the initial outer limit of the clusters), we find that our clusters lie approximately on the borderline of the region where shocks are important. We conclude that the shocks our clusters are experiencing would significantly affect their evolution on long time scales. Indeed, as Fig. 8 illustrates, even the $`30`$ per cent Lagrangian radius of the cluster is affected by the tidal shocks, indicating energy input into the cluster over a wide range of radii. Nevertheless, on a timescale of up to $`10^9`$ yrs, the effects of tidal heating are insufficient to modify the observed $`R_\mathrm{c}`$ of the clusters significantly.

### 3.4 Discussion

The top panel of Fig. 4 illustrates an important point which is worth emphasizing in the context of comparisons between simulated data and observations, namely that $`R_\mathrm{c}`$ and $`r_\mathrm{h}`$ are not, in general, related in a simple way. In particular circumstances, for example at the start of our simulations where the model clusters are in strict dynamical equilibrium, it is easy to show that, for a Plummer model composed of equal mass stars, the two quantities are related by $`r_\mathrm{h}=2R_\mathrm{c}`$. However, Fig. 4 shows that this relation quickly breaks down as the cluster evolves away from its initial conditions. Thus, in the absence of further information, observational estimation of $`R_\mathrm{c}`$ must not be assumed to imply knowledge of the mass distribution in the cluster.

The same is true of the density radius $`r_\mathrm{d}`$ defined in Section 1. Fig. 9 shows the evolution of $`R_\mathrm{c}`$ (solid curve) and $`r_\mathrm{d}`$ (broken curve) for one of the clusters on an elliptical orbit. As the top panel of the figure shows, in a small-$`N`$ simulation, both $`R_\mathrm{c}`$ and $`r_\mathrm{d}`$ are intrinsically noisy quantities: this is emphasised by the middle panel where we have smoothed the curves from the top panel using a window function of width five time bins. The bottom panel of Fig. 9, which shows the ratio of $`R_\mathrm{c}`$ to $`r_\mathrm{d}`$ for the smoothed data, demonstrates that the relation between the two quantities is not a simple proportionality and that the evolutionary trends in the two quantities can be quite different. The complicated relation between light and mass in star clusters is the result of a variety of processes such as stellar evolution and mass segregation which modify the mass and light distributions in different ways. Given the absence of simple relations between the different core radii, great care must be exercised when we interpret observational data through comparison with numerical simulations or theoretical models.

The simulations we have performed contain relatively low numbers of stars and are smaller and of lower total mass than the actual LMC clusters. However, the mean densities of our clusters are very similar to those of the LMC clusters and it is this which determines the extent to which an external tidal field can influence the internal dynamics of a cluster. On this basis, we can be confident that conclusions based on our simulations are relevant to the clusters in the LMC. Our simulations have shown that, on a timescale of $`1`$Gyr, which corresponds to several $`t_{\mathrm{rh}}`$ for our model clusters, tidal forces are not able to produce a spread in core radius similar to that observed among the LMC clusters.

There are two obvious avenues for building on the simulations we have performed. First, simulations of clusters containing a larger numbers of stars $`N`$ could be performed. This would allow us to achieve cluster masses comparable to those of the real LMC clusters. However, increasing $`N`$ in a cluster while holding the physical density of the cluster constant leads to an increase in $`t_{\mathrm{rh}}`$. The expansion seen in $`R_\mathrm{c}`$ for our simulated clusters is driven primarily by two-body processes (which proceed on a time scale $`t_{\mathrm{rh}}`$) rather than by tidal heating and it will therefore be slowed by increasing $`N`$. Further, the disruptive effect of pericentre passages depends on the ability of a cluster to re-fill the region outside its pericentric tidal radius $`r_{\mathrm{t},\mathrm{p}}`$ during the time between close passages. As $`t_{\mathrm{rh}}`$ increases this re-filling is reduced. The cluster will therefore have fewer stars outside $`r_{\mathrm{t},\mathrm{p}}`$ and hence will be less affected by shocks than small-$`N`$ clusters (see Fig. 2 of Baumgardt & Makino (2003)). A cluster on an elliptical orbit and which, as in our simulations, initially filled its pericentric tidal radius would remain truncated at this radius throughout its evolution because there would not be sufficient time between perigalactic passages for the region outside $`r_{\mathrm{t},\mathrm{p}}`$ to be re-filled (Oh & Lin, 1992). Thus shocks should have less impact on the $`R_\mathrm{c}`$ of larger-$`N`$ clusters and the observed expansions of clusters on both circular and elliptical orbits should remain similar.

The second possible improvement would be to use a more realistic model of the LMC potential, for example one based upon the recent mass models of van der Marel et al. (2002). However, the actual LMC mass distribution is more extended than the Keplerian potential we have been using in our simulations and therefore the gradients in its gravitational field are smaller. Consequently, the associated strength of its tidal effects on clusters moving on elliptical orbits is reduced.

Thus, we conclude that neither of the above modifications to our simulations would tend to increase the effect of the tidal field. The fact that the time-varying tidal field in our simulations did not result in changes in the evolution of $`R_\mathrm{c}`$ in the extreme case of a low-mass cluster in orbit about a point-mass galaxy leads us to conclude that tides have not yet had time to play a significant role in the internal evolution of the intermediate-age LMC clusters.

## 4 Models with primordial binaries

### 4.1 Steady tidal field case

In this section we present results for models moving in a steady tidal field but with different numbers of primordial binaries, in order to isolate the impact of large binary fractions on the evolution of the cluster properties. The clusters discussed here move on the same orbits as the circular orbit runs of Section 3.

Based on an analysis of the observed binary sequence in the young cluster NGC 1818, Elson et al. (1998) have shown that the binary fraction $`f_\mathrm{b}`$ in the centre of this cluster is about $`0.35`$. Here $`f_\mathrm{b}=N_\mathrm{b}/(N_\mathrm{s}+N_\mathrm{b})`$,where $`N_\mathrm{s}`$ and $`N_\mathrm{b}`$ are the numbers of single stars and binary systems, respectively. To date, the binary fractions of the older clusters in our sample have not been determined. Thus, at present the binary fractions of these clusters are essentially free parameters. Heggie’s (1975) law implies that only hard binaries are dynamically important in the evolution of a star cluster. A binary is said to be hard when the magnitude of its binding energy $`|E_{\mathrm{bin}}|`$ exceeds the kinetic energy of a typical cluster star (Heggie, 1975). In the case of a binary with component masses $`m_1`$ and $`m_2`$ at a separation $`a`$ moving in a cluster with mean mass $`\overline{m}`$ and velocity dispersion $`\sigma `$ the hardness condition is

$$\frac{Gm_1m_2}{2a}>\frac{1}{2}\overline{m}\sigma ^2.$$
(8)
(We note that the NBODY4 code uses a dispersion of $`\sqrt{2}\sigma `$ in this equation to determine the hardness limit for binaries, in order to take account of the larger velocity dispersion in the core of a centrally concentrated cluster.) In our simulations the boundary between hard and soft binaries occurs at a separation of about $`120`$AU for binaries containing stars of the mean stellar mass in the cluster. Hard binaries act as heat sources because they become more tightly bound as a result of interactions with single stars and thereby impart kinetic energy to the stars they encounter. Heggie & Aarseth (1992) found that the density radii of clusters containing primordial binaries experienced shallower contraction during core collapse and subsequently tended to be larger than those of clusters containing only single stars. This result, based on simulations of single-mass clusters with binary fractions of about $`0.03`$, motivates our investigation of clusters containing large numbers of primordial binaries. Hereinafter, we will use $`f_\mathrm{b}`$ to refer to the fraction of *hard* binaries in the cluster.

As an extreme case, we consider first clusters with initial $`f_\mathrm{b}=0.5`$ and with the hard binaries randomly distributed throughout the cluster. In each binary the mass of the primary star is drawn from the IMF of Kroupa et al. (1993) and the mass of its companion is obtained by randomly choosing a mass ratio for the binary from a uniform distribution constrained by the single star mass limits (for a more detailed discussion, see Hurley et al., 2001). The binary orbital parameters are generated from the distribution of Eggleton et al. (1989). The clusters are otherwise generated exactly as in the single star runs of Section 3. We note that the inclusion of a large number of hard primordial binaries dramatically increases the computation time required for each simulation.

Fig. 10 compares the evolution of $`R_\mathrm{c}`$ for clusters containing $`N=5000`$ single stars and no primordial binaries with that for clusters containing a significant initial binary population. An increased binary fraction produces a noticeable effect on the evolution of $`R_\mathrm{c}`$. After $`1`$Gyr, the core radii of clusters containing $`50`$ per cent primordial binaries are about a factor of two larger than those of clusters with no binaries.

The results presented in Fig. 11 give some insight into the mechanism by which a large binary fraction can drive an expansion of $`R_\mathrm{c}`$. The top panel shows that the rate at which mass is lost from a cluster with initial $`f_\mathrm{b}=0.5`$ is significantly greater than in the case of a cluster with no initial binary population. After $`1`$Gyr this cluster has lost almost $`25`$ per cent more mass than the cluster containing only single stars. The evolution of the energy binding the cluster $`E_\mathrm{b}`$ is shown in the second panel of Fig. 11. In a cluster with a large initial binary fraction, $`E_\mathrm{b}`$ rapidly attains a value significantly higher than that of a cluster of single stars. The enhanced mass-loss rate due to the presence of the binaries decreases the depth of the cluster’s potential-well; the magnitude of the binding energy is further reduced by the heating of the remaining bound stars through weak encounters. As Fig. 12 shows, the energy input due to binary heating is distributed throughout the cluster. All the Lagrangian radii of the cluster rapidly attain larger values when the cluster contains primordial binaries. We note that, as was discussed in Section 3.2, the growth of the clusters containing only single stars is due to the fact that the clusters do not fill their tidal radii at the start of the simulations. We use identical initial conditions for the clusters containing binaries – the excess heating observed in those clusters is due to the presence of the hard binary population.

The third panel of Fig. 11 shows the internal binding energy of the primordial hard binaries $`E_{\mathrm{bin}}`$ as a function of time. In calculating $`E_{\mathrm{bin}}`$ we include the binding energy of those hard binaries which have escaped from the cluster because we are interested in the rate at which encounters lead to the hardening of binaries. This plot demonstrates that, as expected, the hard binaries become harder with time and so impart energy to the cluster. It is well known that energy generation by hard binaries halts core collapse (e.g. Heggie & Aarseth, 1992). However, Fig. 11 demonstrates that even in the very early stages of a cluster’s evolution, the heating effect of a large primordial hard binary population can have a significant impact on the energy budget of the cluster.

We note that the occasional large changes in $`E_{\mathrm{bin}}`$ are the result of stellar evolutionary processes and therefore do not contribute to the direct heating of the cluster. For example, at $`T=34`$Myr the components of one binary system merge to form a single star: prior to the merger the binary is hardened by tidal interactions between its components. The massive star thus produced explodes as a supernova soon after and the velocity kick received by the remnant as a result of the explosion is reflected in the spike in $`E_\mathrm{b}`$. The prominence of this event is due to the relatively small number of stars in the simulation.

Fig. 10 shows that the $`R_\mathrm{c}`$ evolution of clusters with $`f_\mathrm{b}=0.25`$ is very similar to that of clusters with $`f_\mathrm{b}=0.5`$ initially. This suggests that the heating effect of primordial binaries saturates above a certain value of the binary fraction $`f_{\mathrm{b},\mathrm{crit}}`$. For larger values of $`f_\mathrm{b}`$, the $`R_\mathrm{c}`$ evolution is unchanged, although the mass-loss rate is greater for increased $`f_\mathrm{b}`$. From Fig. 10 we see that $`f_{\mathrm{b},\mathrm{crit}}`$ for the clusters we have simulated appears to lie in the range $`0.10.25`$ because an initial binary fraction of $`10`$ per cent has a noticeably smaller impact on $`R_\mathrm{c}`$ than one of $`25`$ per cent. We are currently investigating this effect in more detail, including its dependence on the number of stars in the simulation: the results will be presented elsewhere. Here we note only that the saturation may arise from the fact that at early times only a relatively small number of binaries take part in the heating process. In our simulations only about three per cent of the binaries are dynamically hardened during the first $`100`$Myr. We find that this fraction appears to be independent of the value of $`f_\mathrm{b}`$. Whether a particular binary experiences an encounter which leads to hardening depends on the masses and separation of the components as well as the local stellar density in the neighbourhood of the binary. For a fixed binary mass, however, the hardest binaries have the smallest separations and hence are least likely to interact with other stars. This is reflected in our simulations: for a given binary mass, all the binaries which undergo early hardening are initially among the more weakly bound pairs.

### 4.2 Time-varying tidal field case

In our final set of simulations we examine clusters containing large numbers of primordial binaries and moving on non-circular orbits about the LMC. We have considered the extreme case of $`f_\mathrm{b}=0.5`$ but, as was shown in the previous section, these results also apply to somewhat smaller values of $`f_\mathrm{b}`$. As Fig. 13 shows, the core radius of a cluster on an elliptical orbit and containing primordial binaries exhibits stronger growth than that of a cluster of single stars on a circular orbit: after about $`600`$Myr the core radius of the former cluster is almost a factor of two larger. The increase in core radius is solely due to the presence of binaries. As in the case of clusters containing only single stars (Section 3), a cluster with the same primordial binary fraction but on a circular orbit exhibits almost identical core radius evolution. However, the increased mass-loss rate which results from binary heating leaves the cluster on the elliptical orbit more vulnerable to tidal disruption. After only $`600`$Myr, the cluster has lost almost $`75`$ per cent of its mass and the core radius determination is consequently very uncertain. By $`1`$Gyr it is virtually disrupted and retains only about $`10`$ per cent of its initial mass while the circular orbit cluster still contains about $`45`$ per cent of its initial mass.

### 4.3 Discussion

Our simulations demonstrate that the presence of a large population of hard binaries can produce significant expansion of the core radius of a cluster. In agreement with many previous studies (e.g. Spitzer & Mathieu, 1980; Heggie & Aarseth, 1992; de La Fuente Marcos, 1996) we have thus shown that binaries constitute a dynamically significant population in star clusters. In particular, we have found that primordial binaries can rapidly affect observable cluster properties such as $`R_\mathrm{c}`$. In the context of the LMC clusters, however, it is clear that the presence of binaries alone cannot account for the the observed trends in $`R_\mathrm{c}`$ and in particular for the factor of three difference in core radius between NGC 1831 and NGC 1868. Even in the extreme case in which one cluster has no primordial binaries at all while the other cluster has a hard binary fraction of $`2550`$ per cent, the difference in $`R_\mathrm{c}`$ over the lifetime of the intermediate-age clusters is at most a factor of two. The core radius of a cluster containing $`10`$ per cent binaries expands by only a factor of $`1.5`$ over the same time scale. Further, the similarity of the observed present-day luminosity functions of the clusters in the GO7307 sample (de Grijs et al., 2002b) argues against large inter-cluster variations in the IMF and hence makes significant differences in binary fraction unlikely.

The initial radial distribution of hard binaries has an impact on their heating effect on the cluster stars (e.g. de La Fuente Marcos, 1996) and may also have a bearing on the saturation of the heating effect which we observed in Section 4.1. A binary population which was more centrally concentrated than the single stars as a result of either primordial or dynamical mass segregation would be expected to produce more heating because of the increased interaction rates in the inner parts of the cluster. In our future simulations we will consider alternatives to the uniform distribution used in the present runs. However, based on our current simulations it is possible to estimate the significance of centrally concentrating the binary population. The clusters with $`50`$ per cent binary fractions contain approximately the same number of binaries within $`R_c`$ as a cluster with a $`10`$ per cent overall binary fraction in which all the binaries are located within $`R_c`$. Given that none of our simulations exhibited $`R_c`$ growth of a magnitude sufficient to match the observations, a core binary fraction in excess of $`50`$ per cent would be required to produce greater $`R_c`$ growth. The observed similarity of the degree of mass segregation in the intermediate-age clusters (de Grijs et al., 2002b) render it unlikely that the radial distribution of binaries could vary so greatly between clusters.

As in the case of the simulations in Section 3, it is natural to ask whether our conclusions are also valid for clusters containing larger numbers of stars than those in our current simulations. The simulation of large-$`N`$ clusters containing significant numbers of primordial binaries has only recently become feasible and we are currently investigating the impact of binaries on the evolution of such clusters; the results of this study will be presented elsewhere. The present simulations show that the disruption time for a small-$`N`$ cluster in a time-varying external tidal field is significantly reduced if it contains a large fraction of primordial binaries. It is not trivial to extrapolate these results to larger-$`N`$ clusters. However, binary heating of the cluster is essentially a two-body process and therefore proceeds on the relaxation time scale $`t_{\mathrm{rh}}`$. Given that a larger-$`N`$ cluster with the same mean physical density has a larger $`t_{\mathrm{rh}}`$, the physical rate at which the cluster is heated by its binary population will be reduced. Thus, while larger-$`N`$ clusters would obviously survive for longer than those in our present simulations, their physical rate of $`R_\mathrm{c}`$ expansion would be reduced. Large variations in primordial binary fraction are therefore unlikely to be the explanation of the LMC observations.

## 5 Conclusions

The core radii of star clusters in the LMC display a strong trend with age, namely that the older clusters exhibit a much greater spread in core radius than their younger counterparts. In particular, the $`R_\mathrm{c}`$ of the intermediate-age clusters NGC 1868 and NGC 1831 differ by about a factor of three. It has been suggested (Elson et al., 1989; Mackey & Gilmore, 2003) that the trend represents a physical evolution of the clusters and that the clusters all formed with similar, small core radii which were modified by subsequent dynamical evolution. The possibility that this might be caused by variations in the cluster IMF has been ruled out based on HST photometry of a sample of eight clusters (de Grijs et al., 2002b). In this paper we have used $`N`$-body simulations of star clusters initially containing up to $`N=5000`$ stars to investigate some possible alternative explanations. We have considered two possible dynamical effects which could produce significant $`R_\mathrm{c}`$ evolution: time-varying external tidal fields and a varying population of primordial binaries.

Throughout this paper we have analysed our simulations using the same procedures as were applied by Mackey & Gilmore (2003) to the observed data. We emphasised the importance of such an approach by comparing the evolution of the observational core radius to that of the half-mass radius and the density radius. As has been pointed out by a number of authors (e.g. Portegies Zwart et al., 2001), in general there are no simple relations between these quantities. Thus observations of simulated data sets are essential to ensure that comparisons between observed and simulated data are made in a consistent manner.

We considered clusters on circular and elliptical orbits of equal time-averaged radius about a point mass representing the LMC. The clusters are scaled so that their limiting radii correspond to the pericentric tidal radius of the clusters on elliptical orbits. However, at the start of each simulation the clusters lie well within the tidal radius which corresponds to their initial galactocentric position (apocentre for the clusters on elliptical orbits and the circular radius for the clusters on circular orbits). Initially, the clusters expand due to stellar mass loss but the competing effect of mass segregation means that $`R_\mathrm{c}`$ decreases slightly. At later times, two-body relaxation processes drive a further expansion of the clusters because the clusters are not tidally limited. This expansion is mirrored by a growth in $`R_\mathrm{c}`$. However, the rate of $`R_\mathrm{c}`$ expansion is found to be identical irrespective of the orbit on which the cluster is moving. We conclude that differences in the cluster orbits cannot be responsible for the observed differences in $`R_\mathrm{c}`$ among the LMC clusters.

Given that the binary fractions of all but the youngest clusters are currently unconstrained by observations, we compared the evolution of the $`R_\mathrm{c}`$ of a cluster containing no primordial binaries with that of a cluster with a large initial hard binary fraction. The presence of significant numbers of binaries leads to increased growth of $`R_\mathrm{c}`$ compared with the evolution of clusters containing only single stars. Binary fractions of $`50`$ per cent and $`25`$ per cent produce similar degrees of $`R_\mathrm{c}`$ expansion, while a $`10`$ per cent binary fraction results in noticeably less growth. Even in the most extreme case of a cluster with a binary fraction of $`50`$ per cent on an elliptical orbit, however, $`R_\mathrm{c}`$ increases by at most a factor of two. Not only is this smaller than the increase required to explain the LMC observations but the required spread in primordial binary fraction is already difficult to reconcile with the observed similarity of the cluster stellar luminosity functions (de Grijs et al., 2002b).

Our simulations have demonstrated that neither large variations in primordial binary fraction nor differences in cluster orbits about the LMC can account for the $`R_\mathrm{c}`$–age relation seen in the LMC clusters. The origin of the spread in $`R_\mathrm{c}`$ among the intermediate-age clusters therefore remains an open issue. We now return briefly to the observed data on the LMC cluster system as a whole to see whether there are any further clues to the origin of the trend.

Fig. 14 presents evidence for a systematic difference in the central density of the large-$`R_\mathrm{c}`$ clusters relative to the rest of the cluster population: these clusters have systematically lower central densities for a given cluster mass. This could be indicative of variations in star formation efficiency (SFE) between the clusters: low SFE leads to rapid cluster expansion and consequently a drop in the central density (e.g. Goodwin, 1997). Simulations of cluster expansion due to residual gas expulsion from clusters with low SFE (Goodwin, 1997) do not appear to show sustained core expansion of the magnitude seen in the intermediate-age LMC clusters. However, the core radii quoted in Goodwin (1997) are those of the mass distribution in the simulated clusters which, as was discussed in Section 3, may not accurately reflect the changes in the light distribution. We are currently investigating the long-term evolution of clusters which experience different degrees of primordial gas loss owing to SFEs in the range $`0.250.60`$, required to explain the range of core radii seen in the youngest LMC clusters (Goodwin, 1997). We are also simulating the effects of primordial binaries on larger-$`N`$ clusters. Of particular interest is the possibility that the heating effects of a large binary population could enable a cluster to retain a large initial $`R_\mathrm{c}`$ produced by gas expulsion.

## Acknowledgments

MIW acknowledges financial support from PPARC. JRH acknowledges financial support from a Kalbfleisch Fellowship and thanks the Institute of Astronomy (Cambridge, UK) for supporting a visit during this work. ADM would like to acknowledge the support of a Trinity College ERS grant and a British government ORS award. CAT thanks Churchill College for a fellowship. We are very grateful to Sverre Aarseth for assistance with the NBODY4 code and for numerous useful discussions. We thank Pete Bunclark at the Institute of Astronomy, Cambridge, for handling the installation and setup of the GRAPE-6 hardware. We also thank HongSheng Zhao for helpful discussions and Richard de Grijs and Douglas Heggie for useful discussions and comments on an earlier version of this paper.


text:
# Multiple Andreev Reflections in a Carbon Nanotube Quantum Dot

Present address : \]Cavendish Laboratory, Cambridge, UK

## Abstract

We report resonant multiple Andreev reflections in a multiwall carbon nanotube quantum dot coupled to superconducting leads. The position and magnitude of the subharmonic gap structure is found to depend strongly on the level positions of the single-electron states which are adjusted with a gate electrode. We discuss a theoretical model of the device and compare the calculated differential conductance with the experimental data. (pdf including figures, see: www.unibas.ch/phys-meso/Research/Papers/2003/MAR-MWNT.pdf)

Andreev reflection,quantum dots,carbon nanotubes,electric transport  preprint:

The electronic transport properties of quantum dots coupled to metallic leads have been the object of extensive theoretical and experimental study Kouwenhoven1 . When weakly coupled to its leads, the low-temperature transport characteristics are dominated by size and charge quantization effects, parameterized by the single-electron level spacing $`\mathrm{\Delta }E`$ and the single-electron charging energy $`U_C`$. When the coupling of the quantum dot to the source and drain electrodes is increased, higher-order tunneling processes such as the Kondo effect become important Kouwenhoven2 . New effects are expected when the leads coupled to the quantum dot are superconductors. In that case electron transport at small bias voltages is mediated by multiple Andreev reflection (MAR) Andreev ; Blonder . Unlike conventional S-N-S devices, however, the MAR structure is now expected Yeyati ; Johansson1  to strongly depend on the level positions of the single-electron states of the quantum dot which can be tuned with a gate electrode. The influence of Coulomb and Kondo correlations have been addressed theoretically in Refs. Zaikin .

Because MAR is suppressed rapidly for low-transparency junctions, its observation requires a relatively strong coupling between the leads and quantum dot. Even more so as on-site Coulomb repulsion, which is common to weakly coupled dots, is expected to reduce Andreev processes even further. A quantum dot very weakly coupled to superconducting leads has been studied experimentally by Ralph et al. Ralph . In this case the transport characteristics were indeed dominated by charging effects and MAR was completely suppressed.

The coupling to the leads, expressed in the life-time broadening $`\mathrm{\Gamma }`$ of the levels of the quantum dot, should be compared to the superconducting gap energy $`\mathrm{\Delta }`$. Favourable for the observation of MAR in a quantum dot are coupling strengths of order $`\mathrm{\Gamma }\mathrm{\Delta }`$ and a small charging energy $`U_C<\mathrm{\Delta }`$. Together with the restriction that $`\mathrm{\Gamma }<\mathrm{\Delta }E`$ (for any quantum dot) this leads to the approximate condition $`U_C\mathrm{\Gamma }\mathrm{\Delta }E`$. For most quantum dots typically the opposite is true and $`\mathrm{\Delta }EU_C`$. It has recently been shown Buitelaar1  however, that well-coupled multiwall carbon nanotube (MWNT) quantum dots can have favourable ratio’s of $`\mathrm{\Delta }E/U_C2`$ and $`U_C`$ can be as small as 0.4 meV, comparable to the energy gap $`2\mathrm{\Delta }`$ of a conventional superconductor like Al.

Here we report on the experimental study of resonant MAR in a MWNT quantum dot. The superconducting leads to the MWNT consist of an Au/Al bilayer (45/135 nm). Before investigating the system in the superconducting state, the sample is first characterized in the normal state by applying a small magnetic field. From these measurements relevant parameters such as $`\mathrm{\Delta }E`$, $`U_C`$, and $`\mathrm{\Gamma }`$ are obtained. We then discuss a theoretical model that describes the differential conductance of an individual level in a quantum dot coupled to superconducting electrodes. In the final part of the paper we compare the calculated differential conductance with the experimental data.

Figure 1 shows a greyscale representation of the differential conductance $`dI/dV_{sd}`$ versus source-drain ($`V_{sd}`$) and gate voltage ($`V_g`$) at $`T=50`$ mK when the contacts are driven normal by a small magnetic field. The dotted white lines outline the onset of first-order tunneling and appear when a discrete energy level of the quantum dot is at resonance with the electrochemical potential of one of the leads. ¿From these and other electron states measured for this sample (about $`20`$ in total), we obtain an average single-electron level spacing $`\mathrm{\Delta }E0.6`$ meV and a charging energy $`U_C0.4`$ meV. Since $`U_C=e^2/C_\mathrm{\Sigma }`$ this yields $`C_\mathrm{\Sigma }=400`$ aF for the total capacitance which is the sum of the gate capacitance ($`C_g`$) and the contact capacitances $`C_s`$ (source) and $`C_d`$ (drain). From the data of Fig. 1a we obtain $`C_g/C_\mathrm{\Sigma }=0.0036`$, and $`C_s/C_d=0.45`$. The lifetime broadening $`\mathrm{\Gamma }`$ is obtained from the width of the single-electron peaks at finite source-drain voltage (taking the background conductance into account) and is found to be $`\mathrm{\Gamma }0.35`$ meV.

The high-conductance ridge around $`V_{sd}=0`$ mV in Fig. 1a is a manifestation of the spin-$`1/2`$ Kondo effect occurring when the number of electrons on the dot is odd. As a result, the Coulomb valley in the conductance has disappeared in this region of $`V_g`$ and at $`50`$ mK, which is the base temperature of our dilution refrigerator, only a single peak remains, see Fig. 1b. This is known as the unitary limit of conductance. The appearance of the Kondo effect is an indication that the coupling to the leads is relatively strong. We will not discuss the Kondo effect here, and instead refer to Ref. Buitelaar2 .

When the magnetic field is switched off the leads become superconducting. To calculate the expected differential conductance in the superconducting state of the leads we have used the non-equilibrium Green-function technique Greenfunction . We model the quantum dot as a series of spin-degenerate resonant levels coupled to superconducting electrodes, which are assumed to have a BCS spectral density. Note, that neither electron-electron interaction (Coulomb blockade) nor exchange correlations (Kondo effect) are accounted for in the model, which may, therefore, not explain all details of the actual measurements. However, the interplay between MAR and resonant scattering already leads to strongly nonlinear IV-characteristics and reproduces some of the key features of the data. The details of the model calculation will be presented elsewhere Bruderbelzig . The main parameters entering the calculation are the two tunneling rates $`\mathrm{\Gamma }_{s(d)}=2\pi |t_{s(d)}|^2N_0`$, related to microscopic tunneling amplitudes $`t_{s(d)}`$ via the density of states $`N_0`$ of the source (drain) contact. In the model we account for the gate voltage by a shift of the level, which can be adjusted according to the experimentally observed Coulomb diamonds, see Fig. 1.

The discrete nature of the single-electron states is most pronounced when $`\mathrm{\Gamma }`$ is small. Therefore, before presenting the model calculation that directly compares to the experimental data, we first discuss the transport characteristics of a single spin-degenerate level with a relatively weak and symmetric coupling to the superconducting leads along the lines of Refs. Yeyati ; Johansson1 . The total tunneling rate $`\mathrm{\Gamma }\mathrm{\Gamma }_s+\mathrm{\Gamma }_d`$ is set to $`\mathrm{\Delta }`$. Figure 2 shows the corresponding greyscale representation of the calculated differential conductance $`dI/dV_{sd}`$ versus $`\phi _g:=eV_gC_g/C_\mathrm{\Sigma }`$ and $`eV_{sd}`$. The peak structure in $`dI/dV_{sd}`$ at $`V_{sd}<2\mathrm{\Delta }/e`$ is the result of MAR. In general, Andreev channels become available for transport at voltages $`V_{sd}=2\mathrm{\Delta }/ne`$, where $`n`$ is an integer number. These positions are indicated by the horizontal dashed black lines in Fig. 2a. The appearance and magnitude of the MAR peaks, however, is strongly dependent on the position of the resonant level in the quantum dot with respect to the Fermi energy of the leads. Only those MAR trajectories that connect the resonant level to the leads’ BCS spectral densities give a significant contribution to the current. Consider, for example, the position marked by $`C`$ in Fig. 2a, which corresponds to the schematics of Fig. 2c and indicates the position $`(\phi _g,eV_{sd})=(0,2\mathrm{\Delta }/3)`$. The corresponding second-order Andreev trajectory connects the gap edges of the source and drain electrodes and includes the resonant level which is situated exactly in between the respective Fermi energies. This results in the large peak in $`dI/dV_{sd}`$ seen in Fig. 2b.

A similar peak is absent at $`(\phi _g,eV_{sd})=(0,\mathrm{\Delta })`$, corresponding to point $`D`$ in Fig. 2a. At this voltage, first-order Andreev reflection becomes possible. The corresponding trajectories (see Fig. 2d), however, do not directly connect the resonant level to the leads spectral densities, and therefore do not significantly contribute to the current. Only when the gate voltage is adjusted to align the level with the Fermi energy of one of the leads (indicated by the arrows in Fig. 2a) a peak in $`dI/dV_{sd}`$ is observed. It can be shown (for symmetric junctions) that the subharmonic gap structure at $`V_g=0`$ is suppressed for all voltages $`V_{sd}=2\mathrm{\Delta }/ne`$ with $`n=\text{even}`$ Yeyati ; Johansson1 . When $`V_{sd}`$ is increased beyond $`\mathrm{\Delta }/e`$, peaks are observed either when the level stays aligned with the electrochemical potential of the leads (red dashed lines: $`\phi _g=\pm eV_{sd}/2`$) or when the level follows the gap edges as an initial or final state of an Andreev process (blue dash-dotted lines: $`\phi _g=\pm (\mathrm{\Delta }eV_{sd}/2)`$ or $`\phi _g=\pm (\mathrm{\Delta }3eV_{sd}/2)`$).

We now turn to the actual measurements of the differential conductance when the leads are in the superconducting state. Figure 3 shows a greyscale representation of the measured $`dI/dV_{sd}`$ versus $`V_{sd}`$ and $`V_g`$ at $`B=0`$ mT for the same single-electron state of Fig. 1. A number of differences between the normal state (Fig. 1) and superconducting state (Fig. 3) can be observed. The horizontal high-conductance lines at $`V_{sd}=\pm 0.2`$ mV in Fig. 3, for example, are attributed to the onset of quasi-particle tunneling when $`V_{sd}=2\mathrm{\Delta }/e`$. The subgap structure at $`V_{sd}<2\mathrm{\Delta }/e`$ (i.e. below $`0.2`$ mV) is attributed to MAR. As anticipated, the magnitude and the position (dashed white lines) of MAR peaks depend on $`V_g`$. To allow for comparison with theory, the adjustable parameters of the model are set to the values obtained from the measurement of Fig. 1. The most important parameter is the coupling between the electrodes and the dot which turns out to be $`\mathrm{\Gamma }3.5`$$`\mathrm{\Delta }`$. The voltage division between the two tunnel barriers separating the quantum dot from the leads is $`C_s/C_d=0.45`$. The individual tunneling rates are not exactly known but are not expected to show a strong asymmetry Kondo . The neighboring single-electron states, separated by $`\mathrm{\Delta }E6.5\mathrm{\Delta }`$, are included in the calculation. The finite temperature of the experiment is taken into account and set to $`T=0.1`$$`\mathrm{\Delta }`$.

The resulting calculated greyscale representation of the differential conductance is shown in Fig. 4. The overall appearance clearly resembles the measured data of Fig. 3. For example, both the model and the measured data show a large peak in $`dI/dV_{sd}`$ around $`V_{sd}=0`$ mV when the electron state is at resonance with the electrochemical potential of the leads (i.e. at $`V_g=0`$). When the level is moved away from this position, the linear-response conductance rapidly decays to values below its normal-state value. In contrast, the differential conductance peak at $`V_{sd}=2\mathrm{\Delta }/e`$ shows the opposite behavior (both in the model as in the experiment). At $`V_g=0`$, this peak is much less pronounced then at lower values of $`V_g`$. This can also be seen in the $`dI/dV_{sd}`$ line traces of Fig. 3b-c and Fig. 4b-c (solid curves). These observations are similar to conventional S-N-S structures, such as atomic-sized break junctions Scheer1 : for large transparencies of the junction a peak is observed at $`V_{sd}=0`$ but no structure at $`2\mathrm{\Delta }`$ while for small transparencies a gap is observed around $`V_{sd}=0`$ and a large peak at $`2\mathrm{\Delta }`$ marks the onset of quasi-particle tunneling. For these systems, the transparency depends on the atomic arrangement of the junction, here the effective transparency can be tuned by moving the level position of a single-electron state with a gate electrode. The effective transparency is large if the level is aligned with the Fermi levels of the leads (on resonance) and it is small otherwise (off resonance).

The subharmonic gap structure is clearly visible in the measured data of Fig. 3 and has a similar gate-voltage dependence as in the model calculation of Fig. 4. However, there are several differences. The most dramatic one is the pronounced peak at $`(V_g,V_{sd})=(0,\mathrm{\Delta }/e)`$ in the measurement (Fig. 3c). Because this position corresponds to an even MAR cycle it should be absent based on our previous consideration (see Fig. 2d).

Let us compare theory and experiment by focussing onto the $`dI/dV_{sd}`$ lines traces shown in Fig. 3b-c and Fig. 4b-c (solid curves). In panels b) the dot level is off resonance, while it is on resonance in panels c). For the former case, experiment and theory agree fairly well. The differences in MAR structure between model and experiment are much more pronounced at the resonance position. Whereas the experiment (Fig. 3c) reveals pronounced peaks at $`\mathrm{\Delta }/2`$, $`\mathrm{\Delta }`$ and $`2\mathrm{\Delta }`$, the calculated $`dI/dV_{sd}`$ (Fig. 4c, solid) reveals fine structure for small $`V_{sd}`$ and pronounced peaks at $`2\mathrm{\Delta }/3`$ and $`2\mathrm{\Delta }`$. According to our previous discussion $`dI/dV_{sd}`$ should indeed show a pronounced peak at $`V_{sd}=2\mathrm{\Delta }/3e`$, if the dot level is centered in the middle, i.e. for $`\phi _g=0`$ at point C in Fig. 2a and Fig. 2c. It rather appears in the experiment that, contrary to expectations, the subgap feature at $`2\mathrm{\Delta }/3`$ is missing, while the ‘forbidden’ at $`\mathrm{\Delta }`$ is present. Such behavior would only be expected for very asymmetric junctions having $`C_s/C_d1`$ Johansson2 , which is not the case in the present work.

There are different imaginable scenarios that may account for the observed $`\mathrm{\Delta }`$ peak and the lack of fine structure around $`V_{sd}=0`$ mV in the data of Fig. 3c. Inelastic scattering processes inside the dot, for example, would broaden and obscure higher-order MAR features. Other possible reasons may be found in a broadened BCS spectral density (the superconductor consists of a bilayer of Au/Al Scheer2 ) or a suppression of higher-order MAR due to the on-site Coulomb repulsion.

In a phenomenological approach, we may try to account for the additional broadening by manually introducing larger bare couplings $`\mathrm{\Gamma }_{s,d}`$. Many curves with varying parameters were calculated of which a representative set is displayed in Fig. 4b-c (dashed curves) corresponding to relatively large dot-electrode couplings of $`\mathrm{\Gamma }_s=2.5`$$`\mathrm{\Delta }`$ and $`\mathrm{\Gamma }_d=3.5`$$`\mathrm{\Delta }`$. For the off-resonance position (Fig. 4b) the main effect of the larger $`\mathrm{\Gamma }`$ is the increased magnitude of $`dI/dV_{sd}`$. In contrast, the MAR structures significantly changes for the resonance position (Fig. 4c). Remarkably, at large coupling $`\mathrm{\Gamma }`$, peaks now appear at $`2\mathrm{\Delta }`$, $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }/2`$. These peaks do not originate from the resonant level, but from the two neighboring ones which are off resonance (the dot levels are spaced by $`\mathrm{\Delta }E6.5`$$`\mathrm{\Delta }`$). Though the agreement is now reasonable, there is one remaining problem. We were unable to reproduce the relative peak height between the $`2\mathrm{\Delta }`$ and $`\mathrm{\Delta }`$ peaks. Using any reasonable set of parameters, the $`2\mathrm{\Delta }`$ peaks is always larger than the $`\mathrm{\Delta }`$ peak in the model, while it is the opposite in the experiment. We emphasize that the model does not take into account interaction and correlations. Since a Kondo resonance is observed in the normal state, which need not be suppressed in the superconducting state Buitelaar2 , this may be the origin of the discrepancy. The Kondo resonance changes the spectral density in the leads by adding spectral weight to the center of the gap and removing spectral weight from the gap edges. The former tends to enhance the $`\mathrm{\Delta }`$ peak, while the latter tends to suppress the $`2\mathrm{\Delta }`$ one. This explanation is attractive, but more work both in theory and experiment is needed to substantiate it.

In conclusion, we have investigated the non-linear conductance characteristics of a quantum dot coupled to superconducting electrodes. We find a strong dependence of the MAR structure on the level position of the single-electron states. The experimental data is compared with a theoretical model, assuming a BCS density-of-states in the electrodes and an interaction-free dot. Reasonable agreement is possible, if the tunneling coupling to the leads is enhanced by a factor $`2`$ in the model as compared to the experimental value. There are additional subtle differences which point to the importance of interaction and exchange correlations.

###### Acknowledgements.

We acknowledge contributions by H. Scharf and M. Iqbal. We thank L. Forró for the MWNT material and J. Gobrecht for providing the oxidized Si substrates. This work has been supported by the Swiss NFS, BBW, and the NCCR on Nanoscience.