Manifestation of dynamic Jahn–Teller distortions and surface interactions in scanning tunnelling microscopy images of the fullerene anion C−60

Using scanning tunnelling microscopy (STM), it is possible to observe detailed structure of the molecular orbitals (MOs) of fullerene anions C−60. However, understanding the experimental observations is not straightforward because of the inherent presence of Jahn–Teller (JT) interactions, which (in general) split the MOs in one of a number of equivalent ways. Tunnelling between equivalent distortions means that any observed STM image will be a superposition of images arising from the individual configurations. Interactions with the surface substrate must also be taken into account. We will show how simple ideas involving a symmetry analysis and Hückel molecular orbital theory can be used to understand observed STM images without need for the more usual but more complicated density functional calculations. In particular, we will show that when the fullerene ion is adsorbed with a pentagon, hexagon or double-bond facing the surface, STM images involving the lowest unoccupied molecular orbital (LUMO) can be reproduced by adding together just two images of squares of components of the LUMO, in ratios that depend on the strength of the JT effect and the surface interaction. It should always be possible to find qualitative matches to observed images involving any of these orientations by simply looking at images of the components, without doing any detailed calculations. A comparison with published images indicates that the JT effect in the C−60 ion favours D3d distortions.


Introduction
Scanning tunnelling microscopy (STM) can be used to directly image molecules experiencing distortion due to the Jahn-Teller (JT) effect, and hence to subsequently gather information on the nature of that effect that is difficult to determine by other means. Ions of the fullerene molecule C 60 are ideal candidates in which to observe JT effects in STM: they are large enough to be easily discernible via STM, and their electronic and vibrational states are highly degenerate. The JT effect will lift the electronic degeneracy, and the signature of this will be apparent in observed STM images. If the JT effect is very strong, an isolated C − 60 anion could become locked into a distorted structure of lower symmetry than the icosahedral symmetry of the original molecule. The effect of this on the lowest unoccupied molecular orbital (LUMO) was considered in [1]. However, there will be a set of isoenergetic distortions, referred to as 'orientomers' in [1], which look identical to each other but have different relative orientations. For realistic coupling strengths, tunnelling between these orientomers is expected to take place on a timescale that is fast compared to the millisecond timescale of the STM technique, even at low temperatures [2]. This dynamic JT effect must be taken into account when simulating observed STM images.
When the adiabatic potential energy surface (APES) has large barriers between equivalent distortions, the system can be expected to 'hop' between the configurations and spend negligible time at intermediate distortions. Observed STM images will then be a superposition of images due to the individual orientomers. However, when the barriers are negligibly small, the distortion will rotate from one orientomer to another, in a process known as pseudorotation. The motion for intermediate barrier heights is called a hindered pseudorotation. Observed STM images for both free and hindered pseudorotation will be a temporal average containing contributions 3 from intermediate distortions. Both hopping and pseudorotation will be considered in this paper. This dynamic behaviour is different to the static JT effect considered in [1], where the system can't tunnel between wells so observed STM images will be images due to a single well.
Simulations of images of the neutral C 60 molecule have been published that were obtained using density functional theory (DFT) [3,4] and also using a Hückel molecular orbital (HMO) approach [5]. The agreement between the simulated images and experimental observations is generally good. In [1], the HMO approach was used to describe the LUMO of C − 60 ions in which static JT effects were taken to be the dominant interaction and surface interactions were included as a perturbation. This gave potential matches to the experimental data for C 60 ions on surfaces in which there is a buffer layer between the molecules and the surface in order to reduce the strength of the surface interaction. We will use a similar approach in this paper, but we will take dynamic effects into account and also treat the JT and surface interaction on an equal footing. We will also show how the general features of the results occur due to symmetry considerations in general, irrespective of the model used to describe the surface interactions.
We use the HMO approach in this paper rather than the more usual DFT because, due to its computational simplicity, it can be readily used to explore a whole landscape of different JT and surface interaction parameters and to see the effect of these parameters on the simulated results. While there is no doubt that DFT simulations generally give a good account of experimental observations, it is much more time-consuming to implement and hence the roles of different influencing factors are much harder to extract. Also, when symmetry arguments are used to find the molecular orbitals (MOs), there are no additional complications when dealing with degenerate orbitals. This is not the case for DFT, where problems can easily occur when electronic degeneracies are present [6].
Theoretical and experimental evidence indicates that due to the JT effect alone, C − 60 ions are most likely to distort to D 3d symmetry [7][8][9]. However, this evidence is not conclusive, with the energies of D 5d wells expected to be very similar to those of the D 3d wells [10]. JT theory indicates that distortions to both symmetries are possible, depending upon the values of the quadratic JT coupling constants [11,12] (unless there is no quadratic coupling, in which case there is a continuum of minimum energy points [13][14][15]). D 2h points could also be possible, but only for certain large negative values of the coupling constants [1]. Estimates of the linear JT coupling constants are known, but theoretical values are generally significantly smaller than experimentally determined values. One possible reason for this is the neglect of the contribution made by quadratic factors [9]. In any case, values for the quadratic coupling constants are not known. Therefore, in this paper, we will consider JT distortions to D 3d or D 5d symmetry in general terms. It should be noted that while we will assume that the underlying JT effect will result in distortions to D 3d or D 5d symmetry, interactions with a surface substrate will reduce the symmetry further.
In section 2, we will formulate both the JT problem and the effect of surface interactions. In section 3, we will then show how our formalism can be used to predict STM images of C − 60 ions adsorbed on a surface substrate with either a pentagon, a hexagon or a double-bond facing the surface. We will show how the results can be interpreted in terms of changes in the electronic state associated with each orientomer as a function of the strength of the surface interaction. We will also show that images taking into account either hopping or pseudorotation between orientomers for adsorption in these orientations can always be interpreted in terms of a superposition of either two or three images of components of the LUMO. In section 4, we will show how our results can be extended to incorporate interactions between neighbouring C 60 4 ions in a hexagonal lattice, and how the results can be used to obtain an alternative explanation of published STM images. This in turn gives useful information on the nature of the JT effect. Finally, in section 5, we give some general discussions regarding our results and draw some general conclusions.

Jahn-Teller (JT) and surface interactions
In treating an isolated C − 60 ion adsorbed onto a surface substrate, we must consider the JT effect and the interaction with the surface. We will now give details of how the JT interactions can be formulated and then how surface interactions can be incorporated.

