text:
# May 2003 Preprint MPI-PTh/2003-23 Inversion of the Photon Number Integral

## I Introduction

The “photon number integral” , called $`n`$, is a construction whereby one can define the number of photons radiated by a charged particle following a prescribed trajectory in space-time (minkowski space). For $`n`$ to be finite the trajectory must obey two conditions: smoothness and equality of initial and final trajectories. The integral then returns a real number for any such curve. The idea can be extended from minkowksi to euclidean space where it also becomes possible to consider closed curves. The identification of initial and final velocities (tangents) is then automatic.

An important question left open in  concerns the symmetries of $`n`$. Given $`n(C)`$ for a curve $`C`$ for what other curves $`C^{}`$ do we have $`n(C)=n(C^{})`$? There are the evident invariances under rotations (lorentz transformations), translations, and reflections. Furthermore, since there are no dimensional quantities except the path itself involved in the problem, there is also an evident invariance under scaling, that is rescaling of all coordinates simultaneously, $`x_\mu \lambda x_\mu `$, where $`\lambda `$ is a constant.

Nevertheless the question remains as to further invariances, in particular with respect to conformal transformations. These are closely related to scale transformations, and indeed conformal invariance originally entered physics as a property of electrodynamics, and ought to be expected to apply to photons. Also, in somewhat different contexts the properties of integral expressions like ours have been considered in connection with inversion , or “Mobius” invariance , which is tantamount to conformal invariance. Here we would like present a discussion of inversions for $`n`$.

The expression in question, in minkowski space, is

$$n=𝑑x_\mu \frac{1}{S_{iϵ}^2}𝑑x_\mu ^{}$$
(1)
Since our question here is essentially a mathematical one, we have dropped the (dimensionless) electromagnetic coupling constant appearing in the original expression for photons.

$`S_{iϵ}`$ is the four-distance between the points $`x,x^{}`$ in the following way

$$S_{iϵ}^2=(tt^{}+iϵ)^2(𝐱𝐱^{})^2$$
(2)
The charm of the expression Eq  is that despite the possible singularities it is actually finite. The possible singularity at $`S^2=0`$ is handled by the $`iϵ`$ and the possible divergence at infinity by the “straight line condition” for the equality of initial and final paths, as may be verified from the fact that $`n`$ is zero for the simple straight line . This means it is well defined as it stands and needs no “regularizations” or “subtractions”. In Ref  a regularization is introduced by supersymmetry where a scalar particle cancels the singularity of a vector particle (“gluon”) propagator. In Ref  the singularity is cancelled by subtracting a second term where the straight line distance $`S`$ between $`(x,x^{})`$ is replaced by the arc length along the curve. However as we anticipate from its physical origin, Eq is finite as it stands. Indeed using the identity

$$_\tau _\tau ^{}ln(S^2)=+\frac{4}{S^4}(\mathrm{\Delta }_\mu \frac{dx_\mu }{d\tau })(\mathrm{\Delta }_\nu \frac{dx_\nu }{d\tau ^{}})\frac{2}{S^2}\frac{dx_\mu }{d\tau }\frac{dx_\mu }{d\tau ^{}}.$$
(3)
we were able to rewrite the integral as

$$𝑑x_\mu \frac{1}{S_{iϵ}^2}𝑑x_\mu ^{}=2𝑑x_\mu \frac{\delta _{\mu \nu }\frac{\mathrm{\Delta }_\mu \mathrm{\Delta }_\nu }{S^2}}{S^2}𝑑x_\nu ^{}.$$
(4)
Delta is the vectorial distance between the two points, $`\mathrm{\Delta }_\mu =x_\mu x_\mu ^{}`$. We can rewrite this in a suggestive form if we introduce the “transverse vector”

$$dx_\mu ^T=dx_\mu \frac{\mathrm{\Delta }_\mu (\mathrm{\Delta }dx)}{S^2},$$
(5)
which is $`dx`$ with the “longitudinal part” removed, i.e. has the property $`dx^T\mathrm{\Delta }=0`$. This is suggestive of the transversality property of physical photons and it might be said that our expression is finite because it only contains the radiated and not the “coulomb” photons. In any event, the expression is now manifestly non-singular, as may be seen by expanding the numerator, see Ref .

Observe that this procedure does not introduce any dimensional quantities, so the expression is still scale invariant. Apparently the $`iϵ`$ is harmless in this respect. This manifestly non-singular expression, without the $`iϵ`$, may be used to define $`n`$ in case of doubt and we will use it in the following (mostly in its euclidean version) taking $`n`$ as:

$$n=2𝑑x_\mu ^T\frac{1}{S^2}𝑑x_\mu ^T.$$
(6)
## II Inversion

We shall concentrate on the inversion operation

$$x_i\frac{a^2}{x^2}x_i.$$
(7)
The full conformal group is generated by adding these inversions or “Mobius transformations” to the usual translations and rotations. To keep the physical dimensions in order we have introduced a length constant $`a`$ parameterizing the operation. The inversion is around a point $`𝒪`$ , which we take as the origin. As we shall see below, the cases where the center of inversion is on the curve $`C`$ itself is of special interest. The index $`i`$ means any of the coordinates, and in minkowksi space includes an inversion of the time coordinate $`x_0`$. We shall mainly focus however on the simpler case of euclidean space. There we use boldface notation for vectors, so inversion is $`𝐱\frac{a^2}{x^2}𝐱`$.

We thus begin by considering the formal properties of the euclidean expression

$$n=\frac{\mathrm{𝐝𝐱𝐝𝐱}^{}}{S^2}=\frac{\mathrm{𝐝𝐱𝐝𝐱}^{}}{(𝐱𝐱^{})^2}$$
(8)
under inversion. The (-) sign is the natural choice that makes n positive in euclidean space. The curve over which the integrations are performed may either be a smooth closed curve, or a smooth infinite curve which becomes the same straight line at $`\pm \mathrm{}`$.

## III formal inversion

It is illuminating to start with the relation Eq . We use a differential notation, where $`d_xf(𝐱,𝐱^{})=_xf\mathrm{𝐝𝐱}`$, so that with $`S^2=(𝐱𝐱^{})^2`$

$$d_xlnS^2=\frac{2(𝐱𝐱^{})\mathrm{𝐝𝐱}}{S^2}d_x^{}lnS^2=\frac{2(𝐱𝐱^{})\mathrm{𝐝𝐱}^{}}{S^2}.$$
(9)
Eq  then takes the form of the curious and interesting identity

$$\frac{1}{2}\left[d_xd_x^{}lnS^2+d_xlnS^2d_x^{}lnS^2\right]=\frac{\mathrm{𝐝𝐱}\mathrm{𝐝𝐱}^{}}{S^2}.$$
(10)
In this way we express the integrand of Eq  in terms of certain differentials with simple transformations under inversion. Applying the inversion to $`S^2`$

$$S^2=(𝐱𝐱^{})^2\frac{a^2}{x^2}\frac{a^2}{x^2}S^2$$
(11)
and so

$$lnS^2ln\frac{a^2}{x^2}+ln\frac{a^2}{x^2}+lnS^2$$
(12)
Hence if we insert the substitution Eq  in the lhs of Eq  the first term is unchanged while for the second

$$\frac{1}{2}\left[d_xlnS^2d_x^{}lnS^2\right]\frac{1}{2}\left[d_xlnS^2d_x^{}lnS^2\right]+\frac{1}{2}\left[d_xlnx^2d_x^{}lnx^2+d_xlnS^2d_x^{}lnx^2+d_x^{}lnS^2d_xlnx^2\right]$$
(13)
In other words our fundamental form transforms additively under inversion

$$\frac{\mathrm{𝐝𝐱}\mathrm{𝐝𝐱}^{}}{S^2}\frac{\mathrm{𝐝𝐱}\mathrm{𝐝𝐱}^{}}{S^2}+I$$
(14)
where we call the additional quantity $`I`$. In this way we recover the results of Ref for the gluon propagator.

## IV inversion of the explicitly finite integrand

So far we have ignored the singularity in Eq . We return to the explicitly finite form Eq  as the definition of $`n`$. Since the explicitly non-singular forms were found (see Ref ) by subtracting 1/2 of Eq  from Eq , we do the same here and Eq  becomes

$$d_xd_x^{}lnS^2\frac{1}{2}d_xlnS^2d_x^{}lnS^2=\frac{2\mathrm{𝐝𝐱}\mathrm{𝐝𝐱}^{}}{S^2}\frac{2[\mathrm{𝐝𝐱}(𝐱𝐱^{})][(𝐱𝐱^{})\mathrm{𝐝𝐱}^{}]}{S^4}=\frac{2\mathrm{𝐝𝐱}^T\mathrm{𝐝𝐱}_{}^{}{}_{}{}^{T}}{S^2}.$$
(15)
where the rhs is now the explicitly non-singular integrand of Eq . The difference between the singular and non- singular forms is simply the coefficient of the first term on the left. However, since this term is in any event identically invariant under inversion, Eq  still holds, that is

$$\frac{2\mathrm{𝐝𝐱}^T\mathrm{𝐝𝐱}_{}^{}{}_{}{}^{T}}{S^2}\frac{2\mathrm{𝐝𝐱}^T\mathrm{𝐝𝐱}_{}^{}{}_{}{}^{T}}{S^2}+I,$$
(16)
with the same $`I`$ as in Eq .

## V properties of $`I`$

$`I`$ is defined as

$$I=\frac{1}{2}\left[d_xlnx^2d_x^{}lnx^2+d_xlnS^2d_x^{}lnx^2+d_x^{}lnS^2d_xlnx^2\right],$$
(17)
where the origin is at the center of inversion $`𝒪`$ . If the origin is placed elsewhere, $`x`$ is the distance to $`𝒪`$ . Evidently $`I`$ is symmetric

$$I(x,x^{})=I(x^{},x).$$
(18)
Furthermore we note an important property arising from the fact that two successive applications of the inversion is the identity operation. Applying an inversion to the rhs of Eq  or Eq  again, we should get the original expression. We thus conclude that under inversion

$$II,$$
(19)
which one can also check by directly inserting Eq  into Eq .

$`I`$ is a scalar under rotations, but not under translations (holding $`𝒪`$ fixed) because of the presence of $`x^2`$, the distance from $`𝒪`$ . It is invariant under rescalings $`𝐱\lambda 𝐱`$ because of the presence of the differentials. For the same reason it is independent of the parameter $`a`$ giving the radius of the sphere of inversion. That $`I`$ does not contain the parameter of the inversion suggests that it can only depend on some global properties of the operation.

## VI  integration

To study the behavior of $`n(C)`$ under inversion we proceed according to the following steps. We carry out a change of variables according to Eq  in the integral for $`n(C)`$. This results simply in an identity with a new curve $`C_{inv}`$ in the new variables. The integrand receives an additional term, $`I`$. which is also to be integrated over the new curve $`C_{inv}`$. That is, we have one term with an integral of the desired expression and one with $`I`$, so that

$$n(C)=n(C_{inv})+\underset{C_{inv}}{}I$$
(20)
If we can show that $`\underset{C_{inv}}{}I`$ is zero then $`n(C)=n(C_{inv})`$.

As another expression of the fact that two successive operations with Eq  are the identity we can invert $`n(C_{inv})`$ once more to return to the original curve, giving the relation

$$\underset{C}{}I+\underset{C_{inv}}{}I=0.$$
(21)
This of course just amounts to Eq  if we change variables again to make $`C_{inv}C`$. An evident property following from Eq  is that if $`C_{inv}C`$, as for a circle with the center of inversion in the center, or a closed curve in 3-space on the surface of the sphere of inversion, then

$$\underset{C}{}I=0CC_{inv}.$$
(22)
As opposed to Eq  this is not merely an algebraic identity but involves of course the nature of the curves. It is evidently true by Eq  if $`C`$ and $`C_{inv}`$ are identically the same curve. But it is clearly also true in some more general sense, say if one curve is simply a rotation, translation or rescaling of the other. The wider meaning of “$``$” in $`CC_{inv}`$ is an interesting question and will be further discussed below (section XI E).

## VII  two curves

An amusing generalization of these properties of the integration suggests itself. We mention it although it lies somewhat outside our main topic. The quantity $`I`$ is a functional of one curve $`C`$. But actually it could be regarded as a functional of two curves $`C`$ and $`C^{}`$. We might have different curves in the $`x`$ and $`x^{}`$ spaces. That is, we consider the integral $`_C_C^{}I`$.

Because of Eq  it doesn’t matter to which variable the curve is assigned:$`_C_C^{}I=_C^{}_CI`$. Now let $`C^{}`$ be the inversion of $`C`$ so $`C^{}=C_{inv}`$. Then under inversion according to Eq  the integral goes to minus itself. Therefore