JT interactions
The LUMO of a C − 60 ion is subject to a T 1u ⊗ 8h g JT effect, where a T 1u electronic state couples to eight different sets of h g vibrational levels. However, we will only consider the splitting of the electronic state and ignore the associated displacement of the nuclear coordinates because these displacements are expected to be far too small to be detectable at the resolution currently available using STM. The results will not depend upon the values of the JT coupling constants used to describe the h g mode. It only matters whether they are in the appropriate range to favour D 5d or D 3d distortions. Therefore there will be no difference between treating a single h g mode and the multimode problem, so we will only consider a single mode. The JT interactions can be expressed in terms of a linear coupling Hamiltonian and two linearly-independent quadratic coupling Hamiltonians [13,15]. The two quadratic Hamiltonians arise because the product H ⊗ H contains 2H so is non-simply reducible. The JT Hamiltonian is traditionally defined in terms of a coordinate system that has twofold axes through the centre of a double bond between two hexagons [15]. In this paper, we will call these molecule-fixed axes {X, Y, Z }.
The potential part of the JT Hamiltonian can be written in the form where V 1 is a linear coupling constant, V 2 and V 3 are the quadratic coupling constants, µ is the mass and ω is the frequency of the h g mode. H 1 depends linearly on the five vibrational coordinates Q γ representing the h g mode, and H 2 and H 3 depend upon quadratic powers. We write the Q γ more explicitly as {Q θ , Q , Q 4 , Q 5 , Q 6 }. In terms of d-orbitals, these transform Explicit matrix forms for each of the three Hamiltonians in equation (1) can be found in [13].
For given values of the linear and quadratic coupling constants, there will be certain values of the Q γ s that minimize the energy of the lowest APES to give minima, which are usually of D 5d or D 3d symmetry [13]. There are six D 5d minima which, following [13], we label A to F, and ten D 3d minima which we label a to j. There are also 15 D 2h points at the centres of the C-C double bonds, which we will label A-O as in [16].
The five Q γ can, in general, be written in terms of one distance and four angles [15]. However, on the minimum-energy surface (for both the D 5d and D 3d cases), they can be written in terms of just two angles, θ and φ, having the usual definitions for spherical polar coordinates Pictorial representations of paths between (a) D 5d wells, (b) nearestneighbour D 3d wells and (c) next-nearest neighbour D 3d wells. Non-overbarred letters correspond to the wells in table II of [1] and overbarred ones to diametrically-opposite points.
(the inclination and azimuth), such that These can be obtained from [15] after changing to our notation and converting to our definitions of Q θ and Q . Hence a visual representation of the positions of the minima can be obtained by plotting points {sin θ cos φ, sin θsin φ, cos φ} on the surface of a sphere [13]. It is also useful to include a diagrammatic representation of the structure of the APES by varying the radial coordinate at positions in between the wells in a manner that represents the relative energies of these positions. The result is shown in figure 1. The diametrically-opposite points are also shown, using overbarred labels. The images in figure 1 also show a representation of shortest paths on the APES when moving from one well to another due to hindered pseudorotation between equivalent wells. For D 5d wells, all minima have equal separations, as shown in figure 1(a). However, for D 3d symmetry each minimum has three nearest-neighbours and six next-nearest neighbours. Shortest paths between nearest neighbours are shown in figure 1(b), and shortest paths between next-nearest neighbours are shown in figure 1(c). Later, we will consider pseudorotation along these paths. While a complete calculation would allow for all possible configurations, this will 6 allow us to estimate the effect on STM images of including intermediate configurations due to hindered pseudorotation. Intuitively, it would seem likely that the minimum-energy paths between next-nearest neighbour wells, such as between a and c, will be somewhere in between the shortest paths drawn on figure 1(c) and the path via a nearest-neighbour well, such as from a to c via g. However, we will later find that detailed knowledge of the pseudorotational dynamics is not required.
The electronic T 1u states characterizing each well can also be written in the form c x ψ x + c y ψ y + c z ψ z , where {c x , c y , c z } are constants and {ψ x , ψ y , ψ z } are components of the T 1u state transforming as {T 1ux , T 1uy , T 1uz } respectively (for any chosen Cartesian coordinate system). These states have the same transformation properties as p-orbitals aligned along {x, y, z}. It is found that in the coordinate system {X, Y, Z }, the points {c x , c y , c z } can be represented using the same polar angles θ and φ as used to represent the positions of the wells. Hence figure 1 can be seen as representing either the positions of the wells or the electronic states. When representing electronic states, the overbarred labels correspond to the same electronic states as the unbarred labels but defined with the opposite sign. A pictorial representation of the electronic states will be useful when we interpret our predicted STM images in section 3, although in that section we will represent complete points and not the angular variation with the radial coordinate representing energy. As the states are normalized, the points will therefore lie on the surface of a sphere.

Surface interactions in general terms
In this paper, we will assume that the C 60 ions are adsorbed onto a surface substrate such that the surface perturbs the positions of the atoms of the C 60 molecule closest to it in some way so that its symmetry is reduced from the icosahedral symmetry (I h ) of an isolated molecule by removing the equivalence between an axis perpendicular to the surface and those in the plane of the surface. In this way, we do not need to consider the structure of the surface in detail. Different surfaces will result in different magnitudes for the surface interaction, but if the effect on the symmetry of the molecule is the same then it will lift the degeneracy of the MOs in the same way. This approximation can be expected to be valid for 'simple' surfaces such as Ag or Au, but may not be appropriate for more complex surfaces such as Si(111)-(7 × 7) or WO 2 /W(110), although even in these cases there may be some adsorption sites for which this method is valid. Our approach is similar to that taken by Pascual et al [17], who introduced the symmetry-breaking effect of a surface into their DFT simulations by applying a uniaxial compression in the direction perpendicular to the substrate.
From the above symmetry considerations and the pictorial representation in figure 1 of the electronic states/representation of the positions of D 3d and D 5d wells in coordinate space, it is a simple matter to determine which wells will remain equivalent when a surface interaction is added. However, we first need to define an appropriate coordinate system {x, y, z} for specifying the surface interaction. In general, this is different to the molecule-fixed system {X, Y, Z } used for the JT interaction. In common with [1] and [5], we will define our z-axis to be perpendicular to the surface and through the centre of the C 60 molecule. We take our y-axis to be through the centre of a double bond between two hexagons, which means that y always corresponds to the molecule-fixed axis Y . In this way, the three highest-symmetry orientations can be obtained by rotating the molecule about the y-axis in the X -Z plane, as shown in figure 1 of [1]. Table 1. Groupings due to symmetry considerations of equivalent D 5d , D 3d and D 2h points for pentagon, hexagon and double-bond prone orientations. Groups on the same row between braces are related to each other by reflections in the x-z plane. The groups that can be global minima in our model when 1 > 0 are marked with a + and when 1 < 0 with a * , assuming that 2 = 0 for the double-bond prone orientation (where the 's are defined in section 2.3). The three orientations which result in the highest symmetries upon surface adsorption involve either a pentagon, a hexagon or a double-bond prone to the surface. We will consider these three orientations in detail, although it is fairly straightforward to extend the method to other symmetries, such as when the molecule is tilted away from one of these ideal orientations or when a single bond is prone to the surface, as we will see in section 4. The results of the symmetry analysis for the high-symmetry orientations are collected in table 1. The z-components of the electronic state for all points in a group are the same, and the results for all members in a group are related to each other by rotations about the z-axis by 2π/5, 2π/3 and π respectively. A similar analysis can be performed for the D 2h points using the pictorial representation of their positions in figure 1 of [16]. The results are also shown in table 1. As the timescale of an STM experiment is expected to be orders of magnitude larger than the time taken to tunnel between equivalent configurations [2], all equivalent (isoenergetic) distortions in a set will be imaged together in an STM experiment. In the table, groups on the same line separated by braces are related to each other by reflections in the x-z plane. Therefore, by symmetry both sets can also be expected to be imaged together. This will be confirmed below by our calculations using our specific model. As mentioned previously, when the surface interaction is included, the minimum-energy distortions correspond to lower symmetries than D 5d and D 3d ; our labels indicate the symmetry that the underlying JT effect would produce in the absence of the surface interaction.
If the surface interaction can be treated as very weak, but sufficiently strong to remove the equivalence between all JT wells, predictions of what would be observed in an STM image could be obtained by combining the images of different wells given in [1]. However, in general 8 the surface interaction will both alter the electronic state and change the values of the normal mode coordinates Q γ at which the minima occur. We then need to determine the electronic state for each well in terms of JT and surface interaction parameters. To do this, it is necessary to formulate a Hamiltonian to represent the surface interaction and find new solutions to the combined JT and surface interaction Hamiltonians. We can write down a form for the surface interaction Hamiltonian by considering the reduction in symmetry upon adsorption on a surface. Only certain rotation and reflection operations associated with the axis perpendicular to the surface will survive as symmetry operations when the symmetry is reduced from I h . The actual operations that survive depend on the orientation of the C 60 molecule. The pentagon, hexagon and double-bond prone orientations result in symmetries of C 5v , C 3v and C 2v respectively.
We can use group theory to show that the groups of wells in table 1 will remain equivalent under the action of a surface interaction. Consider one well with electronic state c x ψ x + c y ψ y + c z ψ z . To determine the states for equivalent wells in the C 3v symmetry that applies to the hexagon-prone orientation, for example, we can apply the C 3v group operations to this state. Applying the C 3 rotations yields states c Applying the σ v reflections doesn't generate any additional states. The three states are unique unless c x = c y = 0, in which case they represent the same point. This is consistent with table 1, where states are all in groups of three for this orientation, except for D 3d well c which is a unique point. For the double-bond prone orientation, applying the rotation and two reflections of the C 2v group to a state c x ψ x + c y ψ y + c z ψ z gives four different equivalent-energy states when a x , a y and a z are all non-zero, or two equivalent states when one of the coefficients is zero. This is again consistent with table 1. Similar arguments apply to the pentagon-prone orientation.
This grouping of equivalent wells has fairly fundamental consequences for the appearance of STM images. According to Tersoff-Hamann theory [18], The current I recorded in an STM experiment will be proportional to the sum of squares of degenerate wavefunctions ψ i at the Fermi energy, which is a sum over all equivalent minimum-energy states in our problem. For ease of comparison, we divide this sum by the number of wells so that the result is 'normalized'. This means that we define the current to be Due to the symmetric distributions of the wells, all cross-terms (involving ψ x ψ y etc) cancel in the expression for I for all of the three symmetric orientations we are considering. For example, with the general expressions for three equivalent wells in the hexagon-prone orientation obtained above, the current is . I can therefore be divided into contributions involving squares of the three components of the LUMO, ψ 2 x , ψ 2 y and ψ 2 z , only. In other words, we can write where the coefficients a x , a y and a z are constants for a given set of JT and surface interaction parameters. Furthermore, for the pentagon and hexagon prone orientations we find that the ψ x and ψ y states contribute equally, i.e. a x = a y (as in the example above). This means that the appearance of an STM image can be estimated, at least qualitatively, by examining images of ψ 2 x + ψ 2 y and ψ 2 z alone for the pentagon and hexagon-prone orientations, and images of ψ 2 x , ψ 2 y and ψ 2 z for the double-bond prone orientation. This conclusion does not relate to any specific model used for the surface interaction, arising simply as a result of the symmetry the interaction results in.

Model Hamiltonian for surface interactions
Having determined the reduced symmetry of the molecule when both JT and surface interactions are included, we can now use character tables to determine how a given MO will split due to surface interactions alone, and use the result to construct an appropriate form for the surface interaction Hamiltonian.
With the above definitions, we find that for all three orientations of interest, the ψ z component of the LUMO becomes a singlet transforming as A 1 . For the pentagon-prone case, the ψ x and ψ y components form a doublet transforming as E 1 , and for the hexagon-prone case they form a doublet transforming as E. This separation into a singlet and a doublet seems reasonable because the z direction (perpendicular to the surface) is being treated as different to x and y, whereas x and y are treated in the same way. For the double-bond prone case, ψ x is a singlet B 1 and ψ y is a singlet B 2 . We cannot have a doublet in this case because C 2v symmetry only supports singlet representations.
From the arguments above, an appropriate Hamiltonian for the surface interaction in the three orientations we consider is where 1 and 2 are dimensionless constants determining the strength of the surface interaction, and which can be positive or negative. 2 is zero for the pentagon and hexagonprone cases.
We have written the surface interaction Hamiltonian H sz in terms of basis states relating to the axes {x, y, z} and the JT Hamiltonian H JT in terms of basis states relating to {X, Y, Z }. These bases only coincide for the double-bond prone case. For the pentagon and hexagon-prone cases, we therefore need to either convert H JT to the basis relating to {x, y, z} or convert H sz to the basis relating to {X, Y, Z }. As H JT is the most complex Hamiltonian, we find it easier to convert H sz , although obviously the results will be the same whichever approach we use.
To convert between the bases, we need the matrix to rotate by an angle β in the x-z plane. Therefore the required surface interaction Hamiltonian in the rotated {X, Y, Z } basis is For the pentagon-prone and hexagon-prone cases, The approach we have taken here is different to that used in [1,5], which treated the C 60 molecule as sitting at the centre of a hexagon in the surface, at a site of C 6v symmetry and determined degeneracies by looking for combinations of states which have the same 'characters'. However, as C 6v is not a subgroup of I h , the results are not true characters, and depend upon the choice of basis set chosen. Furthermore, an axis in the plane of the surface is singled out as being the axis relevant for allowed symmetry operations, whereas we would expect the axis perpendicular to the surface to be more important. Also, due to the large size of the C 60 molecule compared to typical lattice spacings, we would not expect every molecule to be able to sit at the centre of a hexagon. In fact, the pattern of splitting in our current approach and that used previously are the same for the pentagon and hexagon-prone cases. In the double-bond prone case, we predict three singlets rather than the doublet and a singlet predicted previously. The form we have presented here also differs from [1] in that it doesn't conserve the centre of energy, although it is a trivial matter to preserve the centre of energy if required.

Pseudorotation and hopping
In the absence of any surface interaction or quadratic JT couplings, the lowest-energy APES is a multi-dimensional trough containing a continuum of equivalent-energy points. The system can move freely between all such points in a free pseudorotation. An STM image from this system would include contributions from all points on the trough equally. When quadratic JT coupling is introduced, the system undergoes a hindered pseudorotation in which the system can be expected to spend more time near to the positions of the wells than at the points in between. If the quadratic coupling is sufficiently strong, the system can be assumed to hop between the wells; intermediate points between the wells become irrelevant and any STM image will be a superposition of images associated with the wells alone. When a surface interaction is also introduced, certain wells become favoured. Different wells will be favoured for different JT coupling constants and different orientations of the molecule on the surface. If the quadratic coupling is weak, there will be a hindered pseudorotation that favours the minimum-energy wells, and if it is strong the system will hop between these wells.
To investigate the differences in the effect of pseudorotation and hopping on observed STM images, we will compare the results of hopping and pseudorotating between two minimumenergy configurations when a very weak surface interaction is present. For the pseudorotation, we will make the assumption that the system will follow the shortest-distance path between wells, and we will treat all points on this path equally. The differences between the images are likely to be greatest under a weak surface interaction because, as we will see in section 3, a strong surface interaction will tend to move the positions of the wells closer together so that less intermediate configurations will be sampled. It also reduces the differences between the electronic states. The dynamics of the real C 60 molecule is expected to be somewhere between our two limits.
We can determine the shortest-distance path between two wells in either electronic or coordinate space because, as mentioned in section 2.1, paths in both spaces can be represented in terms of variation of two angles-the spherical polar angles θ and φ. If the electronic state was not a triplet, e.g. for anions C n− 60 with 1 < n < 5 [15], this would not be possible and it would be necessary to perform the calculation in coordinate space.
As we will assume that the surface interaction is weak, we will take the start and end points on the pseudorotational paths to be at the angles obtained in the absence of any surface interaction. From [13], we can see that the polar angles (θ, φ) for the six D 5d wells A to F are {( π 2 − , π 2 ), ( π 2 + , π 2 ), ( , 0), (− , 0), ( π 2 , ), ( π 2 , − )} respectively, where The shortest path between wells A ↔ B, C ↔ D and E ↔ F (which represents the change in the electronic state when moving between these pairs of wells) can therefore immediately seen to be {( π 2 + a , π 2 ), (a , 0), ( π 2 , a )} respectively, where a varies between −1 and +1. The paths joining all other pairs of D 5d wells are not so obvious as they involve varying both θ and φ. However, applying projection operators to any of the three paths above generates all other paths. Our paths are of C 2h symmetry, passing through a D 2h point midway between the two wells (a = 0). The most significant consequence of this is that the two JT excited states can no longer be degenerate.
We can now repeat the same procedure for paths joining D 3d wells. From [13], we can see that the polar angles for the ten wells a to j are . The shortest path between the nearest-neighbour wells a ↔ b, c ↔ d and e ↔ f can therefore immediately seen to be {(a , π 2 ), ( π 2 − a , 0), ( π 2 , −a )} respectively, where a varies between −1 and +1 again. As before, the paths are of C 2h symmetry passing through a D 2h point midway between the two wells, which again means that the JT excited states cannot be degenerate. Paths between all other nearest neighbours can again be found by the application of projection operators.
If we choose to fix the path in electronic space, we then need to find the minimum possible energy at any point on the path by varying the Qs. Fixing the path in coordinate space yields restrictions on the relative values of the Qs. We then calculate the minimum energy by varying the Qs subject to these constraints. We find that the same results are obtained for both approaches.

Solutions with combined surface and JT interactions
We will start our analysis by determining the minimum-energy points favoured for a given set of JT coupling constants and surface interaction parameters by finding values for the Q γ that minimize the lowest-energy eigenvalue of H JT + H s . This extends the approach used in [1], where H s was treated as a perturbation. We will use the same labels to denote the minimumenergy configurations as those used in figure 1 to identify the positions of the D 5d and D 3d wells when there are only JT interactions present. The labels are such that the configurations evolve to the JT wells when the strength of the surface interaction is reduced to zero, with the wells labelled as in [1]. Results when the parameter 2 (determining the difference in energy between the x and y states) is set to zero are given in table 1 for the pentagon-prone, hexagon-prone and double-bond prone orientations.
Most of the results we will present are for the two cases are dimensionless quadratic coupling constants. These values of V 2 and V 3 are parameters preferring D 5d and D 3d distortions respectively, as shown in the appendix of [1]. In all cases, the results we obtain confirm what was expected from a more general symmetry analysis.

Pentagon-prone orientation
Images of ψ 2 z , ψ 2 x , ψ 2 y and (ψ 2 x + ψ 2 y ) are shown for the pentagon-prone orientation in figure 2 (where, {x, y, z} are axes with z perpendicular to the surface, through the centre of a pentagon).
x + ψ 2 y orbitals of a pentagon-prone C − 60 ion. In this and subsequent figures, the y-axis is horizontal.
We will now use these images to describe the results obtained for JT-distorted C − 60 ions hopping between equivalent minimum-energy configurations.
We will first consider the parameters preferring a D 5d distortion. For 1 > 0, only well C gives a global minimum in energy. As a consequence, the appearance of an STM image does not depend on the strength of the surface interaction. Well C has an associated electronic state ψ z , so an STM image of this state will therefore be equivalent to that shown in figure 2(a). For the (doubly-degenerate) JT excited state, the result will be equivalent to that in figure 2(d). The result for the excited state is virtually identical to the image produced by a set of degenerate T 1u orbitals [5], because adding a contribution from ψ z as in figure 2(a) does not significantly alter the image. This means that this case cannot easily be distinguished from that for the neutral C 60 ion, even though there could be a strong JT effect present.
For 1 < 0, the five wells A, B, D, E and F are all equivalent lowest-energy minima, being distributed in an equivalent manner with respect to the z-axis. The coordinates of their electronic states track towards the x-y plane as the magnitude of 1 increases. Figure 3(a) shows the wells with no surface interaction as solid circles (red online), and plots the electronic coordinates as the magnitude of 1 increases as solid lines starting at the circles (blue online). The effect of the surface interaction on the coefficients a x , a y and a z is shown in figure 3(b). It can be seen that the limit of strong (negative) surface interaction is a doublet {ψ x , ψ y }. This limit is approached asymptotically, which means that different STM images are predicted for different magnitudes of the surface interaction. However, the images are dominated by ψ 2 x + ψ 2 y for most 1 < 0, so will resemble that in figure 2(d). The ψ 2 z contribution is only significant for small negative values of 1 .
For JT parameters preferring a D 3d distortion, results for 1 < 0 are somewhat similar to the case considered above. The minimum points are wells d, e, f, h and j, and their coordinates track asymptotically to the x-y plane, as for the D 5d case, as shown in figure 3(c) (cyan lines online). Again, it is not surprising that these wells all have the same energy from their distributions about the z-axis. As the coordinates are all very close to the x-y plane even with no surface interaction, there is very little change in a x , a y and a z as the magnitude of 1  as shown in figure 3(d). In all cases, any observed STM image will be the same as that in figure 2(d).
To show pictorially how the image equivalent to that of ψ 2 x + ψ 2 y arises in this case, we have presented images for each of the separate orientomers in figure 4. The five peripheral images correspond to the images expected if the system were to be permanently locked into the orientomers d, e, f, h and j. The central image corresponds to the expected STM image if the system is allowed to hop between these orientomers on a timescale faster than that of an STM experiment, and is a superposition of the five outer images.
The averaged image in the centre of figure 4 is much less distinctive than that of the individual orientomers. In fact, as it is an image of ψ 2 x + ψ 2 y , it is the same as the image expected from the doublet component of the LUMO of a neutral molecule split by a surface interaction, as shown in figure 7 of [5]. The increase in dynamic freedom afforded by allowing the system to hop between equivalent orientomers has obscured the ability to observe a unique signature of the JT effect using STM.  figure 2(a)) and will not depend upon the value of 1 as long as it is above the critical value of ≈0.921. However, below that value there is some significant x and y character in the predicted STM images.
We will now examine the case in which the C − 60 ion can pseudorotate between neighbouring wells (rather than hopping). For a D 5d distorted molecule, there is only one minimum-energy well when 1 > 0, so the molecule will not pseudorotate. It is possible to qualitatively determine the effect of the pseudorotation for the other cases by considering figures 1 and 3. Take as an example the case of a positive distortion favouring D 3d wells, where the system pseudorotates between wells {a, b, c, g, i}. The shortest path between nearest neighbours, such as wells a and g, passes through points with slightly larger z-component than the wells themselves. Furthermore, each point contributes a different amount of x and y. This means that when all points on the path are added up, the result is to even further enhance the contribution from z. The result is that the image taking pseudorotation into account has more character like that in figure 2(a) and less of that in figure 2(c). This results in a more rounded central pentagon with a less pronounced hole at its centre compared to the equivalent result that only includes hopping between wells. However, the effect of a surface interaction is also to favour z-character. Hence it is not possible to distinguish between the effect of the surface interaction and pseudorotational JT effects by examining experimentally-observed STM images. Consideration of a negative surface interaction for both D 5d and D 3d distortions leads to a similar conclusion. Quantitative calculations including contributions to the STM current from intermediate points between wells confirm our expectations from the qualitative discussion above.