$$_C_C^{}I=0C^{}=C_{inv}.$$
(23)
Eq  may then be read as the special case where $`C=C_{inv}`$.

A place where Eq  might be useful concerns the question of the additivity of the parts of a curve. In general of course $`I`$ for a single curve cannot broken into two parts such that the total integral is the sum of the integrals for the two parts. $`I`$ is a bilocal object, and there will usually be cross terms between the two parts. However, if the two parts in question are the inversion of each other, then according to Eq , the cross terms are zero and the integrals for the two parts may be simply added.

To develop this idea properly a study of the possible singularities, which may be different than in the single-curve case, would be necessary. We stress that our integrals are always meant over just one curve in the original sense unless explicitly indicated as in Eq .

## VIII  Possible Singularities of $`I`$

There are singularities or potential singularities of $`I`$ which need to be understood. We will come to the conclusion that, despite appearances, $`I`$ is a non-singular object— for a given curve.

We would like to argue that for the investigation of possible singularities at a point it suffices to consider the behavior of $`I`$ for a straight line through that point. Consider the smooth curve $`𝐱(\tau )`$ in the neighborhood of the point $`𝐱(0)`$. It is convenient to introduce the “natural” parameterization where $`\tau `$ is the length along the curve (proper time, in minkowski space). Then $`𝐱𝐱(0)+\dot{𝐱}\tau +\ddot{𝐱}(1/2)\tau ^2+..`$, where $`\dot{𝐱},\ddot{𝐱},\mathrm{}`$ are the first, second,… derivatives with respect to $`\tau `$. Due to the choice of $`\tau `$ as the path length $`(\dot{𝐱})^2=1`$, $`\dot{𝐱}\ddot{𝐱}=0,`$ and so forth. Now $`𝐱(\tau )𝐱(\tau ^{})`$ is odd with respect to interchange of $`\tau ,\tau ^{}`$. This leads to $`(𝐱(\tau )𝐱(\tau ^{}))^2(\tau \tau ^{})^2\left[(\dot{𝐱})^2+(\ddot{𝐱})^2b(\tau ,\tau ^{})+\mathrm{}\right]`$, where $`b`$ is a bilinear expression in $`\tau ,\tau ^{}`$.

Since, with $`(\dot{𝐱})^2=1`$

$$ln\left[(\tau \tau ^{})^2(\dot{𝐱}^2+(\ddot{𝐱})^2b(\tau ,\tau ^{})+\mathrm{})\right]=ln(\tau \tau ^{})^2+ln\left[1+(\ddot{𝐱})^2b(\tau ,\tau ^{})+\mathrm{}\right]$$
(24)
we see that the possibly singular behavior of the logarithm results from the “velocity” or tangent term $`\dot{𝐱}`$, while the dependence on the “acceleration” or curvature is non-singular for $`\tau ,\tau ^{}0`$. Similarly for $`lnx^2`$ near $`x^2=0`$ we have $`lnx^2ln\tau ^2+ln(1+(\ddot{𝐱})^2b(\tau )+\mathrm{})`$. We may thus investigate the possible singularities by looking at the behavior of $`I`$ for a straight line.

We first consider the possible singularity for $`S^2=0`$, while assuming $`x^20`$, that the center of inversion $`𝒪`$ is not on the curve.

Choose the “1” axis such that it is parallel to the tangent $`\dot{𝐱}`$ at the point. Let the distance of the projected tangent from $`𝒪`$ i.e. its “impact parameter”, be $`\delta `$ (Fig 1). Now, for the investigation of the possible singularity we treat the curve as this straight line, as just explained. We then have $`S^2=(x_1x_1^{})^2`$ and $`d𝐱=dx_1`$. After some algebra following the manipulations of Ref  for the straight line we have

$$I_{straightline}=2dx_1dx_1^{}\frac{\delta ^2}{(x_{1}^{}{}_{}{}^{2}+\delta ^2)(x_{1}^{}{}_{}{}^{2}+\delta ^2)},$$
(25)
which is in fact non-singular at $`𝐱𝐱^{}`$. The absence of the singularity may be traced to the fact that $`I`$ is symmetric in $`x,x^{}`$ and so must contain only even powers in $`𝐱𝐱^{}`$. The potential singularity from Eq , however would be odd in this variable, and so is in fact absent.

Concerning the possible singularity at $`x^2=0`$, as occurs for a curve through $`𝒪`$ , first consider $`I(0,x^{})`$, i.e. one variable at some ordinary point and the other near zero. In this case the first and third term of $`I`$, Eq  cancel since $`d_x^{}lnS^2=d_x^{}lnx^2`$ for $`x=0`$. On the other hand the second term of $`I`$ is non-singular, so there is no singularity.

As for $`I(0,0)`$, we examine the straight line through the origin. By direct inspection of $`I`$ one finds

$$I_{straightlinethrough𝒪}=0.$$
(26)
Therefore there are no singularities connected with $`x=0`$.

Eq  may be viewed as the consequence of a general symmetry property since two points along a ray may be interchanged by an inversion. Therefore by Eq  $`I(x,x^{})=I(x^{},x)`$. On the other hand, by Eq  we also have $`I(x,x^{})=+I(x^{},x)`$, so $`I`$ is zero.

Since in Eq  we have zero for the straight line through $`𝒪`$ , we look at the simplest curvature or “acceleration” contribution, that of the circle. Consider a circle passing through $`𝒪`$ and use the relation for the length of a chord $`l=2Rsin\varphi `$, where $`\varphi `$ is the half-angle subtended at the center of the circle. This leads to

$$I=1/2\left[dlnsin^2\varphi dlnsin^2\varphi ^{}+dlnsin^2(\varphi \varphi ^{})dlnsin^2\varphi ^{}+dlnsin^2(\varphi \varphi ^{})dlnsin^2\varphi \right]$$
(27)
that is $`I=2d\varphi d\varphi ^{}\left[cot\varphi cot\varphi ^{}+cot(\varphi \varphi ^{})(cot\varphi ^{}cot\varphi )\right]`$. Using the identity $`cot(\varphi \varphi ^{})=\frac{cot\varphi cot\varphi ^{}+1}{cot\varphi ^{}cot\varphi }`$ we have finally

$$I_{circlethrough𝒪}=2d\varphi d\varphi ^{},$$
(28)
where $`0\varphi \pi `$.

If we have only a segment of a circle passing through $`𝒪`$ the integration will only be over the corresponding angle, and Eq  may be said to approach Eq  in the sense that as the radius of the circle $`R`$ becomes very large this angular segment becomes very small. The very simple form of Eq  suggests that it may sometimes be preferable to use coordinates not centered on $`𝒪`$ but rather on the center of curvature for the curve at $`𝒪`$ .

## IX simplest or “reference” cases

There are four configurations to discuss, according to whether the curve is closed or infinite and whether $`𝒪`$ is on or off the curve. For orientation we discuss a simplest or “reference” curve for each case.

### A circle, not through $`𝒪`$

This is the simplest case of the finite closed curve where the center of inversion is not on the curve. Inversion of a circle produces another circle, and similarly inversion of a general closed curve will produce another closed curve, as long as $`𝒪`$ is not on the curve.

For the circle one finds $`n`$ by directly integrating Eq

$$n_{circle}=2\pi ^2.$$
(29)
Since $`n`$ is purely a property of the geometric figure, this holds for any circle.

We therefore conclude from Eq

$$I_{circlenotthrough𝒪}=0.$$
(30)
As on check on this argument we can show directly that the integral of $`I`$ is zero by Eq , that is by Eq , for those circles going identically into themselves under inversion. This occurs for inversion through the center of a circle for example, or when $`𝒪`$ is in the plane but outside the circle, take $`a^2=dD`$ in Eq  with $`d`$ the closest point of the circle to $`𝒪`$ and $`D`$ the furthest.

### B  infinite straight line not through $`𝒪`$ , or circle through $`𝒪`$

We now turn to the consideration of curves involving infinities. The simplest case is the infinite straight line, not through $`𝒪`$ . Under inversion this becomes a circle through $`𝒪`$ . The general case here refers to curves, which although they become the same straight line  at large distances, have some arbitrary form at finite distances (Fig 3). As with the straight line/circle pair, under inversion these become a closed curve passing through $`𝒪`$ ; and vice-versa. Since they are the images of each other under inversion, we consider the two cases together.

That some subtlety is involved is evident from the fact that on the one hand we have

$$n_{straightline}=0.$$
(31)
But for the inversion of the straight line, a circle crossing the origin (or for that matter any circle), we have not zero but rather

$$n_{circle}=2\pi ^2,$$
(32)
as is found by integrating the definition Eq .

That different values for $`n`$ result is perhaps not entirely surprising since the inversion has produced a basic change in the figure, an infinite curve becoming a finite curve. In Ref  this difference was associated with the “anomaly”.

According to Eq , when starting from the straight line we must integrate $`I`$ over a circle through $`𝒪`$ . Or when beginning with a circle we integrate $`I`$ along an infinite straight line. To be consistent with Eq  and Eq , these integrations ought to produce non-zero and opposite sign contributions.

Indeed, in Eq  we already have $`I`$ for the straight line. Carrying out the integral we see that it is independent of $`\delta `$ and gives

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}I_{straightlinenotthrough𝒪}=2\pi ^2.$$
(33)
Similarly, we use Eq  to integrate over a circle through $`𝒪`$

$$I_{circlethrough𝒪}=2\pi ^2.$$
(34)
These results are in agreement with Eq, Eq  and Eq  and of course Ref . We note the contrast of Eq  with Eq ; or between Eq  and Eq  below. Apparently while $`I`$ is a non-singular object along a given curve, an infinitesimal change in that curve can produce a finite effect in the integral.

### C  straight line through $`𝒪`$

For closed finite curves, as typified by the circle not through $`𝒪`$ , we had neither points at infinity nor points at zero. In the previous subsection we had either one or the other, mapping into each other under inversion. Finally, we consider having both at the same time: curves which are both infinite and pass through $`𝒪`$ . The simplest or reference case is the straight line through $`𝒪`$ . An inversion returns the same straight line through the $`𝒪`$ . Observe that also for the general case of this class, namely an arbitrary curve passing through $`𝒪`$ and becoming the same straight line at $`\pm \mathrm{}`$, inversion produces a curve of this same type (Fig 2).

In Eq  we had $`I_{straightlinethrough𝒪}=0`$. Evidently then

$$I_{straightlinethrough𝒪}=0.$$
(35)
This is of course in agreement with the arguments of section VI since under inversion the straight line through $`𝒪`$ , being on a ray, goes identically into itself.

This case is most liable to be of interest physically in minkowski space where it will correspond to the path of a charge, and the inversion also produces a possible path for a charge.

## X once-integrated expression

We now turn to the explicit integration of $`I`$. To integrate $`I`$, note that it consists of product pairs where one member of each product does not contain both variables. Hence one member of each pair, given as a total derivative, can be explicitly integrated. The limits of integration are also the same in both variables. Therefore it appears useful to interchange the names of $`x,x^{}`$ say in the last term of $`I`$, Eq  to arrive at the expression $`(1/2)d_x[lnx^2+2ln(𝐱𝐱^{})^2]𝑑lnx^2`$, where the $`x`$ integration can be done immediately. However, this has the disadvantage that it is no longer explicitly symmetric in $`x,x^{}`$, which introduces a singularity of the integrand for $`𝐱𝐱^{}`$. Although this is only apparent and integrates to zero, it is perhaps more comfortable to let $`ln[S^2]ln[S^2+\eta ^2]`$ so as to make the integrand explicitly finite and then let $`\eta 0`$ at the end of the calculation. It may be verified that $`I`$ is completely regular under this procedure. After the relabeling of variables we have an explicitly finite expression (for $`𝒪`$ not on the curve) and can now perform the $`x`$ integration along an arbitrary curve between limits $`𝐀`$ and $`𝐁`$ to obtain

$$\underset{𝐁}{\overset{𝐀}{}}I=_𝐁^𝐀𝑑lnx^2\left(\frac{1}{2}ln\frac{A^2}{B^2}+ln\frac{(𝐱𝐀)^2}{(𝐱𝐁)^2}\right)\text{C}notthrough𝒪.$$
(36)
Observe that there are no singularities connected with this expression. At an endpoint, for $`𝐱𝐀`$ we have $`^𝐀𝑑lnx^2ln(𝐱𝐀)^2^1𝑑x(lnA^2+ln(x1))`$ which is non-singular; while $`x0`$ is excluded by assumption.

However, we shall also need a once-integrated expression when $`𝒪`$ is on the curve, where $`x0`$ must be considered. To see what replaces Eq  note that Eq  or Eq  tell us that there is in fact no singularity associated with $`𝒪`$ . Therefore we can excise an infinitesimal region along the curve around $`x=0`$ without affecting the value of the integral:

$$\underset{𝐁}{\overset{𝐀}{}}I=\left(_𝐁^ϵ+_ϵ^𝐀\right)\left(_𝐁^ϵ+_ϵ^𝐀\right)I,ϵ0$$
(37)
Now again writing $`I=(1/2)d_x[lnx^2+2ln(𝐱𝐱^{})^2]𝑑lnx^2`$ and carrying out the $`d_x`$ we obtain

$$\underset{𝐁}{\overset{𝐀}{}}I_{Cthrough𝒪}=\left(_𝐁^ϵ+_ϵ^𝐀\right)dlnx^2\left(\frac{1}{2}ln\frac{A^2}{B^2}+ln\frac{(𝐱𝐀)^2}{(𝐱𝐁)^2}+ln\frac{(𝐱+ϵ)^2}{(𝐱ϵ)^2}\right).$$
(38)
We first consider the last term, $`_𝐁^ϵ+_ϵ^𝐀𝑑lnx^2ln\frac{(𝐱+ϵ)^2}{(𝐱ϵ)^2}`$. As $`ϵ0`$, the integrand vanishes for any non-infinitesimal value of $`x`$. Thus it suffices to evaluate the integral for a straight line in the vicinity of $`𝒪`$ . Introducing the variable $`y=x/ϵ`$, this leads to the integral $`_{\mathrm{}}^1+_1^{\mathrm{}}\frac{2dy}{y}ln(\frac{y+1}{y1})^2=8_1^{\mathrm{}}\frac{dy}{y}ln(\frac{y+1}{y1})=2\pi ^2`$ , so we may write

$$_𝐁^ϵ+_ϵ^𝐀𝑑lnx^2ln\frac{(𝐱+ϵ)^2}{(𝐱ϵ)^2}=2\pi ^2ϵ0,$$
(39)
for any smooth curve through $`𝒪`$ , with A and B on opposite sides of $`𝒪`$ (compare Eq \[2.15\] of Ref).

As for the remaining part of Eq , we verify that it is non-singular. For the endpoints, say $`𝐱𝐀`$, it is finite for the same reason given above with regard to Eq . For $`x0`$, we write $`ln\frac{(𝐱𝐀)^2}{(𝐱𝐁)^2}=ln\frac{A^2}{B^2}+ln(12\frac{\mathrm{𝐀𝐱}}{A^2}+\frac{x^2}{A^2})ln(12\frac{\mathrm{𝐁𝐱}}{B^2}+\frac{x^2}{B^2})`$. The constant terms $`_𝐁^ϵ+_ϵ^𝐀(1/2)ln\frac{A^2}{B^2}=(1/2)ln\frac{A^2}{B^2}(ln\frac{ϵ^2}{B^2}+ln\frac{A^2}{ϵ^2})=(1/2)(ln\frac{A^2}{B^2})^2`$ are nonsingular. For the $`x`$-dependent terms we can write $`ln(12\frac{\mathrm{𝐀𝐱}}{A^2}+\frac{x^2}{A^2})2\frac{\mathrm{𝐀𝐱}}{A^2}`$, and similarly for the $`B`$ term, leading to an integral of the type $`x𝑑lnx`$ which is also non-singular at $`x=0`$.

Thus the expression is $`ϵ`$ independent and well defined, and we introduce the symbol $`𝒫`$ for this principal value-like integral: $`𝒫_𝐁^𝐀=_𝐁^ϵ+_ϵ^𝐀`$ and can thus finally write for Eq

$$_𝐁^𝐀I=2\pi ^2+𝒫_𝐁^𝐀𝑑lnx^2\left(\frac{1}{2}ln\frac{A^2}{B^2}+ln\frac{(𝐱𝐀)^2}{(𝐱𝐁)^2}\right)Cthrough𝒪.$$
(40)
We wish to use Eq  and Eq  in the next section to extend the simple results found for straight lines and circles to general curves, but first we check that these formulas give the expected results in these simple cases:

1)For the circle, not through $`𝒪`$ , we expect zero, according to Eq , which is indeed what results from setting $`𝐀=𝐁`$ in Eq .

2)For the infinite straight line, not through $`𝒪`$ , we expect $`2\pi ^2`$ according to Eq . Using Eq , we consider the straight line at a distance $`\delta `$ from $`𝒪`$ , as in Eq  and introduce the variable $`y=x/\sqrt{AB}`$. In the limit $`A,B\mathrm{}`$ such that $`A/B1`$ and $`\delta /\sqrt{AB}0`$ one obtains the integral $`4_1^1\frac{dy}{y}ln\frac{1y}{1+y}`$, which is indeed $`2\pi ^2`$.

3)Turning now to $`𝒪`$ on the curve and Eq , for the circle we expect $`2\pi ^2`$, according to Eq . This is what we obtain upon setting $`𝐀=𝐁`$ in Eq .

4)For the final example, the straight line through $`𝒪`$ , we expect zero according to Eq . To evaluate the principal value integral, we repeat the arguments just given for the infinite straight line not through $`𝒪`$ , with the difference that in place of $`\delta `$, we now have $`ϵ/\sqrt{AB}0`$. This leads to $`8_0^1\frac{dy}{y}ln\frac{1y}{1+y}=2\pi ^2`$ again and Eq  is zero.

We emphasize that the complications of this section arise from our desire to bring $`I`$ into the once-integrated form and the resulting asymmetric treatment of the variables; $`I`$ itself is perfectly well behaved for a given curve. Perhaps if another method could be found for the problem of general curves as discussed in the next section these complications could be avoided.

## XI Integral of $`I`$ as an invariant

Perhaps the most remarkable property of $`I`$ is that it is a type of invariant, having the same value for all curves of a given configuration. In our language, this was what was essentially concluded in Ref  for certain field theoretic amplitudes (for our cases 1-3). We shall show this using Eq  and Eq  of the previous section. Since the quantity in parenthesis in Eq  and Eq  does not depend on the particular curve, the nature of the actual curve in question appears only in the remaining single integration over $`x`$. This feature that makes the once-integrated expressions useful for the examination of general curves.

### A  closed curves not through $`𝒪`$

For our first case we take the generalization for section IX A : inversion of a general, finite closed euclidean curve, with $`𝒪`$ not on the curve. The inversion operation produces another finite closed curve. Therefore we can set $`𝐀=𝐁`$ and the integrand of Eq  is zero.

$$I_{closednotthrough𝒪}=0.$$
(41)
This is naturally in agreement with the example Eq , and in view of Eq  we can finally conclude that for this case $`n`$ is inversion invariant:

$$n(C)=n(C_{inv})$$
(42)
where $`C_{inv}`$ is the inversion of any finite closed euclidean curve $`C`$, with the center of inversion not on the curve $`C`$.

### B  infinite curves not through $`𝒪`$

We now consider the generalization of the infinite straight line. By these we mean curves which although they become the (same) straight line at large distances, have some arbitrary form at finite distances (Fig 3).

We employ the once-integrated Eq  again, where the limits A and B limits are to be sent to $`\pm \mathrm{}`$. The quantity in parenthesis in Eq  would be the same for any curve and in particular for the straight line. Since the curve under consideration only differs from the straight line in a finite region, say between points $`𝐅`$ and $`𝐆`$, let us add and subtract a straight line contribution from $`𝐅`$ to $`𝐆`$ in Eq . This leads to

$$_𝐁^𝐀𝑑lnx^2(\mathrm{})=_{straightline}𝑑lnx^2(\mathrm{})+_{finite}𝑑lnx^2(\mathrm{}).$$
(43)
The $`(\mathrm{})`$ stands for the parenthesis in Eq . $`_{straightline}`$ is from B to A, and as A, B are sent to $`\pm \mathrm{}`$ it becomes the integral for the infinite straight line, $`2\pi ^2`$. The second term, $`_{finite}`$, stands for the integral from $`𝐅`$ to $`𝐆`$ along the actual curve in question minus the integral $`𝐅`$ to $`𝐆`$ along a straight line.

We now argue that $`_{finite}`$ goes to zero as $`𝐀,𝐁\pm \mathrm{}`$. This is because $`x`$ in $`_{finite}`$ is confined to finite values as say $`𝐀\mathrm{}`$. Thus we can write $`ln\left[(𝐱𝐀)^2\right]=lnA^2+ln\left[12\frac{\mathrm{𝐀𝐱}}{A^2}+\frac{x^2}{A^2}\right]lnA^22\frac{\mathrm{𝐀𝐱}}{A^2}`$ to leading order in $`1/A`$. The $`2\frac{\mathrm{𝐀𝐱}}{A^2}`$ will lead to a contribution to the integral vanishing as $`1/A`$. Then $`_{finite}`$ becomes the difference of the integral of a total derivative over two paths between the same endpoints and so is zero.

( There might appear to be some difficulty with the argument when $`𝒪`$ is located such that the straight line from $`𝐅`$ to $`𝐆`$ passes through it, i.e when the straight lines at infinity lie along a ray from $`𝒪`$ . We can deal with this by replacing the straight line as the reference curve by a curve where a semi-circle avoids $`𝒪`$ . In this case the semi- circle contribution goes to zero as $`𝐀\mathrm{}`$, while the straight line integrals can be evaluated in this limit to again give $`2\pi ^2`$.)

We thus conclude that the integral $`I`$ over any curve, not passing through $`𝒪`$ , differing from the straight line in a finite region is the same as that for the simple straight line, not passing through $`𝒪`$ :

$$I_{infinitenotthrough𝒪}=2\pi ^2.$$
(44)
### C closed curves through $`𝒪`$

Inversion now leads to the generalization of the circle through $`𝒪`$ , the arbitrary closed curve through $`𝒪`$ . We use Eq , and setting A= B, $`I`$ is the same as for the circle:

$$I_{closedthrough𝒪}=2\pi ^2.$$
(45)
And for $`n`$ we can say

$$n(C)=n(C_{inv})+2\pi ^2.$$
(46)
Where $`C`$ is a generalized circle and $`C_{inv}`$ is its inversion, a generalized infinite straight line.

### D infinite curve through $`𝒪`$

We come to our last case, the generalization of the straight line through $`𝒪`$ . By this we mean an infinite curve with our usual condition that it becomes the same straight line at infinity, but now also passing through $`𝒪`$ at some finite point (Fig 2). This class inverts into itself.

We now use Eq  where we must evaluate $`𝒫_𝐁^𝐀𝑑lnx^2(\mathrm{})`$. In this integral, adding and subtracting a straight line piece between the finite points $`𝐅`$ and $`𝐆`$ (As in Fig 3 but with $`𝒪`$ on the solid curve), we have again

$$𝒫_𝐁^𝐀𝑑lnx^2(\mathrm{})=_{straightline}𝑑lnx^2(\mathrm{})+𝒫_{finite}𝑑lnx^2(\mathrm{}),$$
(47)
where we have dropped the $`𝒫`$ in $`_{straightline}`$ since it does not pass through $`𝒪`$ . Since $`_{straightline}=2\pi ^2`$ as $`𝐀,𝐁\mathrm{}`$, the $`2\pi ^2`$ from Eq  is cancelled and we are left with

$$𝒫_{finite}𝑑lnx^2(\mathrm{})=\left(_𝐅^𝐆+𝒫_𝐅^𝐆\right)dlnx^2\left(\frac{1}{2}ln\frac{A^2}{B^2}+ln\frac{(𝐱𝐀)^2}{(𝐱𝐁)^2}\right),$$
(48)
where the first integral is along the straight line piece and the second integral is along the curve in question.

As before, for large $`𝐀,𝐁`$ the $`(\mathrm{})`$ goes to an $`x`$ independent piece and terms vanishing as $`𝐀,𝐁\mathrm{}`$. Since $`_𝐅^𝐆𝑑lnx^2=𝒫_𝐅^𝐆𝑑lnx^2`$ for the two different paths, Eq  goes to zero and we can conclude

$$I_{infinitethrough𝒪}=0$$
(49)
for arbitrary curves. This implies finally for $`n`$

$$n(C)=n(C_{inv})$$
(50)
where $`C_{inv}`$ is the inversion of an infinite euclidean curve $`C`$, with the center of inversion on the curve.

### E summary

We can summarize the results of this section as follows. We have two types of curves: “finite” curves and “infinite” curves. The first are closed curves, the second open curves becoming the same straight line at large distances. There are also two possibilities for the placement of the origin of inversion $`𝒪`$ : “on” and “off” the curve. Labeling the configurations from 1 to 4 we can exhibit their properties in a table:

| CONFIGURATION | CURVE | $`𝒪`$ | RESULT | $`I`$ |
| --- | --- | --- | --- | --- |
| 1 | finite | off | 1 | $`0`$ |
| 2 | finite | on | 3 | $`2\pi ^2`$ |
| 3 | infinite | off | 2 | -$`2\pi ^2`$ |
| 4 | infinite | on | 4 | 0 |

The column “Result” refers to the configuration resulting from the inversion. Thus: configuration 2 is a closed finite curve with $`𝒪`$ on the curve. The integral of $`I`$ over this curve has the value $`2\pi ^2`$. The inversion produces configuration 3, which is an infinite curve with $`𝒪`$ off the curve.

A notable feature of the table is that those configurations, namely 1 and 4, which invert into themselves have $`I=0`$. This is a type of generalization of Eq , and appears to answer the question raised in sect VI as to the wider meaning of “$``$”. Apparently two curves should be considered “equivalent” when they belong to the same configuration in the sense of the table. Depending on the position of $`𝒪`$ , a curve may map into its own configuration or not. When it does, the “self conjugate” property $`I=0`$ obtains.

## XII  infrared/ultraviolet duality

We remarked in  that it appeared as if the finiteness of $`n`$ both at short distances (ultraviolet) and at long distances (infrared) could in a sense be attributed to the same thing, namely that the average velocity $`U(x,x^{})`$ and the instantaneous velocity $`u`$ become equal. In euclidean space $`U`$ and $`u`$ refer to the chord and tangent of the curve respectively. At short distances as $`xx^{}`$we have $`Uu`$ because the curve is taken to be smooth, so for small enough intervals the chord and the tangent become the same. At large distances the curve has by assumption a constant and equal slope at $`\pm \mathrm{}`$, so that $`Uu`$ again.

It is interesting to note how the configurations of the table where $`𝒪`$ is “on”, that is where “Result” is an infinite curve, represents just this situation. In inversion a point on the curve is projected along the (directed) ray connecting it to the origin. The point goes to large distances if it was close to the origin and goes to small distances if it was far from the origin. If $`𝒪`$ is directly on the curve, the points approaching the origin along the curve from one side will be sent to a straight line at $`+\mathrm{}`$, while points approaching $`𝒪`$ from the other side will be sent to a straight line at $`\mathrm{}`$ in the opposite direction. That the same, straight line at $`\pm \mathrm{}`$ results is a consequence of the presumed smoothness of the curve at $`𝒪`$ .

Similarly, a curve coming from large distances and finally going to large distances as the same straight line will be mapped into a curve smooth at the origin. If the infinite slopes had been different there would be a kink at the origin.

Inversion thus makes it clear how our two assumptions needed to make $`n`$ finite both in the infrared and the ultraviolet are related: under an inversion whose origin is on the curve, curves which are locally smooth become straight and parallel at $`\mathrm{}`$ and vice-versa.

## XIII minkowski space

In minkowksi space a new aspect enters in that Eq  can make points that are time-like separated into points that are space-like separated. That is, in minkowski space, where $`S^2=(x_0x_0^{})^2(𝐱𝐱^{})^2`$ the purely algebraic relation Eq  is still true

$$S^2\frac{a^2}{x^2}\frac{a^2}{x^2}S^2,$$
(51)
but for an initially physical time-like curve, with $`𝒪`$ “off” there will be points on it space-like to $`𝒪`$ . This gives $`x^2`$ negative in Eq  and a time like-like separation $`S^2>0`$ may be mapped into a space-like separation $`S^2<0`$ .

On the other hand for a time-like path with $`𝒪`$ on the curve another possible physical path results. For $`𝒪`$ on the path, $`x^2`$ is always positive and $`S^2`$ in Eq  cannot change sign and all relatively time-like pairs of points remain time-like. This operation is like configuration 4 of the euclidean problem where, according to the table, $`n`$ is inversion invariant. It might be interesting to investigate if the inversion method could be helpful in analyzing certain practical radiation problems.

## Acknowledgments

I am grateful to David Gross for stressing the interest of studying the inversion, for several discussions, and for bringing the references  and  to my attention. I would also like to thank E. De Rafael for conversations concerning dilogarithms and related problems, as well as E. Seiler for discussions and a reading of the manuscript.


format: text

text:
# Field-doping of C60 crystals: Polarization and Stark splitting

## I Motivation

The proposal of doping C60 crystals in a field-effect device (FET) and the possibility of metallic conduction and even superconductivity in such devices had raised wide-spread interest. While the revelation of dishonest data handling in some casesbeasley  led to a severe damping of the initial enthusiasm, fundamental aspects of field effect doping remains a timely and interesting problem. For field-effect transistors made with self-assembled monolayers this question was addressed in Ref. SAMFET,  and for the reported enhancement of the superconducting transition temperature in C60 crystals intercalated with haloform molecules in Refs. haloform, . Attempts to observe the field-effect in graphite were reported in Ref. graphite, . Here we address, from a theoretical point of view, the question how strongly C60 can be doped in an electric field before its electronic structure is substantially changed and how this structure changes in even stronger fields. This is relevant for understanding the fundamental features of field-effect devices based on C60 and involves a number of interesting physical problems.

It appears that doping C60 crystals in a field-effect device (FET) is very hard to achieve in practice, one of the reasons being the exceptionally strong fields that are required. Very strong fields, however, not only induce charge carriers, but also polarize the molecules, and, due to the Stark effect, in general, lift degeneracies. These effects are of particular importance, as C60 is quite polarizable, and as its molecular levels are highly degenerate. The term “field-doping” naively implies that these effects are small, such that the main effect of the external field is inducing charge carriers into electronic levels which are essentially unaffected by the field. It is clear that if the external field is strong enough, the electronic structure of the crystal will be strongly modified by the field, and one thus can no longer speak of doping. A fundamental question about field-effect devices is therefore connected with the doping levels achievable, before the electronic structure of the active material is substantially altered. Moreover we discuss in detail how the electronic structure is changed, when the external field is beyond this “doping limit”. In this regime the commonly used analogy of the field-doped C60 with the alkali-doped fullerenes no longer holds. Nevertheless a field-effect device with high carrier concentration would be an interesting device in its own right. One might, e.g., speculate that for C60 in a strong electric field the electron-phonon coupling is enhanced compared to unperturbed C60, as for molecular orbitals of lower symmetry less couplings are forbidden by symmetry.

For a first estimate of the fields involved in field-doping, consider a simple capacitor. Given a charge per area $`\sigma `$ on the plates, the electric field between the plates is $`E=4\pi \sigma `$. Assuming that in field-doped C60 the induced charge resides in the top-most layer,art:wehrli1  inducing $`n`$ elementary charges per molecules requires an external field (originating from the gate electrode) of $`E_{\mathrm{ext}}=2\pi n/A_{\mathrm{mol}}`$, where $`A_{\mathrm{mol}}`$ is the area per molecule in the top-most layer of the crystal. For a C60 crystal with lattice constant $`a`$ 14 Å typical areas per molecule are $`A_{(111)}=\sqrt{3}a^2/4`$ for the (111)-plane and $`A_{(001)}=a^2/2`$ for the (001)-plane. Even though these areas are quite large, the external fields necessary for field doping are substantial, being of the order of 1 V/Å per induced elementary charge per molecule. This is, however, not the field experienced by a molecule. As C60 is highly polarizable ($`\alpha 83`$ Å3), in the solid, the external field at the site of a molecule is screened by the polarization of the neighboring molecules: A monolayer of dipoles $`p`$ centered on the lattice sites $`R_i`$, generates a field $`E_{\mathrm{dip}}=p_{i0}|R_iR_0|^3`$ at $`R_0`$, where the sum is over all sites in the monolayer, except $`R_0`$. For the (111) layer, this sum is about $`31.2/a^3`$, for the less dense (001) layer it is about $`25.6/a^3`$. As the dipole moments $`p`$ of the molecules are induced by the screened field $`E_{\mathrm{scr}}=E_{\mathrm{ext}}+E_{\mathrm{dip}}`$ at the site of the molecule ($`p=\alpha E_{\mathrm{scr}}`$), we find, by self-consistent solution, that the external field is reduced by about a factor of two.

Inducing charge and polarizing the molecules is, however, not the only effect of the external field. In addition it also leads to a splitting of the molecular levels — the Stark effect. As the molecular orbitals of C60 have a definite parity, a homogeneous field splits the levels only in second order. Thus for low fields the splitting is quite small, but increases quickly for larger fields. We can expect that the splitting of the molecular levels disrupts the electronic structure of the crystal only when it becomes of the order of the band width, which is about $`1/2`$ eV in C60. Calculations indicate that the splittings of the $`T_{1u}`$\- and $`H_u`$-levels in a homogeneous field are surprisingly small, being less than $`1/2`$ eV up to $`E_{\mathrm{scr}}2`$ V/Å (cf. figure 2). Given the crudeness of the argument, it thus seems entirely feasible that doping of a few elementary charges per molecule could be achieved before the electronic structure of the C60 is substantially changed.

In the argument above we have described the C60 molecules as polarizable points. It is then natural to refine the model by considering also higher multipoles. Such a multipole expansion is particularly suitable for C60 as the molecules are nearly spherical. The approach is then the following: First we calculate the response of a C60 molecule to external multipole fields. Then we self-consistently solve the electrostatic problem for a lattice of molecules in a homogeneous external field. This provides us with the multipole expansion of the screened field acting on each molecule in the solid. Given that screened field, we can then determine the splitting of the molecular levels as a function of the induced charge.

The organization of the paper reflects this approach: In section II we describe the density functional calculations for determining the multipole response and the Stark splitting for a C60 molecule in a multipole field for several symmetrical configuration. Using group theory, in section III, the irreducible parameters for the multipole response are determined. This allows the calculation of the polarization for arbitrary configurations. In section IV we give an analogous treatment for the Stark splittings and explicitly show how the splitting changes as the molecule is rotated relative to the external field. In section V we use these ingredients to self-consistently solve for the screened field seen by a molecule in a charged monolayer. The splitting of the molecular levels in this self-consistent multipole field and the effect of this splitting on the density of states is presented in section VI. Our conclusions are presented in section VII. The methods for calculating the response and splitting for an arbitrarily oriented external field from the results of the density functional calculations that were performed only for special orientations are described in the appendices. Appendix A gives an example of how to calculate the response of a molecule using the irreducible parameters derived in section III. In appendix B we derive the coupling matrices needed for the calculating the level splitting when the molecule is rotated in the external field. Finally, appendix C gives the derivation of the matrix describing the field generated by a lattice of identical multipoles at the origin, which is needed for finding the self-consistent electrostatic field.

## II Response of a C60 molecule

To determine the response of a C60 molecule to external multipole fields we have performed all-electron density functional calculations using Gaussian-orbitals.NRLMOL  The basis set comprises 5s4p3d for carbonbasis  and we use the Perdew-Burke-Ernzerhof functional.PBE

In our calculations we apply an external multipole field and study the change in the multipole moments of the charge density and the splitting of the molecular levels as a function of the strength of the external field. To take advantage of the molecule’s symmetry we consider multipole fields with the $`z`$-axis along the 2-, 3-, and 5-fold axis of the molecule (cf. figure 1). As these axes are each contained in a mirror plane, which we chose to be the $`x`$-$`z`$-plane, we can treat the thus oriented molecule as having symmetry group $`C_{2v}`$, $`C_{3v}`$, and $`C_{5v}`$ respectively. Applying external multipole fields with $`l>0`$, this symmetry is maintained if the fields are proportional to the real part of the spherical harmonic $`Y_{lm}`$, where $`m`$ is an integer multiple of the order ($`n=`$2, 3, or 5) of the symmetry axis. Likewise, the response of the charge density will only have multipole components proportional to $`\text{Re}(Y_{lm})`$, with $`m`$ an integer multiple of $`n`$. Calculations were done for such symmetry conserving multipole fields up to $`l=6`$. For the 3-fold axis oriented along $`z`$, we have, in addition, calculated the response to external fields proportional to $`\text{Im}(Y_{5,3})`$, $`\text{Im}(Y_{6,3})`$, and $`\text{Im}(Y_{6,6})`$, i.e., with a symmetry lowered to $`C_3`$. As we are interested in the linear response, we have considered small multipole fields and made sure that the calculated response of the charge density is indeed proportional to the strength of the external field. We find that the linear response of a C60 molecule is very similar to that of a metallic sphere of radius 4.4 Å. This effective sphere radius shows a slight increase with $`l`$. In addition there are weak off-diagonal terms. To judge the accuracy of our calculation, we have checked how well the selection rules, that are not already imposed by the $`C_{nv}`$ symmetry, are fulfilled for these off-diagonal terms. From the calculated response, we have determined the irreducible linear response coefficients, which will be given in the next section.