Hexagon-prone orientation
We now turn our attention to the hexagon-prone orientation. As with the pentagon-prone case, we find that STM images arising from the LUMO when hopping between wells is taken into account can be composed from a linear combination of images of ψ 2 z and (ψ 2 x + ψ 2 y ), where z now refers to a C 3 -axis perpendicular to the surface and through the centre of a hexagon. The results for ψ 2 z , ψ 2 x , ψ 2 y and (ψ 2 x + ψ 2 y ) are all shown in figure 5. For a system with JT coupling parameters preferring D 5d distortions, wells A, B and D are favoured when 1 < 0. These wells are all close to the equator even in the absence of any surface interaction, so increasing the magnitude of 1 has only a very small effect on ψ tot . This is shown in figure 6(a) (cyan lines online) and 6(b).
For 1 > 0, the behaviour is rather different to that seen previously. The favoured wells are C, E and F. As 1 increases, these wells are pulled towards the poles, but at 1 ≈ 1.035 there is a sharp change in which the minimum-energy configuration suddenly jumps to the pole (i.e. the wavefunction becomes ψ z ). This is shown in figure 6(a) (blue lines online) and figure 6(b). This is because up to 1 ≈ 1.5456, there are two sets of local minima, one solution just involving ψ z and one set of three solutions involving ψ x , ψ y and ψ z . The Qs for each set are quite different. The energies of the two solutions are very similar, especially for 1 > 0.8. The solid lines on figure 6(b) are for the solution set that is lowest, with the dashed lines showing the result of the other solution.
For a system with JT coupling parameters preferring D 3d wells, well c, which is aligned along the z-axis, is the global minimum for all 1 > 0. This means that the results are the same whatever the magnitude of the surface interaction. For 1 < 0, the favoured wells are {a, b, e, f, h, j}. These are all pulled towards the equator as | 1 | increases, until a critical value of 1 is reached at which point the wells reach the equator, as shown by the curved paths (blue online) in figure 6(c). They then remain there as the strength of the surface interaction is increased further. This is illustrated in figure 6(d). The most noticeable difference between these results and the previous ones is that the paths are not around great circles. The shortest path round the sphere between two wells is also shown in figure 6(d) (dashed black line). The energy of this path is slightly higher than that of the curved path. Figure 6(c) also shows a path from two of the D 3d  ) shows the shortest path, which has slightly higher energy. The two dashed lines which do not track to a well (yellow online) show a path to a double bond, which also has slightly higher energy.
wells towards a single bond and on to a D 5d well (dashed lines, yellow online). Near to the D 3d wells, the energies of the two types of direct path are not much higher than the (curved) minimum-energy path. Competition between the two different shortest paths results in neither having a global minimum in energy.
As for the pentagon-prone case, the JT ground state when there is only one minimumenergy well (namely for well c for a D 3d distortion with 1 > 0) produces an image that is different to the image arising from degenerate T 1u orbitals. This is shown in figure 7(a). The image for the equivalent JT excited state is shown in figure 7(b). The images where there are degenerate minimum-energy wells are again very similar to those for degenerate T 1u orbitals, although two slightly different patterns emerge. The result for the JT excited state for hopping between wells C, E and F of a D 5d -distorted ion with 1 > 0, and for the JT ground state with