In addition to polarizing the molecule, the external field also splits the degenerate molecular levels of the C60 molecule. As the unperturbed molecular orbitals have a definite parity, for a multipole field with odd $`l`$ there will be no first-order splitting — the quadratic Stark effect. On the other hand, a multipole field with even $`l`$ can couple states of like parity, so in that case the splitting is linear. This is shown in figure 2. In section IV it will be demonstrated how the splitting of the HOMO- and LUMO-orbitals that were calculated for high-symmetry geometries can be extended by group theory to arbitrary orientations of the molecule relative to the external multipole field.

All calculations have been performed for the equilibrium geometry of the unperturbed, neutral C60 molecule. To estimate the effect of an external field on the shape of the molecule, we have relaxed the structure in the presence of homogeneous external fields of up to 1 V/Å. We find only small changes (up to about 0.005 Å) in the lengths of the bonds (1.40 Å for the short and 1.45 Å for the long bonds). Likewise, the polarizability changes by less than 1.5%.

Finally, we have calculated the total energy of the isolated C60 ion as a function of its charge $`q`$ (spin unpolarized calculation with relaxed geometries) and extracted the second order term $`1/2U_0q^2`$. To compare with previous calculations,Mark  we find for the polarizability (multipole field with $`l=1`$) of the neutral molecule $`\alpha =`$9.3 bohr3/atom, and $`U_0=`$3.2 eV.

## III Irreducible response parameters

The polarizability $`\alpha `$ of a molecule describes the linear dependence of the induced dipole moment $`𝐩=\alpha 𝐄`$ on the applied electric field $`𝐄`$. For a multipole expansion, $`\alpha `$ becomes a matrix $`\alpha _{l_1m_1l_2m_2}`$ describing the response to all multipole fields.book:Stone  To fix the notation (which follows Ref. book:Stone, ) we briefly review the definition of the multipole response matrix.

The solutions of the Laplace equation $`^2V(𝐫)=0`$ are given by

$$V(𝐫)=V^\text{e}+V^\text{i}=\underset{lm}{}V_{lm}R_{lm}^{}(𝐫)+Q_{lm}I_{lm}^{}(𝐫),$$
(1)
where the two terms denote the external potential ($`V^\text{e}`$), and the induced potential ($`V^\text{i}`$) due to a charge distribution $`\rho `$ located around $`𝐫=0`$. Note that both, the Laplace equation as well as the expansion of $`V^\text{i}(𝐫)`$ into multipoles only holds for $`𝐫`$ which lie outside the charge distribution. We have introduced the regular and irregular spherical harmonics,book:Stone

$`R_{lm}(𝐫)`$ $`=`$ $`r^l\sqrt{{\displaystyle \frac{4\pi }{2l+1}}}Y_{lm}(\mathrm{\Omega }),`$ (2)
$`I_{lm}(𝐫)`$ $`=`$ $`{\displaystyle \frac{1}{r^{l+1}}}\sqrt{{\displaystyle \frac{4\pi }{2l+1}}}Y_{lm}(\mathrm{\Omega }).`$ (3)
Special cases for the regular spherical harmonics are $`R_{00}=1`$ and $`R_{10}=z`$, hence, the external field, $`V_{00}=V_0`$ corresponds to a constant shift, and $`V_{10}=E_z`$ is the $`z`$-component of the electric field. For the irregular spherical harmonics we have $`I_{00}=1/r`$ and $`I_{10}=z/r^3`$, thus for the induced potential, $`Q_{00}=q`$ gives the monopole charge while $`Q_{10}=p_z`$ is the dipole moment. Generally, the coefficients $`Q_{lm}`$ are the multipole moments of the charge distribution $`\rho `$

$`Q_{lm}`$ $`=`$ $`{\displaystyle d^3r\rho (𝐫)R_{lm}(𝐫)}.`$ (4)
Decomposing the charge distribution $`\rho =\rho _0+\mathrm{\Delta }\rho `$ into the unperturbed charge density and the change in the charge density due to the external potential, we obtain a decomposition of the multipole moments $`Q_{lm}=Q_{lm}^0+\mathrm{\Delta }Q_{lm}`$. Within linear response, the coefficients $`\mathrm{\Delta }Q_{lm}`$ of the induced multipole moments depend linearly on the coefficients $`V_{lm}`$ of the external potential, which defines the linear-response matrix $`\alpha _{l_1m_1l_2m_2}`$:

$$\mathrm{\Delta }Q_{l_1m_1}=\underset{l_2m_2}{}\alpha _{l_1m_1l_2m_2}V_{l_2m_2},$$
(5)
where the sign takes into account that the induced fields oppose the external fields. Then $`\alpha _{1m_1\mathrm{\hspace{0.17em}1}m_2}`$ gives the dipolar response tensor, while $`\alpha _{\mathrm{00\hspace{0.17em}00}}`$ is the self-capacitance $`U_0`$.

The interaction energy of the molecule with the external potential is $`E_{\text{ext}}=d^3r\rho (𝐫)V^\text{e}(𝐫)`$, which, using the previous definitions, reduces to

$$E_{\text{ext}}=\underset{lm}{}V_{lm}Q_{lm}^{}.$$
(6)
Therefore, $`V_{lm}`$ and $`Q_{lm}`$ are pairs of conjugate variables and the total energy of the molecule as a function of the external field is given by

$$E_{\text{tot}}=\frac{1}{2}\underset{\genfrac{}{}{0pt}{}{l_1m_1}{l_2m_2}}{}\alpha _{l_1m_1l_2m_2}V_{l_1m_1}V_{l_2m_2}^{}+\underset{lm}{}V_{lm}(Q_{lm}^0)^{}.$$
(7)
Since this is a quadratic form, we see that the matrix $`\alpha `$ is hermitian. We can make $`\alpha `$ real and symmetric by unitarily transforming to a real basis (using $`\sqrt{2}\text{Re}(Y_{lm})`$ and $`\sqrt{2}\text{Im}(Y_{lm})`$ instead of $`Y_{lm}`$ and $`Y_{l,m}`$ for $`m0`$).

The structure of the response matrix $`\alpha `$ depends on the symmetry of the molecule. For a metallic sphere of radius $`R`$ the response is isotropic, i.e., $`\alpha `$ is diagonal in the basis of the spherical harmonics: $`\alpha _{l_1m_1l_2m_2}=\delta _{l_1l_2}\delta _{m_1m_2}R^{2l_1+1}`$. Lowering the symmetry to icosahedral, $`I_h`$, introduces some anisotropy. To understand the response matrix for C60 we have to consider how the irreducible representations (IR) of the rotation group $`SO(3)`$ split into IRs of the $`I_h`$ (see table 1). An external multipole field of angular momentum $`l`$ can only give rise to a response with angular momentum $`l^{}`$, if both IRs of the $`SO(3)`$ contain a common IR of the $`I_h`$. In particular, because of parity, fields with even (odd) $`l`$ can only give rise to responses with even (odd) $`l^{}`$. Furthermore, as the irreducible representations with $`l2`$ are also irreducible with respect to the $`I_h`$, for $`l2`$ we have $`\alpha _{lmlm^{}}=\alpha _l\delta _{mm^{}}`$. Thus, restricting the multipole expansion to $`l2`$, the response of C60 is isotropic, with $`\alpha _08.1\text{bohr}`$, $`\alpha _1556\text{bohr}^3`$, and $`\alpha _244100\text{bohr}^5`$.

For $`l>2`$ the space spanned by the spherical harmonics $`Y_{lm}`$ is no longer irreducible with respect to the $`I_h`$. Thus we need to find linear combinations of the spherical harmonics that span the irreducible representations of the icosahedral group. We call them $`Y_{lxk}`$ where $`l`$ and $`x`$ denote the IR of the $`SO(3)`$ and $`I_h`$, respectively, while the index $`k`$ labels the functions within an irreducible representation of the $`I_h`$. If in the decomposition an IR $`x`$ should occur several times, we would have to introduce an additional multiplicity label. However, as can be seen from table 1, up to $`l=6`$ each IR appears at most once. We therefore suppress the multiplicity label here. Explicit expressions for the basis functions $`Y_{lxk}`$, can, e.g., be found in Ref. book:Butler, , chapter 16, or Ref. dresselhaus, , table 4.2.

In the new basis, the matrix $`\alpha `$ is built of blocks of diagonal matrices

$$\alpha _{l_1x_1k_1l_2x_2k_2}=\alpha _{l_1l_2}(x_1)\delta _{x_1x_2}\delta _{k_1k_2},$$
(8)
where $`\alpha _{l_1l_2}(x_1)`$ constitute the minimal set of parameters, and, $`\alpha `$ being real symmetric, $`\alpha _{l_1l_2}(x_1)=\alpha _{l_2l_1}(x_1)`$. The matrix elements of $`\alpha `$ were calculated up to $`l=6`$ using the results of the density functional calculations described in section II. The $`\alpha _{l_1l_2}(x_1)`$ are listed in table 2. From this minimal set of independent parameters we can determine the response for arbitrary orientations of the C60 molecule relative to the external multipole field. An example of how to do this is given in appendix A. The practical advantage of this procedure is clear: We only need to perform density functional calculations for a number of highly symmetric configurations, for which the numerical effort is much reduced. Using group theory the response for arbitrary configurations can then be determined from these special cases.

Note that the group theoretical approach presented in this section is particularly elegant in the case of a neutral molecule which has icosahedral symmetry. Upon charging, orbitals become partially filled and the symmetry is reduced which leads to a higher number of irreducible response coefficients. Furthermore, the symmetry of a charged molecule depends on how the additional charge arranges in the degenerate orbitals. This is a subtle question in the case of an isolated C60 molecule and involves Jahn-Teller effects and Coulomb interaction in competition art:Auerbach . In the present work we restrict the analysis to neutral molecules.

## IV Level splitting

In this section, a minimal set of parameters describing the Stark effect in a neutral molecule is deduced using the Wigner-Eckart theorem for the icosahedral symmetry book:Butler . This is achieved in the framework of density-functional perturbation theory art:baroni1  where an external perturbation $`V`$ induces a change in the effective potential (self-consistent field) $`V^{\text{eff}}`$ which couples the degenerate energy levels. In linear order the change is given by

$`\mathrm{\Delta }V^{\text{eff}}(𝐫)=V(𝐫)+\mathrm{\Delta }V^\text{i}(𝐫)+{\displaystyle \frac{dv_{\text{xc}}}{dn}}|_{n=n(𝐫)}\mathrm{\Delta }n(𝐫),`$ (9)
where $`\mathrm{\Delta }V^\text{i}(𝐫)=e^2𝑑𝐫^{}\mathrm{\Delta }n(𝐫^{})|𝐫𝐫^{}|^1`$ is the change in the induced potential due to the linear change $`\mathrm{\Delta }n(𝐫)`$ in the electron distribution. The last term is the exchange-correlation potential. Within linear response, the change in the charge distribution $`\mathrm{\Delta }n`$ as well as the induced potential $`\mathrm{\Delta }V^\text{i}`$ have the same symmetry (with respect to $`I_h`$) as the external perturbation $`V^\text{e}`$. Furthermore, if $`n(𝐫)`$ is the unperturbed charge distribution, then the factor $`\frac{dv_{\text{xc}}}{dn}|_{n=n(𝐫)}`$ in the last term of (9) has $`A_g`$ symmetry and does not change the symmetry of $`\mathrm{\Delta }n`$. Consequently, $`\mathrm{\Delta }V^{\text{eff}}`$ has the same symmetry as $`V^\text{e}`$ and, using the notation of the previous section, can be written as

$$\mathrm{\Delta }V^{\text{eff}}(𝐫)=\underset{lxk}{}V_{lxk}f_{xk}^l(𝐫)+𝒪[V_{lxk}^2],$$
(10)
where $`f_{xk}^l(𝐫)`$ are partner functions of the IR $`x`$ of $`I_h`$. Note that $`f_{xk}^l`$ contains spherical harmonics $`Y_{l^{}xk}`$ with all $`l^{}`$ allowed by table 1. The coupling of the degenerate levels is given by the matrix elements of $`\mathrm{\Delta }V^{\text{eff}}`$ with respect to the eigenstates of the unperturbed molecule which are denoted by $`|nxk`$, where $`n`$ is the quantum number differentiating between orbitals with the same IR $`x`$. In this context, the functions $`f_{xk}^l`$ play the role of tensor operators of $`I_h`$ and the Wigner-Eckart theorem can be used to write the matrix elements as