Double bond-prone orientation
Images of ψ 2 x , ψ 2 y , ψ 2 z and (ψ 2 x + ψ 2 y ) are shown for the double-bond prone orientation in figure 8. For this orientation, we have the additional complication that the surface interaction involves two parameters, 1 and 2 , and splits the LUMO into three singlets. However, as mentioned previously, we can expect | 2 | | 1 |. We will therefore start by considering the special case of 2 = 0. For 1 < 0, the favoured wells are E and F for parameters preferring a D 5d distortion, and e and f for D 3d parameters. All of these wells are already at the equator even with no surface interaction. Hence there is no dependence on the magnitude of 1 , as shown in figure 9. Predicted STM images show characteristic combinations of x and y character. For 1 > 0, the favoured wells are C and D for D 5d parameters, and a and b for D 3d parameters. In both cases, the wells track to the z-axis (in the x-z and y-z planes respectively). The most noticeable feature here is that the coefficients a x and a y will not be equal unless 1 is large and positive. This is because of the distribution of the wells with respect to the x and y axes.
The results for both D 5d and D 3d JT distortions and for both signs of 1 (with 2 = 0) are shown in figure 10. For 1 < 0 and with parameters favouring D 3d distortions, figure 9 shows that the contributions to the overall image, which do not depend on the magnitude of 1  dominated by ψ 2 y but with a small contribution from ψ 2 x . This can be seen in the resultant image ( figure 10(b)), where the lobes adjacent to the vertical centre line take on a 'horseshoe' shape, but don't quite join to give a pentagon as they would if there were more x-character.
When 2 is non-zero, some additional features become apparent in the results. The results for 1 < 0 and parameters favouring a D 5d distortion depend on the value of 2 , so are no longer unaffected by the strength of the surface interaction (which was the case when 2 = 0). The portion of the equatorial line on figure 11(a) near to the x-axis (cyan online) shows that for this case, the two minimum-energy wells (E and F) track towards the x-axis in the x-y plane as 2 becomes more negative, corresponding to a decrease in the contribution from ψ 2 y and a corresponding increase in ψ 2 x , as seen in figure 11(b). For small positive values of 2 , the same two wells are still minima, and their positions start to track towards the y-axis. This is It should be noted that we have specified the surface interaction in terms of 1 which alters the energy of ψ z relative to ψ x , and 2 which alters the energy of ψ y relative to ψ x . It is obviously possible to use alternative definitions, e.g. to specify energies relative to ψ z . The figures would then be different, but could still be interpreted in a similar manner.
For 1 > 0, images for this orientation allowing for free pseudorotation are more markedly different to equivalent images allowing for hopping than was the case for the pentagon-prone and hexagon-prone orientations. This is because the shortest path between wells C and D (for a D 5d distortion) or between wells a and b (for a D 3d distortion) both have the z-axis at their centre, as shown in figure 9, so intermediate points between the wells include significantly more z-character. Figure 12 shows images for the JT ground state involving wells C and D of a D 5ddistorted ion with a positive surface interaction ( 1 > 0) and with no splitting between x and y states ( 2 = 0). In figure 12, image (a) shows well C and (b) shows well D. (c) is the result taking into account hopping between the two wells; in other words, this is the image obtained by superimposing the separate images for wells C and D. Image (d) shows the result at an equivalent current allowing for pseudorotation between wells C and D, and, due to the intermediate points, includes contributions from ψ 2 z , as shown in figure 8(a). This is effectively a superposition of images obtained at all points on the shortest path between these two wells. While the images in (c) and (d) are visually different, it is still the case that a surface interaction could have the same effect, so the conclusion that pseudorotational effects cannot be distinguished from surface interaction effects still holds.
The origin of figure 12(d) is illustrated in the animated gif included as supplementary information to this paper (online supplementary data available from stacks.iop.org/NJP/14/ 083038/mmedia). This shows the pseudorotation from well C ( figure 12(a)) to well D ( figure 12(b)) via intermediate points on the pseudorotational path. Any observed STM image will be a superposition of all frames in this animation. The animation is what would be seen if it were possible to carry out an STM experiment on a faster timescale than that of the pseudorotation. Although it can't be conceived that such an experiment would be possible in practice, it does illustrate how the contribution from ψ 2 z becomes apparent in the image with pseudorotation but not in that with hopping.
For 1 < 0, the shortest path between wells E and F, relevant to a D 5d -distorted ion, passes though the x-axis. Thus equivalent arguments to above lead to the conclusion that an STM image taking into account pseudorotation will have a stronger x-like character than the equivalent hopping case. Similarly, images taking into account pseudorotation between wells e and f of a D 3d -distorted ion will have a stronger y-like character. Again, pseudorotational effects cannot be distinguished from surface interaction effects.

Approximate method
The method we have used in the previous sections has been to calculate values for the normal mode coordinates Q γ that result in minimum-energy configurations for any given set of JT and surface interaction parameters, which requires minimizing the energy of the lowest eigenvalue of a 3 × 3 matrix (involving H JT + H s ) with respect to five different variables. A much simpler but approximate method is to calculate the Q γ in the absence of any surface interactions, and then calculate new electronic states by finding the lowest eigenvalue of the 3 × 3 matrix but without recalculating the Q γ .
The positions of the wells obtained using this approximate method follow the same paths on a sphere as those in figures 3, 6, 9 and 11, with the exception of figure 6(c), for D 3d wells in the hexagon-prone orientation, where there are two sets of competing minima with similar energies but very different Q γ values. The exact calculation resulted in curved paths tracking to the equator, whereas the approximate method results in lines heading directly to the equator. The plots of the coefficients a x , a y and a z are also very similar in all cases, with both the same initial behaviour and same strong surface-interaction limits. The main difference is in the rate at which the strong surface-interaction limits are attained. The exact calculation tends to give coefficients that vary almost linearly with the strength of the surface interaction before attaining the strong surface interaction values, whereas the approximate method tends to result in the coefficients attaining their strong surface-interaction limits asymptotically. Nevertheless, the coefficients cover the same ranges of values, so if it is determined that a certain set of coefficients is needed to match some experimental data, then matches will be found using both methods. The only difference will be that the two methods predict the match to occur at different values of the surface interaction constants.
While the approximate method will not be explored any further in this paper, it will be useful in future studies of JT effects in higher charge states n of the fullerene ion C n− 60 , where the matrix to be diagonalized is of a higher dimension, namely 6 for the low-spin states with n = 2 and 4, and 8 for n = 3 [15]. The approximate method could be used to quickly explore the whole parameter space in order to identify potential matches to experimental data, and then the more exact method could be used to confirm those results.