$`n_2x_2k_2|\mathrm{\Delta }V^{\text{eff}}|n_1x_1k_1=`$
$`{\displaystyle \underset{lxk}{}}V_{lxk}{\displaystyle \underset{\lambda }{}}t_\lambda (n_1x_1n_2x_2;lx)𝒞_{k_2k_1}^k(\lambda ;x_2x_1;x).`$
The coefficients $`𝒞_{k_2k_1}^k(\lambda ;x_2x_1;x)`$ denote the 3jm symbols (or Clebsch-Gordan coefficients) of $`I_h`$ and are entirely determined by the icosahedral symmetry. In the present work they were taken from Ref. book:Butler,  and are discussed in detail in appendix B. In order for $`𝒞_{k_2k_1}^k(\lambda ;x_2x_1;x)`$ to be non-zero, the IR $`x`$ needs to be present in the decomposition of $`x_1x_2`$. If $`x`$ occurs more than once in $`x_1x_2`$ then the multiplicity index $`\lambda `$ is required. From this selection rule follows again that even potentials couple linearly whereas odd potentials couple only in second order (quadratic Stark effect). Finally, the factors $`t_\lambda (n_1x_1n_2x_2;lx)`$ are the coupling constants which constitute the minimal set of parameters describing the level splitting.

We analyze in more detail the level splitting of the LUMO ($`x=T_{1u}`$) and HOMO ($`x=H_u`$) for $`l=1`$ and $`l=2`$ external potentials which correspond to an electric field and a quadrupole potential. We will see in Section V that these two multipoles are dominant in the case of a charged C60 layer exposed to an electric field. Within the icosahedral symmetry $`I_h`$, the two potentials form partners for the IR $`T_{1u}`$ and $`H_g`$ respectively (see table 1). Details of the calculation are given in appendix B. The coupling constant derived from the level splitting calculated by DFT are given in table 3. In the case of the odd $`l=1`$ potential, only second-order coupling to closest-by orbitals was considered ($`T_{1g}`$ for the LUMO, $`H_g`$ and $`G_g`$ for the HOMO, see figure 3), which, however, as shown below, gives very satisfactory results. For the splitting of the HOMO under the $`l=2`$ ($`x=H_g`$) potential, there are two coupling constants because $`H_g`$ occurs twice in the product $`H_uH_u=A_gT_{1g}T_{2g}2G_g2H_g`$. In appendix B the coupling matrix $`H_{T_{1u}}^{}`$ ($`H_{H_u}^{}`$) describing the level splitting of the LUMO (HOMO), due to an applied $`(lm)=(10)`$ and $`(lm)=(20)`$ potential along an arbitrary direction of the molecule, is calculated using perturbation theory. To a very good approximation, the result can be cast into the form

$`H_{T_{1u}}^{}(\theta )`$ $`=`$ $`\left(V_{10}^2c_1+V_{20}c_2\right)𝒞_{T_{1u}}(\theta ),`$ (12)
$`H_{H_u}^{}(\theta )`$ $`=`$ $`\left(V_{10}^2d_1+V_{20}d_2\right)𝒞_{H_u}(\theta ),`$ (13)
where $`c_1`$, $`c_2`$, $`d_1`$ and $`d_2`$ are constants depending on the coupling constants of table 3 and the energies given in figure 3. $`𝒞_{T_{1u}}(\theta )`$ and $`𝒞_{H_u}(\theta )`$ are matrices within the LUMO and HOMO subspace respectively and depend on the angle $`\theta `$ between the $`z`$-direction of the $`(lm)=(10)`$$`(20)`$ potentials and the 5-fold axis of the molecule. (for more details see appendix B). In figure 4 and figure 5 the splitting of the LUMO and HOMO is shown using the previous relations along with the points calculated by DFT. The group-theoretical fit is very satisfactory. Note that the splitting of the LUMO is independent of the orientation of the molecule when only coupling among the LUMO or to the $`T_{1g}`$ is considered (see appendix B). Furthermore, relations (12) and (13) are quite remarkable, as they imply that the contributions of the $`l=1`$ and $`l=2`$ potential lift the degeneracy of the molecular levels in the same way, and thus the total splitting is given by the sum of the $`l=1`$ and $`l=2`$ splittings.

## V Molecules in a Layer

In this section a simple model of a C60 FET is considered. As mentioned above a FET can be understood as capacitance where one plate is the gate and the other plate is the material (here C60) which is investigated. An analysis of this device, in particular the calculation of the charge distribution, was done in a previous work art:wehrli1 . It was found, that the charge concentrates on the first layer in the high doping regime. In what follows, we take this as a motivation to consider a single layer of C60 molecules which acts as a plate of a capacitance. We use the response of a *neutral* molecule (table 2) to describe the electrostatic behavior of the molecules in the layer. The doped charge is taken care of by adding a monopole on every site. The molecules are then exposed to the electric field arising from the gate as well as the monopole fields of the neighboring molecules. In order to simplify the calculation, we assume a perfect lattice, either square, for the (001), or triangular, for the (111) plane of the fcc lattice formed by the C60 molecules in the bulk. In order for the sites to be equivalent we neglect the non-spherically symmetric part of the response of the C60-molecule and use the averaged response given by

$$\alpha _{l_1m_1l_2m_2}^{\text{spher.-sym.}}=\delta _{l_1l_2}\delta _{m_1m_2}\underset{xk}{}\frac{\alpha _{l_1xkl_1xk}}{2l_1+1},$$
(14)
where the sum is taken over all components of the $`l_1`$ subspace. Since the response of C60 for multipoles $`l2`$ is isotropic (cf. section III), this averaging is exact for $`l=1,2`$.

Because of translational invariance, the total potential is given by

$$V_{\text{tot}}(𝐫)=V^\text{e}(𝐫)+\underset{𝐑_i}{}V^\text{i}(𝐫𝐑_i),$$
(15)
where the sum is taken over all lattice sites $`𝐑_i`$. At a given site, say $`𝐑_i=0`$, the total potential can also be decomposed as $`V_{\text{tot}}(𝐫)=V^{\text{scr}}(𝐫)+V^\text{i}(𝐫)`$ where the screened potential $`V^{\text{scr}}(𝐫)=V^\text{e}(𝐫)+_{𝐑_i0}V^\text{i}(𝐫𝐑_i)`$ contains the external potential $`V^\text{e}(𝐫)`$ as well as the sum of all induced potentials of the other molecules. This sum depends linearly on the induced potential $`V^\text{i}`$ and, hence, the coefficients of $`V^{\text{scr}}`$ are given by

$$V_{l_1m_1}^{\text{scr}}=V_{l_1m_1}+\underset{l_2m_2}{}\beta _{l_1m_1l_2m_2}Q_{l_2m_2}.$$
(16)
The matrix $`\beta `$ is entirely given by geometry and discussed in appendix C. On the other hand, the screened potential $`V^{\text{scr}}`$ induces a potential $`\mathrm{\Delta }V^\text{i}`$ as given in equation (5). Therefore the total induced potential $`V^\text{i}=V^{\text{i,0}}+\mathrm{\Delta }V^\text{i}`$ is

$$Q_{l_1m_1}=Q_{l_1m_1}^0\underset{l_2m_2}{}\alpha _{l_1m_1l_2m_2}V_{l_2m_2}^{\text{scr}}.$$
(17)
Equations (16) and (17) can be combined by eliminating the coefficients $`Q_{lm}`$ which yields

$`V_{l_1m_1}^{\text{scr}}`$ $`=`$ $`{\displaystyle \underset{l_2m_2}{}}[1+\beta \alpha ]_{l_1m_1l_2m_2}^1V_{l_2m_2}^{\text{bare}},`$ (18)
$`V_{l_1m_1}^{\text{bare}}`$ $`=`$ $`V_{l_1m_1}+{\displaystyle \underset{l_2m_2}{}}\beta _{l_1m_1l_2m_2}Q_{l_2m_2}^0,`$ (19)
where the $`V_{l_1m_1}^{\text{bare}}`$ describes the bare potential arising from the external potential (the electric field of the gate) and intrinsic moments of the molecules (the induced charge, i.e. monopoles). The square (triangular) lattice in the presence of the electric field has the rotational symmetry $`C_{4v}`$ ($`C_{6v}`$). As a consequence, only components with these symmetries are non-zero and therefore they are given by $`\text{Re}(Y_{lm})`$ with m a multiple of 4 (6). Using relation (18) and (19) the screened potential can be calculated. The non-zero components entering (19) are the monopole charge $`Q_{00}^0=q`$ and electric field $`V_{10}=E_{\text{Gate}}`$. As the FET is overall neutral $`E_{\text{Gate}}=2\pi q/A_{\text{mol}}`$ with $`A_{(001)}=a^2/2`$ for the square lattice and $`A_{(111)}=\sqrt{3}a^2/4`$ for the triangular lattice. The results are given in table 4 and graphically depicted in figure 6. The components $`(lm)=(10),(20)`$ are most dominant and higher ones are at least one order of magnitude smaller. This justifies a posteriori the assumption of spherical symmetry because (14) is exact for $`(lm)=(10),(20)`$. From figure 6 it can be seen that the electric field is efficiently screened within the layer. Note that decrease of the electric field yields negative sign of $`v_{20}`$ in table 4. We also have checked the influence of adjacent layers of C60 and of a close-by dielectric (with a dielectric constant $`\epsilon =10`$). The effects on the parameters in table 4 were less than 2 %. The reason is that the field inhomogeneities induced by a 2D lattice of multipoles decay exponentially outside the lattice.

## VI Splitting in self-consistent multipole field

We are now in the position to estimate the effect of the external field on the electronic structure of the C60 molecules in the monolayer that carries charge. To do so we have performed density functional calculations for a molecule in the self-consistent multipole fields as determined in the previous section (cf. table 4). Figure 7 shows the splitting of the molecular levels in the self-consistent field for a (001) monolayer, where the molecule is oriented with one of its two-fold axes pointing in the direction of the external field. The maximum energy difference between split states is given in table 5 and compared to the result from the perturbative formula (12) and (13), which is in good agreement for $`|q|2`$. As figure 8 demonstrates, similar results are obtained for other geometries. Qualitatively, the results are also in agreement with an approximative tight-binding calculation published earlier in Ref. art:swehrli02, . As expected we find that the stronger the external field, i.e., the larger the induced charge, the stronger the splitting. We notice also a pronounced asymmetry in the splitting: when the monolayer is charged with electrons the splitting is different from when it is charged with holes. Again, the reason is parity: Because of parity an external homogeneous field, and more generally any multipole potential with $`l`$ odd, leads to a quadratic Stark effect. Hence, for odd $`l`$ the splitting is independent of the sign of the field. For $`l`$ even, however, the levels split already in first order (linear Stark effect), so the splitting changes sign with the external field. Thus the asymmetry originates from the multipoles with even $`l`$. Moreover, because of the first-order versus second-order effect, even though the largest even multipole (20) is significantly smaller than the largest odd multipole (the screened external field), it contributes considerably to the splitting.

We have seen in section IV that for the HOMO and LUMO of C60 the splittings caused by (10) and (20) multipole potentials are essentially additive (cf. equation (IV)). Thus we expect them to add up for one sign of the external field, while they should compensate for the opposite sign. In fact, for electron doping the splittings of the HOMO seem to almost perfectly cancel, while upon hole doping the splitting of the HOMO is essentially doubled compared to the splitting caused by the screened homogeneous field alone. Thus for the HOMO, the splittings happen to add up for an external field that induces charge carriers into that orbital — a situation that is particularly unfavorable for hole doping. For the LUMO the situation is similar: For an external field that induces electrons in the LUMO the splitting is enhanced. So it turns out that the contributions of the higher multipoles conspire to enhance the splitting of the orbital that carries the induced charge. For both, HOMO and LUMO, the splitting becomes comparable to the respective band width for a field that corresponds to about two charge carriers per molecule.

The calculations reported above have been performed for an uncharged molecule. Considering instead the splitting for a molecule that carries the proper induced charge, the splittings are substantially reduced. This is easily understood: Due to the splitting, the electrons will fill only the lowest of the split levels. But since the interaction between electrons in the same orbital is larger than the interaction of electrons in different orbitals, the occupied levels will be shifted upwards, compared to the levels that were left empty – thereby reducing the splitting. We are, however, not interested in the splitting per se, but in the effect of the splitting on the band structure of a monolayer. Thus allowing electrons only in the energetically lowest levels would mean that in the lattice, the electrons are not allowed to hop to energetically higher ones of the split levels. This implies that the original band structure would already be separated into a set of bands originating from occupied, and another set of bands originating from the empty orbitals. To eliminate the undesired differences in the interaction between electrons in the split orbitals, we therefore work with the splittings obtained for a neutral molecule.

To estimate the effect of the Stark splitting on the density of states (DOS), we have performed tight-binding calculations for the (001) monolayer (square lattice), assuming the unidirectional structure (two-fold axis of the molecules pointing in the direction of the external field). The basis for the tight-binding Hamiltonian and the hopping matrix elements were taken from the parametrization given in Ref. art:satpathy92, . The splittings shown in figure 7 were then used to derive an on-site coupling between the different orbitals. In the case of the LUMO, the on-site coupling is diagonal and reduces to orbital dependent on-site energies. The DOS for the LUMO ($`T_{1u}`$) and the HOMO ($`H_u`$) bands calculated with this model for different charging are shown in figure 9 and 10, respectively. We find that already for an induced charge of one carrier per molecule the change with respect to the unperturbed DOS is sizable. For $`q=2`$ one of the $`T_{1u}`$-bands is already completely separated from the other two. Also the HOMO density of states shows for $`q=2`$ hardly any resemblance to the original DOS, and for $`q=3`$ also the $`H_u`$-bands fall into two groups. We thus conclude that beyond filling $`|q|=2`$ the electronic structure is distorted so much compared to the unperturbed monolayer that one can no longer speak of doping.

The calculation of the density of states in a minimal tight-binding basis involves, of course, approximations: First, in the lattice, not only hopping between LUMO (HOMO) levels is allowed, but also hopping via energetically close-by levels. An orbital at $`\mathrm{\Delta }\epsilon `$ away will give a contribution to the hopping of about $`t^2/\mathrm{\Delta }\epsilon `$, where $`t`$ is the hopping matrix element from the orbital, that we consider explicitly, to the orbital at $`\mathrm{\Delta }\epsilon `$. The influence of this effect on the hopping between molecules was studied in Ref. art:satpathy92,  and changes of the order of 5 % were found. More importantly, due to the deformation of the molecular orbitals in the field, the hopping matrix elements between the HOMO or LUMO orbitals will change. For a simple estimate, we have performed tight-binding calculations of a C60 molecule in an external homogeneous field and determined the average hopping matrix element between the $`t_{1u}`$-orbitals following the approach of Ref. hopp, . We find that the change in the hopping matrix elements depends strongly on the orientation of the molecules. Typical changes are of the order of 10–20%.

## VII Conclusions

We have analyzed the changes in the electronic structure of a C60 monolayer in which charge carriers are induced by the application of an external homogeneous electric field. We find that the effective field seen by each molecule in the monolayer is strongly screened, but that there are additional higher multipole potentials. Although these components are considerably weaker than the screened homogeneous field, for even $`l`$, they give a significant contribution to the level splitting as they are of first order. In addition the $`l=1`$ (homogeneous field) and $`l=2`$ potentials split the HOMO and LUMO in almost the same way, so the splittings they produce add up or counteract, depending on the sign of the external field. For both, the HOMO and the LUMO, the signs turn out to be such, that the splitting is enhanced when the charge is induced in the respective level. Thus the level that carries the field-induced charge is strongly changed by the effective field — a particularly unfavorable situation if one wants to achieve doping, i.e., filling of a level without substantially changing its electronic structure.

There are, of course, some effects that have been neglected in our analysis: The polarizability of C60 increases when charge is put in the $`T_{1u}`$ orbital. Thus when filling the LUMO the screening of the external field should become somewhat more efficient. Second, due to the electric field the molecular orbitals are deformed and thus the hopping matrix elements between neighboring molecules change. This effect depends sensitively on the orientation of the molecules but it typically leads to a slight decrease of the band width with applied field, making the Stark splitting even more important. Furthermore we have neglected the effect of electron-phonon coupling, which also tends to narrow the bands. This should be particularly important for the $`H_u`$-band.elph  Finally, we have not considered correlation effects. It is known that, at integer filling, the alkali-doped Fullerenes are close to a Mott transition and that they are metallic because of orbital degeneracy.c60mott  Due to the reduced coordination in two dimensions, a doped monolayer of C60 should be even more strongly correlated, even if the orbitals are still degenerate. Lifting the degeneracy, e.g., by the Stark effect significantly increases correlations even further.manini  It therefore seems inevitable that for an integer number of induced charges a C60 monolayer will be a Mott insulator, even though it is not clear what effect the gate oxide might have on the Coulomb repulsion $`U`$ between two electrons on a C60 molecule.kirchb03

To summarize, it seems safe to conclude that the physics of field-effect devices based on C60 as active material should be quite different from that of the alkali-doped fullerenes, at least at doping levels beyond two carriers per molecule.

###### Acknowledgements.

The authors would like to thank O. Gunnarsson, T.M. Rice, B. Batlogg, and C. Helm for fruitful discussions. This work has been supported by the Swiss Nationalfonds.

## Appendix A Response for general orientation

To illustrate the use of the response coefficients given in table 2, we show, how to calculate the response for a multipole potential with $`(lm)=(30)`$ for different orientations of the $`z`$-axis with respect to the C60 molecule. This is the first non-trivial case, as the response for multipoles with $`l2`$ is isotropic. To start with, we need the basis functions $`Y_{lxk}`$ spanning the irreducible representation of the $`I_h`$, that were introduced in section III. For specific orientations, these can be found, e.g., in the Ref. book:Butler, , chapter 16, or Ref. dresselhaus, , table 4.2. For arbitrary orientations, they have to be derived by explicitly finding the basis functions that span the irreducible representations of the icosahedral group. For the sake of the example, we consider the response for $`z`$ along the 5-fold and the 3-fold axis of the C60 molecule. The corresponding basis functions are then both found in Ref. book:Butler, . Since we are interested in the response to an external multipole with $`(lm)=(30)`$, we have to identify those basis functions that contain $`Y_{30}`$. They are shown in table 6. It turns out that for $`z`$ parallel to the 5-fold axis the $`(lm)=(30)`$ potential corresponds to a pure $`T_{2u}`$ potential. Thus the response is given by $`\mathrm{\Delta }Q_{30}=\alpha _{33}(T_{2u})V_{30}`$. For $`z`$ along the 3-fold axis the situation is more complicated, as the potential is now a mixture of partner functions of $`G_u`$ and $`T_{2u}`$: As can be seen from table 6, it mixes with the $`Y_{3c3}=\sqrt{2}\text{Re}(Y_{33})`$ component. By construction, the response matrix $`\alpha `$ in the subspace spanned by $`\{Y_{30},Y_{3c3}\}`$ is diagonal, with diagonal elements $`\alpha _{33}(H_u)`$ and $`\alpha _{33}(T_{2u})`$, when written in the basis functions listed in table 6. However, to obtain the multipole response we have to use the basis $`(Y_{30},Y_{3c3})`$, for which there are off-diagonal elements:

$$\alpha =\frac{1}{9}\left(\begin{array}{cc}8\alpha _{33}(H_u)+\alpha _{33}(T_{2u})& \sqrt{8}[\alpha _{33}(T_{2u})\alpha _{33}(H_u)]\\ \sqrt{8}[\alpha _{33}(T_{2u})\alpha _{33}(H_u)]& \alpha _{33}(H_u)+8\alpha _{33}(T_{2u})\end{array}\right).$$
(20)
We thus find

$`\mathrm{\Delta }Q_{30}`$ $`=`$ $`{\displaystyle \frac{1}{9}}\left[8\alpha _{33}(H_u)+\alpha _{33}(T_{2u})\right]V_{30},`$
$`\mathrm{\Delta }Q_{3c3}`$ $`=`$ $`{\displaystyle \frac{\sqrt{8}}{9}}[\alpha _{33}(T_{2u})\alpha _{33}(H_u)]V_{30}.`$
## Appendix B Coupling Matrices

In this section we discuss the calculation of the level splitting for arbitrary directions within perturbation theory using the coupling constants of table 3. As discussed above, we restrict the analysis to $`l=1`$ and $`l=2`$ external potentials, which corresponds to $`x=T_{1u}`$ and $`x=H_g`$ potentials in the icosahedral symmetry $`I_h`$. In first order perturbation theory, the splitting of the levels in the degenerate subspace $`_{nx}`$ is given by the eigenvalues of the matrix

$$H_{k_1k_2}^{(1)}(nx)=nxk_2|V^{\text{eff}}|nxk_1.$$
(21)
This matrix vanishes in the case of an odd potential and the splitting is given by the second order expression

$$H_{k_1k_2}^{(2)}(nx)=\underset{\genfrac{}{}{0pt}{}{(n^{}x^{})}{(nx)}}{}\underset{k^{}}{}\frac{nxk_2|V^{\text{eff}}|n^{}x^{}k^{}n^{}x^{}k^{}|V^{\text{eff}}|nxk_1}{E_{nx}E_{n^{}x^{}}}.$$
(22)
The matrix-elements in (21) and (22) are given in (IV) and involve the icosahedral Clebsch-Gordan coefficients $`𝒞_{k_2k_1}^k(\lambda ;x_2x_1;x)`$. In order for the $`k`$-indices to be defined we consider the molecule oriented with the 5-fold axis parallel to the $`z`$-axis (see figure 1). This allows to label the states within a multiplet unambiguously with its $`C_5`$ index $`k`$. The ordered basis of a $`T_{1u}`$ subspace has $`k`$-indices $`(0,1,1)`$ whereas the ordered basis of an $`H_u`$ subspace is given by $`(0,1,1,2,2)`$. Note that in the case of applied $`l=1`$ or $`l=2`$ potential, the $`k`$ index corresponds to the $`m`$ index of the spherical harmonics. For a detailed discussion we refer to Ref. book:Butler, . We will present the coefficients $`𝒞_{k_2k_1}^k(\lambda ;x_2x_1;x)`$ as matrices with respect to the indices $`k_1`$ and $`k_2`$. In order to reduce the number of matrices, we will give the coupling matrices for $`(lm)=(10),(20)`$ potentials, which are rotated around the $`y`$-axis by an angle $`\theta `$. The resulting matrices are then given by

$$𝒞^\theta (\lambda ;x_2x_1;x)=\underset{k}{}R_{0k}^l(\theta )𝒞^k(\lambda ;x_2x_1;x),$$
(23)
where $`R_{k^{}k}^l(\theta )`$ is the rotation matrix of the spherical harmonics in a given $`l`$-subspace. Using the previous relations and equation (IV), the coupling matrix for an even $`l=2`$ potential is given by

$$H^{(1)}(nx,\theta )=V_{20}\underset{\lambda }{}t_\lambda (nxnx;H_g)𝒞_\lambda ^\theta (xx;H_g).$$
(24)
Note that the multiplicity label $`\lambda `$ is only relevant for the HOMO $`H_u`$ because $`H_g`$ occurs twice in the product $`H_uH_u=A_gT_{1g}T_{2g}2G_g2H_g`$. The coupling matrix for the odd $`l=1`$ potential is

$$H^{(2)}(nx,\theta )=V_{10}^2\underset{\genfrac{}{}{0pt}{}{(n^{}x^{})}{(nx)}}{}\frac{t(nxnx^{};T_{1u})^2}{E_{nx}E_{n^{}x^{}}}𝒞^\theta (xx^{};T_{1u})^\text{T}𝒞^\theta (xx^{};T_{1u}).$$
(25)
In the following we restrict the sum over subspaces $`_{n^{}x^{}}`$ closest in energy to $`_{nx}`$ which is the $`T_{1g}`$ subspace in the case of the LUMO and the $`H_g`$ and $`G_g`$ subspaces in the case of the the HOMO (see Fig 3). Below, the coupling matrices $`𝒞_\lambda ^\theta (x_1x_2;x)`$ which are needed to calculate the splitting of the HOMO and LUMO are given. They are traceless $`\text{Tr}𝒞^\theta =0`$ and normalized such that $`\text{Tr}(C^\theta )^\text{T}𝒞^\theta =1`$. The coupling matrices which describe the splitting of the LUMO are

$`𝒞^\theta (T_{1u}T_{1g};T_{1u})`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\begin{array}{ccc}0& \mathrm{sin}\theta & \mathrm{sin}\theta \\ \mathrm{sin}\theta & \sqrt{2}\mathrm{cos}\theta & 0\\ \mathrm{sin}\theta & 0& \sqrt{2}\mathrm{cos}\theta \end{array}\right),`$ (29)
$`𝒞^\theta (T_{1u}T_{1u};H_g)`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(\begin{array}{ccc}4\sqrt{\frac{2}{3}}+2\sqrt{6}\mathrm{sin}^2\theta & \sqrt{3}\mathrm{sin}2\theta & \sqrt{3}\mathrm{sin}2\theta \\ \sqrt{3}\mathrm{sin}2\theta & 2\sqrt{\frac{2}{3}}\sqrt{6}\mathrm{sin}^2\theta & \sqrt{6}\mathrm{sin}^2\theta \\ \sqrt{3}\mathrm{sin}2\theta & \sqrt{6}\mathrm{sin}^2\theta & 2\sqrt{\frac{2}{3}}\sqrt{6}\mathrm{sin}^2\theta \end{array}\right).`$ (33)
The eigenvalues of these two matrices are independent of $`\theta `$ and given by $`(0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})`$ and $`(\frac{2}{\sqrt{6}},\frac{1}{\sqrt{6}},\frac{1}{\sqrt{6}})`$ which implies that the splitting is independent of the orientation of the molecule with respect to the direction of the applied $`l=1`$ and $`l=2`$ potentials. The coupling of the HOMO ($`H_u`$) to the lower lying $`H_g`$ and $`G_g`$ levels is given by