Comparison with experimental results
We will now show how our theoretical results can be used to help explain experimentallyobserved images. We can determine a form for the current I in terms of ψ x , ψ y and ψ z needed to obtain the best match to any experimental image. We can then relate these relative contributions to our figures in the previous section in order to deduce information on possible orientations and ranges for the JT coupling constants and surface interaction parameters. The majority of STM images of C 60 involve molecules that are intrinsically neutral. However, the states that are involved in the transmission of current, and hence the states that are imaged, are not necessarily those of the neutral molecule. Transfer of charge from the STM tip to the molecule could result in states of the negatively-charged ion being imaged. Alternatively, a metallic surface could provide some charge [19], although it seems unlikely that charge transfer from the surface would result in each molecule being singly charged. As the charge state may not be obvious, it is worthwhile examining images in papers that don't explicitly purport to arise from C − 60 ions, especially those of C 60 molecules adsorbed on a buffer layer so that surface interactions are less likely to dominate over JT effects. [1] attempted to find matches to configurations arising from a static JT effect in C 60 ions adsorbed on a buffer layer, without considering the dynamics that result in equivalent configurations being imaged together. The results could also be explained if the JT effect is intrinsically dynamic but there are other perturbations that result in one well being favoured over others. It suggested that the most likely candidates for observing JT effects are the images of C 60 molecules on an alkylthiol self-assembled monolayer (SAM) on a Au(111) substrate in [20]. The experimental images from [20] are reproduced in figure 13(a). We will now look whether a good match to the experimental images can be obtained using our dynamic JT model. This will provide an alternative to the model of neutral C 60 molecules proposed in [20].
The left-hand image in 13(a) is rather distinctive so we will concentrate on this first. Examining the images for ψ 2 x , ψ 2 y and ψ 2 z for the pentagon, hexagon and double-bond prone orientations suggests that the image arises from predominantly the ψ 2 z component of the double-bond prone orientation, with a small contribution from ψ 2 y . However, it is clear that the experimental images involve molecules with a small tilt of their C 2 z-axis away from the axis perpendicular to the surface. In [20], the authors estimated the tilt to be around 1.9 • from the highest carbon atom, which is at the upper corner of a pentagon. In terms of our angles, this corresponds to a tilt of around 9.3 • from the C 2 z-axis in the x-z plane, towards a fivefold axis. We will investigate the results that would be obtained with small angles of tilt in either the x-z plane or the y-z plane.
The tilt has implications for the modelling process. Due to the tilt, the only remaining molecular symmetry operation is a reflection in the plane of the tilt, indicating a reduction in symmetry to the group C s . As a result, cross-terms in the expression for the STM current will not necessarily cancel. For example, with a tilt in the x-z plane, the application of the σ h reflection operation of the C s group (which reflects in the x-z plane) to a general electronic state c x ψ x + c y ψ y + c z ψ z produces an equivalent-energy state c x ψ x − c y ψ y + c z ψ z , such that the sum of squares of these states will contain an additional term in ψ x ψ z when c x and c z are nonzero. This term wasn't present in the higher-symmetry C 2v situation of the double-bond prone orientation with no tilt because the additional rotation and reflection produced two more wells (when the coefficients are all non-zero), such that when summing over all four wells the crossterms all cancel. Similarly, a tilt in the y-z plane could result in cross-terms in ψ y ψ z .
Before starting our calculations, we note that interactions between neighbouring ions are likely to be important in the images of [20]. In fact, the authors suggest that these interactions dominate and surface interactions are negligible. Although we haven't included neighbour interactions in our model, our results can be extended to cover this situation also for the doublebond prone and tilted orientations. In the left-hand image in 13(a), the C 60 molecules sit in a hexagonal lattice with either our x-z plane or our y-z plane aligned with the diagonal of a hexagon. We will first consider the case where the molecules are adsorbed in pure doublebond prone orientations (i.e. with no tilt). Although the actual symmetry should be described by a space group, the point group symmetry operations that can be applied with neighbour interactions but no surface interactions are C 2 rotations about the x, y and z axes and reflections in the x y, yz and zx planes (plus inversion), resulting in a symmetry of D 2h . In this symmetry, {ψ x , ψ y , ψ z } transform as {B 3u , B 2u , B 1u } for a tilt in the x-z plane or {B 2u , B 3u , B 1u } for a tilt in the y-z plane. Hence, although the point group symmetry is D 2h rather than C 2v , the effect of both the surface interaction and the interactions with neighbours is to separate the T 1u orbital into the same three triplets. Hence the same Hamiltonian can be used to describe the neighbour interactions as we have used to describe surface interactions. When the molecules are adsorbed with a tilt of the C 2 axis, the symmetry is again C s , the same as with a surface interaction, so again both effects can be described by the same Hamiltonian.
It is a simple matter to perform our calculations with different angles of tilt (in either the x-z plane or the y-z plane) and look for a match with a current containing contributions from ψ 2 x , ψ 2 y , ψ 2 z and either ψ x ψ z or ψ y ψ z as appropriate. After a detailed investigation, we find that the best match to the experimental image in the left-hand side of figure 13(a) is obtained when the C 60 molecule is tilted in the y-z plane so that its molecular C 2 Z -axis is ≈0.1 rad (≈5.7 • ) from the axis perpendicular to the surface. The corresponding current is I ≈ 0.13ψ 2 y + 0.87ψ 2 z with respect to the molecular C 2 Z -axis. The result of our simulation with this current is shown in figure 13(c). This is a closer match to the experimental image than the authors own simulation, shown as the inset to the left-hand image in figure 13(a), where the central feature is much more rounded, as it was taken to come from a pentagonal face of the C 60 molecule rather than a double bond. We note that our best-fit angle of tilt both has a different value and is in a different plane to that in [20], who (using our terminology) assumed a tilt in the x-z plane of β ≈ 9.3 • . However, it is not surprising that we need a different angle of tilt because the central feature has a different origin.
Obviously, as we are only looking for a qualitative match to the experimental image, the numerical values in I are not fixed precisely. However, as the image is so distinctive, the parameters can't be varied significantly from these values. With respect to a tilted C 2 z-axis, the lowest limit on the ψ 2 y contribution comes from I ≈ 0.075ψ 2 y + 0.925ψ 2 z ; for smaller contributions, the two leftmost lobes have become too small. Conversely, the upper limit on the ψ 2 y contribution comes from I ≈ 0.2ψ 2 y + 0.8ψ 2 z ; for larger contributions, the two right most lobes round the central feature have become too prominent. It is also relevant to note that for a current corresponding to ψ 2 z alone, the uppermost and lowermost features are crescent-shaped rather than each being two distinct lobes, as can be seen in figure 8(a) for the non-tilted case. Finally, we note that we do not want any cross-terms in the expression for I to get a match to the experimental image. This means that it is not possible to get our match from a model of neutral C 60 , as the current from any state that is any linear combination of ψ y and ψ z will necessarily contain unwanted cross-terms (in ψ y ψ z ). It is only when two equivalent configurations are imaged together that the cross-terms cancel. No combination of degenerate states could give the required form for I either.
The next step in the identification of the origins of the experimental image is to determine what values of the JT and surface interaction parameters could produce the required expression for I . This can be done using our graphs of a x , a y and a z in section 3 for the pure doublebond prone orientation. Although those graphs are without a tilt, they indicate which ranges of parameters are likely. These situations can then be explored further with explicit calculations including the tilt. The graphs confirm that the only situation that can result in differing amounts of ψ 2 x and ψ 2 y are from the double-bond prone orientation, which is consistent with our conclusion that the image arises from a small tilt from this orientation.
Before we can proceed further however, we need to consider a further implication of the reduction in symmetry to C s due to a tilt in the y-z plane. In this case, ψ x transforms as A , but ψ y and ψ z both transform as A . As ψ y and ψ z can't be distinguished on group theory grounds, it is not possible to write down a unique surface interaction Hamiltonian from group theory alone. One possible form of the surface interaction Hamiltonian is simply that in equation (2.3), in which case the expression for the current would be exactly the same as we obtained above for the double-bond prone orientation, but the image would appear different because it would be viewed from a different angle. However, it is not clear why interactions with the surface and neighbours would produce a Hamiltonian of this form. Another form would be diagonal in the basis of a z-axis perpendicular to the surface and x and y axes in the plane of the surface. This would then be rotated by 0.1 rad to apply to a molecular C 2 Z -axis. A more general form of Hamiltonian would have two mutually orthogonal but otherwise arbitrary linear combinations of ψ y and ψ z taken to be at different energies, although this introduces an additional parameter in the model.
To get an image that predominantly involves ψ y and ψ z , the figures show that it must originate from D 3d wells a and b. We have already noted that, while it is not known for certain what distortion symmetry a C 60 molecule prefers, D 3d symmetry appears to be the most likely [7][8][9]. From the positions of the wells as shown in figure 1, we can see that with a tilt in the y-z plane away from the C 2 z-axis, wells a and b will have slightly different energies. In order to determine whether the difference in energy is significant, we will need to consider the JT model with specific parameters. As one state is of the form c y ψ y + c z ψ z and the other is of the form −c y ψ y + ca z ψ z , there will be no cross-terms to consider if the two states are imaged together.
If the current only involves ψ y and ψ z , we can reduce the 3 × 3 matrix representation of the total Hamiltonian to the 2 × 2 block involving ψ y and ψ z only. It is easy to see that this block can be extracted from the full matrix if Q 5 = Q 6 = 0 (which is what we expect for a solution involving D 3d wells a and b [13]) and the surface/neighbour interaction puts ψ x at a higher energy than ψ y and ψ z . Also, this means that it is only necessary to consider the separation = 1 − 2 between ψ y and ψ z , rather than the values of the two parameters themselves. Investigating the dependence of the simulated image on V 2 , V 3 and shows that, as long as wells a and b can both be imaged together, good matches can be obtained with all values of V 2 and V 3 that favour D 3d wells and with a very small value of . For typical parameters of V 2 = 0, V 3 = 0.5 and = 0.01 and assuming a surface/neighbour interaction Hamiltonian that is diagonal in a basis with z perpendicular to the surface, our calculations show that the difference in energy between the two wells is 0.0013V 2 1 /(µω 2 ). In [9], the JT energy was estimated to be 57.94 meV. As the JT energy is V 2 1 /(5µω 2 ), this suggests that the difference in energy between wells a and b is around 0.38 meV. [20] doesn't give any details of the type of tip they used, so it is difficult to determine the energy resolution in their experiments. However, the energy resolution of an STM experiment at temperature T with an ideal sharp metal tip is ≈5.4k B T [21]. At the 5 K temperature used to obtain the images in [20], this would suggest that the energy resolution would be around 2.3 meV, which is an order of magnitude larger than the difference in energy between the two wells. As non-sharp tips are likely to have a lower energy resolution, we can probably safely conclude that if the images in [20] do involve the JT system, they will be of wells a and b imaged together.
Further exploring the dependence of our results on indicates that can be positive or negative, but must have a sufficient magnitude for wells a and b to be sufficiently separated in energy from the other wells that only these two wells will be imaged. As the expression for I is not very sensitive to the value of , we find that any value of in the range −0.03 < < 0.16 is acceptable. Note that the optimal value of ≈0 is not in the middle of this range, and only has a very weak dependence on V 2 and V 3 .
In the absence of surface interactions and nearest-neighbour interactions, D 3d wells are obtained if V 3 > 3V 2 / √ 5, as shown in the appendix of [1]. We must have V 2 < −5/4 √ 2 and V 3 < 3 √ 5/4 √ 2 for the system to be stable. As the optimal value of is very small, the acceptable ranges of V 2 and V 3 are still very close to this in the current problem.
Obviously, as we don't have a unique form for the surface interaction in this case, there could be other matches to different JT parameters and well combinations, but the key result is that we have shown that it is possible to get matches with a JT model.
We will now turn our attention to the image from [20] reproduced in the right-hand half of figure 13(a). The authors propose that their image is from molecules adsorbed with a C 2 axis perpendicular to the viewing plane, and our calculations confirm this. Hence we do not need to consider the additional complications due to tilt that were necessary to explain their left-hand image. However, the molecules are rotated in the x-y plane on their hexagonal lattice positions, by slightly less than 90 • with respect to the orientation of the molecules in the left-hand image. Due to this lowering of symmetry, the only point group operations that survive are a C 2 rotation about the z-axis and a reflection in the x-y plane (plus inversion), resulting in a point group symmetry of C 2h to describe neighbour-interactions. In this symmetry, ψ z transforms as A u , and ψ x and ψ y both transform as B u . Hence this is another case where we can't use group theory alone to write down a unique Hamiltonian.
As the molecule is adsorbed in the double-bond prone orientation, the current must be a linear combination of ψ 2 x , ψ 2 y and ψ 2 z with no cross-terms. However, this image is rather generic and simulated images obtained for a wide range of different numerical values for the coefficients cannot be distinguished. It is probably an image of predominantly ψ 2 x , and indeed an image of ψ 2 x alone does closely resemble the experimental image, as shown in 13(d). This is a lower current (and rotated) version of the simulation in figure 8(b). However, adding in even a fairly large contribution from ψ 2 z may not significantly alter the image (although this does depend on the resolution, i.e. the magnitude of the current), to the extent that the possibility of the image being from 0.5(ψ 2 x + ψ 2 z ) cannot be ruled out. However, we can conclude that the contribution from ψ 2 y must be smallest. An upper limit on the y contribution is ≈0.2ψ 2 y ; for larger contributions, the two lobes start to become too crescent-shaped, and new lobes at the top and bottom of the image start to appear. The fact that the ψ 2 x and ψ 2 y contributions are different confirms that the image must arise from the double-bond prone orientation.
Images obtained from C 60 ions in a K 5 C 60 monolayer [22] or multilayer [23] should involve JT distortions of C 5− 50 ions to D 5d or D 3d symmetry that are the same as those in C − 60 as the two JT problems are isoelectronic. We could therefore expect them to be described using the results in this paper. However, current images in the literature are not sufficiently distinctive for an unequivocal identification to be made, and also the results would need to be described in terms of different (i.e. multi-electron) states. Hence no match is attempted here.