$`𝒞^\theta (H_uH_g;T_{1u})`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{5}}}\left(\begin{array}{ccccc}0& \sqrt{3}\mathrm{sin}\theta & \sqrt{3}\mathrm{sin}\theta & 0& 0\\ \sqrt{3}\mathrm{sin}\theta & \sqrt{2}\mathrm{cos}\theta & 0& \sqrt{2}\mathrm{sin}\theta & 0\\ \sqrt{3}\mathrm{sin}\theta & 0& \sqrt{2}\mathrm{cos}\theta & 0& \sqrt{2}\mathrm{sin}\theta \\ 0& \sqrt{2}\mathrm{sin}\theta & 0& 2\sqrt{2}\mathrm{cos}\theta & 0\\ 0& 0& \sqrt{2}\mathrm{sin}\theta & 0& 2\sqrt{2}\mathrm{cos}\theta \end{array}\right),`$ (39)
$`𝒞^\theta (H_uG_g;T_{1u})`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{5}}}\left(\begin{array}{ccccc}\sqrt{3}\mathrm{sin}\theta & 2\sqrt{2}\mathrm{cos}\theta & 0& \frac{1}{\sqrt{2}}\mathrm{sin}\theta & 0\\ \sqrt{3}\mathrm{sin}\theta & 0& 2\sqrt{2}\mathrm{cos}\theta & 0& \frac{1}{\sqrt{2}}\mathrm{sin}\theta \\ 0& \sqrt{2}\mathrm{sin}\theta & 0& \sqrt{2}\mathrm{cos}\theta & \frac{3}{\sqrt{2}}\mathrm{sin}\theta \\ 0& 0& \sqrt{2}\mathrm{sin}\theta & \frac{3}{\sqrt{2}}\mathrm{sin}\theta & \sqrt{2}\mathrm{cos}\theta \end{array}\right).`$ (44)
The coupling of the HOMO among themselves is given by

$`𝒞^\theta (1;H_uH_u;H_g)`$ $`=`$ $`{\displaystyle \frac{1}{4\sqrt{5}}}\left(\begin{array}{ccccc}8+12\mathrm{sin}^2\theta & \sqrt{\frac{3}{2}}\mathrm{sin}2\theta & \sqrt{\frac{3}{2}}\mathrm{sin}2\theta & \sqrt{\frac{3}{2}}\mathrm{sin}^2\theta & \sqrt{\frac{3}{2}}\mathrm{sin}^2\theta \\ \sqrt{\frac{3}{2}}\mathrm{sin}2\theta & 23\mathrm{sin}^2\theta & 3\mathrm{sin}^2\theta & 3\mathrm{sin}2\theta & 3\mathrm{sin}^2\theta \\ \sqrt{\frac{3}{2}}\mathrm{sin}2\theta & 3\mathrm{sin}^2\theta & 23\mathrm{sin}^2\theta & 3\mathrm{sin}^2\theta & 3\mathrm{sin}2\theta \\ \sqrt{\frac{3}{2}}\mathrm{sin}^2\theta & 3\mathrm{sin}2\theta & 3\mathrm{sin}^2\theta & 23\mathrm{sin}^2\theta & 3\mathrm{sin}2\theta \\ \sqrt{\frac{3}{2}}\mathrm{sin}^2\theta & 3\mathrm{sin}^2\theta & 3\mathrm{sin}2\theta & 3\mathrm{sin}2\theta & 23\mathrm{sin}^2\theta \end{array}\right),`$ (50)
$`𝒞^\theta (2;H_uH_u;H_g)`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(\begin{array}{ccccc}0& \sqrt{\frac{3}{2}}\mathrm{sin}2\theta & \sqrt{\frac{3}{2}}\mathrm{sin}2\theta & \sqrt{\frac{3}{2}}\mathrm{sin}^2\theta & \sqrt{\frac{3}{2}}\mathrm{sin}^2\theta \\ \sqrt{\frac{3}{2}}\mathrm{sin}2\theta & 23\mathrm{sin}^2\theta & \mathrm{sin}^2\theta & \mathrm{sin}2\theta & \mathrm{sin}^2\theta \\ \sqrt{\frac{3}{2}}\mathrm{sin}2\theta & \mathrm{sin}^2\theta & 23\mathrm{sin}^2\theta & \mathrm{sin}^2\theta & \mathrm{sin}2\theta \\ \sqrt{\frac{3}{2}}\mathrm{sin}^2\theta & \mathrm{sin}2\theta & \mathrm{sin}^2\theta & 2+3\mathrm{sin}^2\theta & \mathrm{sin}2\theta \\ \sqrt{\frac{3}{2}}\mathrm{sin}^2\theta & \mathrm{sin}^2\theta & \mathrm{sin}2\theta & \mathrm{sin}2\theta & 2+3\mathrm{sin}^2\theta \end{array}\right).`$ (56)
In equation (25) the product $`𝒞^\theta (xx^{};T_{1u})^\text{T}𝒞^\theta (xx^{};T_{1u})`$ enters. For the matrices given in (29), (39) and (44) these products can be expressed in terms of the $`l=2`$ coupling matrices (33), (50) and (56):

$`𝒞^\theta (T_{1u}T_{1g};T_{1u})^\text{T}𝒞^\theta (T_{1u}T_{1g};T_{1u})`$ $`=`$ $`{\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{\sqrt{6}}}𝒞^\theta (T_{1u}T_{1u};H_g),`$ (57)
$`𝒞^\theta (H_uH_g;T_{1u})^\text{T}𝒞^\theta (H_uH_g;T_{1u})`$ $`=`$ $`{\displaystyle \frac{1}{5}}+{\displaystyle \frac{\sqrt{5}}{10}}𝒞^\theta (1;H_uH_u;H_g){\displaystyle \frac{3}{10}}𝒞^\theta (2;H_uH_u;H_g),`$ (58)
$`𝒞^\theta (H_uG_g;T_{1u})^\text{T}𝒞^\theta (H_uG_g;T_{1u})`$ $`=`$ $`{\displaystyle \frac{1}{5}}+{\displaystyle \frac{\sqrt{5}}{10}}𝒞^\theta (1;H_uH_u;H_g)+{\displaystyle \frac{3}{10}}𝒞^\theta (2;H_uH_u;H_g).`$ (59)
Using these relations, the total coupling matrix (neglecting the constant terms in the previous relations) due to the applied $`l=1`$ and $`l=2`$ potential, is given by

$`H^{}(T_{1u},\theta )`$ $`=`$ $`\left(V_{10}^2c_1+V_{20}c_2\right)𝒞^\theta (T_{1u}T_{1u};H_g),`$ (60)
$`H^{}(H_u,\theta )`$ $`=`$ $`V_{10}^2d_1\left[\mathrm{cos}\delta _1𝒞^\theta (1;H_uH_u;H_g)+\mathrm{sin}\delta _1𝒞^\theta (2;H_uH_u;H_g)\right]+`$ (62)
$`V_{20}d_2\left[\mathrm{cos}\delta _2𝒞^\theta (1;H_uH_u;H_g)+\mathrm{sin}\delta _2𝒞^\theta (2;H_uH_u;H_g)\right],`$
where $`c_1=\frac{1}{\sqrt{6}}\frac{t(T_{1u}T_{1g};T_{1u})^2}{E_{T_{1u}}E_{T_{1g}}}`$ and $`c_2=t(T_{1u}T_{1u};H_g)`$. Similarly we have $`d_1\mathrm{cos}\delta _1=\frac{\sqrt{5}}{10}\left[\frac{t(H_uH_g;T_{1u})^2}{E_{H_u}E_{H_g}}+\frac{t(H_uG_g;T_{1u})^2}{E_{H_u}E_{G_g}}\right]`$, $`d_1\mathrm{sin}\delta _1=\frac{3}{10}\left[\frac{t(H_uH_g;T_{1u})^2}{E_{H_u}E_{H_g}}+\frac{t(H_uG_g;T_{1u})^2}{E_{H_u}E_{G_g}}\right]`$ and $`d_2\mathrm{cos}\delta _2=t_1(H_uH_u;H_g)`$, $`d_2\mathrm{sin}\delta _2=t_2(H_uH_u;H_g)`$. Equation (60) implies that the contributions of the $`l=1`$ and $`l=2`$ to the splitting of the LUMO add up trivially. Using the values in table 3 and the energies of figure 3 yields the values $`\delta _1=0.064`$ and $`\delta _2=0.037`$, which are almost equal when compared to $`\pi `$. This can be understood by the following remarks: $`\delta _1=0`$ in the case of $`t(H_uH_g;T_{1u})=t(H_uG_g;T_{1u})`$ and $`E_{H_g}=E_{G_g}`$. Furthermore, it can be shown that $`\delta _2=\mathrm{arctan}(1/\sqrt{125})0.090`$ assuming that the angular dependence of the HOMO is given by $`l=5`$ spherical harmonics. Taking an average value of $`\delta =0.050`$ yields the approximate relation

$`H^{}(H_u,\theta )`$ $``$ $`\left(V_{10}^2d_1+V_{20}d_2\right)\left[\mathrm{cos}\delta 𝒞^\theta (1;H_uH_u;H_g)+\mathrm{sin}\delta 𝒞^\theta (2;H_uH_u;H_g)\right],`$ (63)
which shows that, to a good approximation, the contributions of the $`l=1`$ and $`l=2`$ potential to the splitting of the HOMO add up trivially.

## Appendix C Calculation of the matrix $`\beta `$

In this section it is shown how to calculate the matrix $`\beta `$ appearing in equation (16). The second term on the right side of this equation describes the coefficients of the term $`_{𝐑_i0}V^\text{i}(𝐫𝐑_i)`$ which enters the screened potential $`V^{\text{scr}}`$ and which describes the potential induced by all neighboring sites. Using the definition (1), we can rewrite this expression as

$$\underset{𝐑_i0}{}V^\text{i}(𝐫𝐑_i)=\underset{lm}{}Q_{lm}\underset{𝐑_i0}{}I_{lm}^{}(𝐫𝐑_i)=\underset{lm}{}Q_{lm}\underset{𝐑_i0}{}(1)^mI_{lm}(𝐫𝐑_i).$$
(64)
The function $`I_{lm}(𝐫_1𝐫_2)`$ can be decomposed using the translation formula (for $`r_1<r_2`$)

$$I_{LM}(𝐫_1𝐫_2)=\underset{\genfrac{}{}{0pt}{}{l_1,l_2=0}{l_2l_1=L}}{\overset{\mathrm{}}{}}(1)^{l_2}\sqrt{\frac{(2l_2+1)!}{(2L+1)!(2l_1)!}}\underset{m_1,m_2}{}C_{l_1m_1l_2m_2}^{LM}R_{l_1m_1}(𝐫_1)I_{l_2m_2}(𝐫_2),$$
(65)
where $`C_{l_1m_1l_2m_2}^{LM}`$ denote the Clebsch-Gordan coefficients. This formula can be found in different forms in the literature, see for example Ref. book:QToAngularMomentum, ; art:epton94, . Substituting (65) in the sum (64) with $`𝐫_1=𝐫`$ and $`𝐫_2=𝐑_i`$ yields

$$\underset{𝐑_i0}{}V^\text{i}(𝐫𝐑_i)=\underset{\genfrac{}{}{0pt}{}{l_1m_2}{l_2m_2}}{}\beta _{l_1m_1l_2m_2}Q_{l_2m_2}R_{l_1m_1}(𝐫),$$
(66)
where the matrix $`\beta `$ is given by

$`\beta _{l_1m_1l_2m_2}`$ $`=`$ $`(1)^{m_2+l_2}\sqrt{{\displaystyle \frac{[2(l_1+l_2)]!}{(2l_1)!(2l_2)!}}}C_{l_1m_1l_2m_2}^{l_1+l_2m_1m_2}{\displaystyle \underset{𝐑_i0}{}}I_{l_1+l_2m_1m_2}(𝐑_i)`$
$`=`$ $`(1)^{m_2+l_2}\sqrt{{\displaystyle \frac{(l_1+l_2+m_1m_2)!(l_1+l_2m_1+m_2)!}{(l_1+m_1)!(l_1m_1)!(l_2+m_2)!(l_2m_2)!}}}{\displaystyle \underset{𝐑_i0}{}}I_{l_1+l_2m_1m_2}(𝐑_i)`$
In the last equality in (C) the explicit form of the Clebsch-Gordan coefficients was used book:QToAngularMomentum . One verifies, that the matrix $`\beta `$ is complex conjugate under the exchange of all indices. The remaining sums over the lattice sites $`𝐑_i`$ in (C) are easily performed by computer.