Summary and discussion
In this paper, we have used HMO theory and a symmetry analysis to investigate the influence of a combination of surface interactions and the dynamic JT effect in a C − 60 ion on the appearance of the ion in STM images. We have obtained results for when the C − 60 ion is adsorbed onto a surface substrate with a pentagon, hexagon or double-bond facing the surface. We found that any STM image will be a superposition of images of just ψ 2 z and (ψ 2 x + ψ 2 y ) for pentagon-prone and hexagon-prone molecules, and images of two of ψ 2 z , ψ 2 x and ψ 2 y for double-bond prone molecules, with no cross-terms possible. This gives a simple way of interpreting STM images observed experimentally, without the need for more complex DFT calculations.
We have also included interactions with neighbours for the double-bond prone orientation. This allowed us to provide an alternative explanation of the images in [20], which involves C − 60 anions subject to dynamic JT effects, adsorbed on the surface substrate in different orientations, interacting with their neighbours and also possibly subject to interactions with the surface. As the molecules have different orientations in the two cases, the neighbourinteraction Hamiltonian will be different, which is why more than one type of image is possible. It is unlikely that this situation could be described by a surface interaction but no neighbour interactions though, because both sets of images involve a double bond or nearly double-bond prone orientation, so both are described by very similar surface interaction Hamiltonians. Also, we can't get our match to the image on the left-hand side of figure 13(a) without a JT effect (e.g. for neutral C 60 ), as it is a particular combination of ψ 2 y and ψ 2 z with no cross-terms. Surface/neighbour interactions alone could result in states that are a combination of ψ y and ψ z , but there would be unwanted cross-terms, as it is only when we superimpose images of equivalent wells that such cross-terms cancel. Similarly, we can't get a match assuming a static JT effect in C − 60 ions as any such image would again contain unwanted cross-terms. The only way we could get a match to the right-hand image in figure 13(a) would be if the lowest-energy state is just ψ x , or the ψ x and ψ z states are accidentally degenerate (where the current would contain equal combinations from ψ 2 x and ψ 2 z ). Our alternative explanation of the images in [20] is not decisive proof that JT effects are present and the model in [20] is not valid. However, it is important to note that the two possible explanations have very different physical implications and different potential matches to observed images. Our model is not simply an alternative explanation of the same match. Therefore care must be taken in analysing such images in the knowledge that there may be more than one potential explanation for any given result. It should also be noted that it could still be possible that STM images of fullerene ions in different experimental setups to that in [20] should be explained using a static JT model, as was done in [1]. Even if the intrinsic JT coupling in C 60 ions would lead to dynamic effects, other factors such as interactions between C 60 ions could remove the equivalence between JT distortions and hence suppress the tunnelling between wells.
A consideration of a combination of surface and neighbour interactions for molecules adsorbed in the pentagon or hexagon prone orientations is beyond the scope of this paper. However, a method of treating such interactions by writing down a Hamiltonian to reflect the symmetry reduction due to the interactions with neighbours could be followed in the manner we have presented for the double-bond prone orientation. This could potentially explain the images of C 60 molecules on an alkylthiol SAM shown in figure 2 of [24], where the results cannot be explained by a combination of ψ 2 z , ψ 2 x and ψ 2 y , although an alternative explanation is that the images are a signature of a static JT effect [1]. If the surface substrate is of a lower symmetry than we have assumed, such as if there are surface corrugations, then the equivalence between JT distortions could also be removed.
We have considered the case where the system is assumed to 'hop' between equivalent minimum-energy distortions without spending any appreciable time in intermediate configurations, and also examined the effect of pseudorotation, where the paths connecting the distortions is also taken into account. We found that while there are some subtle differences in the results for the two cases, the major features of any STM image are unaltered, and that the pseudorotational effects cannot be separated from those of the surface interaction.
As well as splitting MOs, the JT effect displaces the carbon nuclei. Our treatment has ignored the effect of these displacements on STM images, as these are expected to be negligibly small. Significant displacement of the nuclei during pseudorotation would 'blur' the images presented here, as in the case of an E ⊗ e system [2].
We have seen that even a limited amount of freedom to tunnel between equivalent distortions is sufficient to prevent identification of features that might be ascribed to the JT effect. It is possible that other perturbations could act on such a system in order to lock it into one particular distortion. This means that it could be possible for STM images in the literature to match theoretical images of individual distortions. Also, conversely it is possible that images previously identified as being due to neutral C 60 could actually arise from JT-distorted ions, where dynamic effects have removed the obvious characteristics of the JT effect.