text:
# Useful Alternative to the Multipole Expansion of 1/𝑟 Potentials

## I Introduction

For pairwise Coulomb or gravitational potentials, one often expands the inverse distance between two points $`𝐱`$ and $`𝐱^{}`$ in the standard multipole form ref1

$$\frac{1}{|𝐱𝐱^{}|}=\frac{1}{\sqrt{rr^{}}}\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\left(\frac{r_<}{r_>}\right)^{\mathrm{}+\frac{1}{2}}P_{\mathrm{}}(\mathrm{cos}\gamma ),$$
(1)
where $`r_<`$ $`(r_>)`$ is the smaller (larger) of the spherical distances $`r`$ and $`r^{}`$, and $`P_{\mathrm{}}(\mathrm{cos}\gamma )`$ is the Legendre polynomial ref2  with argument

$$\mathrm{cos}\gamma \widehat{𝐱}\widehat{𝐱}^{}=\mathrm{cos}\theta \mathrm{cos}\theta ^{}+\mathrm{sin}\theta \mathrm{sin}\theta ^{}\mathrm{cos}(\varphi \varphi ^{}).$$
(2)
The set of six coordinates {$`𝐱,𝐱^{}`$} may be viewed either as defining two points relative to the origin or as the coordinates of a three-body system once the motion of the center of mass has been separated. In the “body frame”, three out of the six coordinates are dynamical variables, the potential energy depending only on them. Of various choices for these variables, one is the set of three separation distances, another the triad ($`r_<`$, $`r_>`$, $`\gamma `$) as in Eq. (1). With respect to a space-fixed “laboratory frame,” three more angles constitute the full set of six coordinates, the choice in Eq. (2) of ($`\theta `$, $`\theta ^{}`$, $`\varphi \varphi ^{}`$) being suited to the spherical polar coordinates of the individual vectors; thus, $`𝐱`$: ($`r\mathrm{sin}\theta \mathrm{cos}\varphi `$, $`r\mathrm{sin}\theta \mathrm{sin}\varphi `$, $`r\mathrm{cos}\theta `$).

In this paper, we present an alternative expansion to Eq. (1) based on cylindrical (azimuthal) symmetry which should be of wide interest in physics and astrophysics. This expansion arose in a recent investigation by two of us ref3  of gravitational potentials in circular cylindrical coordinates $`𝐱`$ = ($`R`$, $`\varphi `$, $`z`$):

$$\frac{1}{|𝐱𝐱^{}|}=\frac{1}{\pi \sqrt{RR^{}}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}Q_{m\frac{1}{2}}(\chi )e^{im(\varphi \varphi ^{})},$$
(3)
with $`Q_{m\frac{1}{2}}`$ a Legendre function of the second kind of half-integer degree ref4 , and $`\chi `$ defined as

$$\chi \frac{R^2+R^^2+(zz^{})^2}{2RR^{}}=\frac{r^2+r_{}^{}{}_{}{}^{2}2rr^{}\mathrm{cos}\theta \mathrm{cos}\theta ^{}}{2rr^{}\mathrm{sin}\theta \mathrm{sin}\theta ^{}}.$$
(4)
Although this expansion has been recorded in many places ref5 ; ref6 , its full significance has not been appreciated. We have traced its earliest occurrence to the work of E. Heine in the mid-nineteenth century ref5  and will, therefore, call it the “Heine identity.” At its most general, it takes the form

$$\frac{1}{\sqrt{v\mathrm{cos}\psi }}=\frac{\sqrt{2}}{\pi }\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}Q_{n\frac{1}{2}}(v)e^{in\psi },$$
(5)
reducing to Eq. (3) when applied to the inverse distance between two points. We have found it to be a more efficient approach for problems with a cylindrical geometry and have used it for compact numerical evaluation of gravitational potential fields of several axisymmetric and nonaxisymmetric mass distributions ref3 . We now set Eq. (3) in a broader context, together with new associated addition theorems and a novel application in quantum physics, hoping to encourage wider use of this expansion throughout physics.

## II The Alternative Expansions

The expansion in Eq. (1) disentangles the dynamics contained in the radial variables from symmetries, particularly under rotations and reflections, pertaining to the angle $`\gamma `$. Whereas the three variables {$`r_<,r_>,\gamma `$} at this stage are joint coordinates of $`𝐱`$ and $`𝐱^{}`$, a further disentangling in terms of the independent coordinates so as to handle permutational and rotational symmetry aspects of the problem is often useful and achieved through the addition theorem for spherical harmonics ref7 :

$`P_{\mathrm{}}(\mathrm{cos}\gamma )`$ $`=`$ $`{\displaystyle \frac{4\pi }{2\mathrm{}+1}}{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}Y_\mathrm{}m(\theta ,\varphi )Y_\mathrm{}m^{}(\theta ^{},\varphi ^{})`$ (6)
$`=`$ $`{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{\Gamma }(\mathrm{}m+1)}{\mathrm{\Gamma }(\mathrm{}+m+1)}}P_{\mathrm{}}^m(\mathrm{cos}\theta )P_{\mathrm{}}^m(\mathrm{cos}\theta ^{})e^{im(\varphi \varphi ^{})},`$
where $`Y_\mathrm{}m`$ are the usual spherical harmonics ref7 , $`\mathrm{\Gamma }(z)`$ is the gamma function, and $`P_{\mathrm{}}^m(z)`$ is the integer-order, integer-degree, associated Legendre function of the first kind ref2 . Using this to replace $`P_{\mathrm{}}(\mathrm{cos}\gamma )`$ in Eq. (1), we obtain the familiar Green’s function multipole expansion in terms of all six spherical polar coordinate variables $`𝐱`$ and $`𝐱^{}`$,

$$\frac{1}{|𝐱𝐱^{}|}=\frac{1}{\sqrt{rr^{}}}\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\left(\frac{r_<}{r_>}\right)^{\mathrm{}+\frac{1}{2}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\frac{\mathrm{\Gamma }(\mathrm{}m+1)}{\mathrm{\Gamma }(\mathrm{}+m+1)}P_{\mathrm{}}^m(\mathrm{cos}\theta )P_{\mathrm{}}^m(\mathrm{cos}\theta ^{})e^{im(\varphi \varphi ^{})}.$$
(7)
Apart from the first factor with dimension inverse-distance formed from the geometric mean of the two lengths $`r`$ and $`r^{}`$, this expression involves only four combinations: $`(r_</r_>,\theta ,\theta ^{},\varphi \varphi ^{})`$ of the six coordinates $`𝐱`$ and $`𝐱^{}`$. In spite of widespread familiarity with the multipole expansion, this reduction from six to four essential variables has not been appreciated fully. We were led to it by the parallel investigation below of the alternative expansion, and note that it is a natural consequence of the separation distance being independent of the orientation of that separation in the laboratory frame, and thus independent of two angles serving to specify that orientation.

This multipole expansion in terms of spherical harmonics is very broadly utilized across the physical sciences. For example, with $`\mathrm{}`$ and $`m`$ interpreted as the quantum numbers of orbital angular momentum and its azimuthal projection, respectively, a whole technology of Racah-Wigner or Clebsch-Gordan algebra is available ref8  for handling all angular (that is, geometrical or symmetry) aspects of an $`N`$-body problem, the dynamics being confined to radial matrix elements of the coefficients $`(r_</r_>)^{\mathrm{}+\frac{1}{2}}`$ in Eq. (1). Although many other systems of coordinates have been studied for problems with an underlying symmetry that is different from the spherical, Eq. (1), in combination with the addition theorem for spherical harmonics, has gained such prominence as to have become the Green’s function expansion of choice even for non-spherically symmetric situations.

But, consider now the expansion in Eq. (3), which may be viewed either, in analogy with Eq. (1), as an expansion in Legendre functions, now of the second kind in the joint variable $`\chi `$ of the whole system, with coefficients $`(RR^{})^{1/2}e^{im(\varphi \varphi ^{})}`$, or as a Fourier expansion in the variable $`(\varphi \varphi ^{})`$ with the $`Q`$’s as coefficients. In this latter view, a further step allows us to develop a new addition theorem for these Legendre functions. Interchanging the $`\mathrm{}`$ and $`m`$ summations (Fig. 1) in Eq. (7), we obtain

$$\frac{1}{|𝐱𝐱^{}|}=\frac{1}{\sqrt{rr^{}}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}e^{im(\varphi \varphi ^{})}\underset{\mathrm{}=|m|}{\overset{\mathrm{}}{}}\left(\frac{r_<}{r_>}\right)^{\mathrm{}+\frac{1}{2}}\frac{\mathrm{\Gamma }(\mathrm{}m+1)}{\mathrm{\Gamma }(\mathrm{}+m+1)}P_{\mathrm{}}^m(\mathrm{cos}\theta )P_{\mathrm{}}^m(\mathrm{cos}\theta ^{}).$$
(8)
Comparing with Eq. (3), we deduce:

$$Q_{m\frac{1}{2}}(\chi )=\pi \sqrt{\mathrm{sin}\theta \mathrm{sin}\theta ^{}}\underset{\mathrm{}=|m|}{\overset{\mathrm{}}{}}\left(\frac{r_<}{r_>}\right)^{\mathrm{}+\frac{1}{2}}\frac{\mathrm{\Gamma }(\mathrm{}m+1)}{\mathrm{\Gamma }(\mathrm{}+m+1)}P_{\mathrm{}}^m(\mathrm{cos}\theta )P_{\mathrm{}}^m(\mathrm{cos}\theta ^{}).$$
(9)
This is a new addition theorem for the Legendre function of the second kind. Note that $`Q_{m\frac{1}{2}}=Q_{m\frac{1}{2}}`$ as per Eq. (8.2.2) in ref2 .

Similarities and contrasts between the pairs of equations, Eqs. (1) and (7) and Eqs. (3) and (8), are worth emphasizing. Of the four variables $`(r_</r_>,\theta ,\theta ^{},\varphi \varphi ^{})`$, the first pair of equations, Eqs. (1) and (7), expresses the inverse distance as a series in powers of the first variable with coefficients Legendre polynomials of the first kind in $`\gamma `$, which is a composite of the other three variables and decomposable in terms of them through the addition theorem as in Eq. (6). The second pair of equations, Eqs. (3) and (8), on the other hand, expands in Eq. (3) the inverse distance in terms of the variable $`\varphi \varphi ^{}`$, with expansion coefficients Legendre functions of the second kind in $`\chi `$, which is a composite of the other three variables $`(r_</r_>,\theta ,\theta ^{})`$ and decomposable in terms of them through the addition theorem in Eq. (9). For this comparison, it is useful to recast Eq. (1) in the more suggestive form,

$$\frac{1}{|𝐱𝐱^{}|}=\frac{1}{\sqrt{rr^{}}}\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}P_{\mathrm{}}(\mathrm{cos}\gamma )e^{(\mathrm{}+\frac{1}{2})(\mathrm{ln}r_>\mathrm{ln}r_<)}.$$
(10)
Whereas this expansion has half-integers in the exponents and integer degree Legendre polynomials of the first kind as coefficients, Eq. (3) has integer $`m`$’s in the exponents and half-integer degree Legendre functions of the second kind as coefficients.

Yet another alternative to Eq. (3) follows upon casting the square root in the expression for the distance in terms of $`r,r^{},`$ and $`\gamma `$ in the form Eq. (5) through the definition

$$v\frac{1}{2}\left(\frac{r_<}{r_>}+\frac{r_>}{r_<}\right)=\frac{r^2+r_{}^{}{}_{}{}^{2}}{2rr^{}}.$$
(11)
This gives the expression

$$\frac{1}{|𝐱𝐱^{}|}=\frac{1}{\pi \sqrt{rr^{}}}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}Q_{n\frac{1}{2}}(v)e^{in\gamma },$$
(12)
now a Fourier expansion in $`\gamma `$ instead of the angle $`(\varphi \varphi ^{})`$ of Eq. (3), with $`Q_{n\frac{1}{2}}(v)`$ as the coefficients. In terms of hyperspherical coordinates, widely used in atomic and nuclear study of three (or more) bodies ref9 , the variable $`v`$ is $`\mathrm{csc}2\alpha `$, $`\alpha `$ being a “hyperangle”. In an Appendix, we present a number of alternative expressions for the functions $`Q_{m\frac{1}{2}}`$, some of them new, which we have found useful while working with the expansions in Eqs. (3) and (12).

## III Two-electron integrals

One familiar application of expressions for the inverse distance is to the Coulomb interaction between two charges. We contrast usage of the alternative expansions in Eqs. (1) and (3) for calculating the electrostatic interaction as it appears in atomic, molecular and condensed matter physics. Thus, consider the so-called “direct” part of this interaction between two electrons in the $`3d^2`$ configuration,

$$V_{ee}^D=𝑑𝐱𝑑𝐱^{}\psi _{3d}^{}(𝐱)\psi _{3d}^{}(𝐱^{})\frac{1}{|𝐱𝐱^{}|}\psi _{3d}(𝐱)\psi _{3d}(𝐱^{}).$$
(13)
The standard treatment ref10  uses Eqs. (1) and (8), carries out all the angular integrals through Racah-Wigner algebra, leaving behind radial “Slater integrals” $`F^k(dd),k=0,2,4,`$ and yielding (for illustrative purposes, all $`m`$ values have been set equal to zero)

$$V_{ee}^D=F^0(dd)+(4/49)F^2(dd)+(36/441)F^4(dd),$$
(14)
where the coefficients are evaluated in terms of Wigner $`3j`$-symbols or are available in tables ref10 . The Slater integrals,

$$F^k(dd)=r^2𝑑rr^2𝑑r^{}(r_<^k/r_>^{k+1})R_{3d}^2(r)R_{3d}^2(r^{}),$$
(15)
remain for numerical evaluation. In this example, upon evaluation with hydrogenic radial functions, we obtain $`V=0.092172`$ in atomic units.

The alternative calculation through Eq. (3) involves only the $`m=0`$ term and thereby the integral

$$V_{ee}^D=\frac{1}{\pi }R^{1/2}𝑑RR^{1/2}𝑑R^{}𝑑z𝑑z^{}Q_{\frac{1}{2}}(\chi )|\psi _{3d}(𝐱)|^2|\psi _{3d}(𝐱^{})|^2.$$
(16)
The integrand is a function of $`z`$ and $`R`$ variables alone and our numerical evaluation of this integral reproduces the value cited above.

As a second example, we computed an exchange integral for the $`3d4f`$ configuration again setting, for simplicity, all $`m`$ equal to zero:

$$V_{ee}^E=𝑑𝐱𝑑𝐱^{}\psi _{3d}^{}(𝐱)\psi _{4f}^{}(𝐱^{})\frac{1}{|𝐱𝐱^{}|}\psi _{3d}(𝐱^{})\psi _{4f}(𝐱).$$
(17)
The standard method through exchange Slater integrals $`G`$ and Wigner coefficients gives ref10

$$V_{ee}^E=(9/35)G^1(df)+(16/315)G^3(df)+(500/7623)G^5(df)$$
(18)
and, again through hydrogenic radial functions, gives the value $`V_{ee}^E=0.0082862`$. We reproduce the same result upon directly computing Eq. (17) with Eq. (3), again involving a single four-dimensional integral as in Eq. (16) with $`Q_{\frac{1}{2}}`$.

As the orbital angular momenta involved of the two electrons increase, the number of terms in expressions such as Eqs. (14) and (18) also grows, necessitating the computing of more Wigner coefficients and Slater integrals. By contrast, only a single term of the expansion in Eq. (3) and a single integral is necessary in our suggested alternative, the $`\varphi `$ integrations setting $`m=0`$ for direct terms and $`m`$ equal to the difference in the $`m`$ values of the two orbitals for exchange terms. This same selection rule imposed by the $`\varphi `$ integrations means that even in a calculation with several configurations and the imposition of antisymmetrization, such as in a multi-configuration Hartree-Fock scheme, matrix elements of $`|𝐱𝐱^{}|^1`$ between each term in the bra and in the ket gets a contribution from only one $`m`$ value in the expansion in Eq. (3). This is a significant economy.

Although the evaluation of the four-dimensional integrals in Eq. (16) is computationally more demanding than the two-dimensional integrals of Eq. (15), the preparatory work of the Wigner–Racah algebra is avoided. A more efficient alternative to directly evaluating the four-dimensional integrals in expressions such as Eq. (16) is to use the fact that they can be viewed as electrostatic interaction energies between charge densities given by the product of wave functions. The charge density $`\rho (𝐱^{})`$ is then used to find the potential at $`𝐱`$ and the resulting potential integrated with the appropriate charge density $`\rho (𝐱)`$ to compute the integral ref11 . This same approach applies to the decomposition according to Eq. (3) for each $`m`$-th component of the potential. A two-dimensional Poisson equation in the primed variables in Eq. (16) is first solved and then Eqs. (16) and (17) evaluated as two-dimensional integrals over the unprimed variables ref12 . As already noted, direct integrals will involve only the $`m=0`$ component; exchange ones involve a single $`m`$ value equal to the difference in the azimuthal quantum numbers of the two orbitals.

## IV Summary

The inverse distance between two points $`𝐱^{}`$ and $`𝐱`$ is intimately involved in Coulomb and gravitational problems. Its expansion in terms of Legendre polynomials $`P_{\mathrm{}}`$ of the angle between the vector pair or a further double-summation expansion involving the individual polar angles of the vectors are well known and widely used in physics and astronomy. We have discussed an alternative in terms of cylindrical coordinates, a single summation in terms of Legendre functions $`Q_{m\frac{1}{2}}`$ of the second kind in a pair variable $`\chi `$ or double summations involving the individual coordinates. These expansions are better suited to problems that are decomposable in azimuthal symmetry as shown by applications in ref3  and by an illustration here for very common electron-electron calculations throughout many-electron physics. Further variants are possible for other coordinates such as ring or toroidal, parabolic, bispherical, cyclidic and spheroidal ref13 , and we hope to return to them in future publications. Connections to the theory of Lie groups will also be of interest ref14 .

###### Acknowledgements.

We thank Prof. Ken Schafer for his remark that better ways of evaluating two-electron exchange integrals would be useful. ARPR thanks the Alexander von Humboldt Stiftung and Profs. J. Hinze and F. H. M. Faisal of the University of Bielefeld for their hospitality during the course of this work. This work has been supported, in part, by NSF grant AST-9987344 and NASA grant NAG5-8497. Work performed at LLNL is supported by the DOE under Contract W7405-ENG-48.

*

## Appendix A Alternative expressions

We present in this appendix a number of alternative expressions for the functions $`Q_{m\frac{1}{2}}`$ which are useful in calculations using expansions such as Eqs. (3) and (13). Setting $`\theta =\theta ^{}=\pi /2`$ in Eq. (9) and using Eq. (8.756.1) of ref15  gives

$$Q_{m\frac{1}{2}}(v)=\pi \underset{\mathrm{}=|m|}{\overset{\mathrm{}}{}}\left(\frac{r_<}{r_>}\right)^{\mathrm{}+\frac{1}{2}}\frac{\mathrm{\Gamma }(\mathrm{}m+1)}{\mathrm{\Gamma }(\mathrm{}+m+1)}\frac{\pi 2^{2m}}{[\mathrm{\Gamma }(1+\frac{\mathrm{}m}{2})\mathrm{\Gamma }(\frac{1\mathrm{}m}{2})]^2},$$
(19)
which can be rewritten as

$$Q_{m\frac{1}{2}}(v)=\pi e^{(m+\frac{1}{2})\eta }\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}2^{12m4\mathrm{}}\left(\begin{array}{c}2\mathrm{}\\ \mathrm{}\end{array}\right)\left(\begin{array}{c}2\mathrm{}+2m1\\ \mathrm{}+m\end{array}\right)e^{2\mathrm{}\eta },$$
(20)
where we have defined $`r_</r_>e^\eta `$, $`v=\mathrm{cosh}\eta `$. Although the $`\mathrm{}`$-th term of these series is in different form from what one obtains through the more familiar formula for $`Q`$ as a hypergeometric function ref16 , namely,

$$Q_{m\frac{1}{2}}(v=\mathrm{cosh}\eta )=\frac{\sqrt{\pi }\mathrm{\Gamma }(m+\frac{1}{2})}{\mathrm{\Gamma }(m+1)}e_2^{(m+\frac{1}{2})\eta }F_1(\frac{1}{2},m+\frac{1}{2};m+1;e^{2\eta }),$$
(21)
their equivalence follows from straightforward algebra. Also, another standard expansion for $`Q`$ in powers of $`v`$ as in Eq. (8.1.3) of ref2 ,

$$Q_{\nu \frac{1}{2}}(v)=\frac{\sqrt{\pi }\mathrm{\Gamma }(\nu +\frac{1}{2})}{\mathrm{\Gamma }(\nu +1)}(2v)_2^{\nu \frac{1}{2}}F_1(\frac{\nu }{2}+\frac{3}{4},\frac{\nu }{2}+\frac{1}{4};\nu +1;\frac{1}{v^2}),$$
(22)
is equivalent. However, the results directly in powers of $`r_</r_>`$ in Eqs. (19), (20), and (21) are more convenient in many applications. Among specific features worth noting in these alternative expansions are that only even powers of $`(r_</r_>)`$ occur in the sum in Eq. (20) and that for any $`m`$, the sum in Eq. (19) runs over all $`\mathrm{}`$ values compatible with it, $`\mathrm{}|m|`$, as per their interpretation as angular momentum quantum numbers.

In the multipole expansion in Eq. (1), $`\gamma `$ is an angle formed out of the set $`(\theta ,\theta ^{},\varphi \varphi ^{})`$ and, therefore, $`\mathrm{cos}\gamma `$ in the functions $`P_{\mathrm{}}`$ has range of variation from –1 to 1. On the other hand, in the expansions in Eqs. (3) and (5), the arguments $`v`$ and $`\chi `$ of the Legendre functions of the second kind range from 1 to infinity and, therefore, can be written in terms of hyperbolic functions as $`\mathrm{cosh}\eta `$ and $`\mathrm{cosh}\xi `$, respectively. From Eq. (4) we have the link between them,

$$\mathrm{cosh}\eta =\mathrm{cos}\theta \mathrm{cos}\theta ^{}+\mathrm{sin}\theta \mathrm{sin}\theta ^{}\mathrm{cosh}\xi .$$
(23)
This disentanglement of $`v`$ (or $`\eta `$) in terms of a triad is the counterpart of Eq. (2) and may be used with addition theorems given in the literature such as ref5 ; ref6

$$Q_{m\frac{1}{2}}(\mathrm{cosh}\eta )=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}(1)^n\frac{\mathrm{\Gamma }(mn\frac{1}{2})}{\mathrm{\Gamma }(n+m\frac{1}{2})}Q_{m\frac{1}{2}}^n(\mathrm{cos}\theta )P_{m\frac{1}{2}}^n(\mathrm{cos}\theta ^{})e^{n\xi }.$$
(24)
An alternative to writing $`\chi `$ as $`\mathrm{cosh}\xi `$ is to set $`\chi =\mathrm{coth}\zeta `$ in Eq. (23). In that case, it follows that

$$\mathrm{exp}\zeta \sqrt{\frac{\mathrm{cosh}\eta \mathrm{cos}(\theta +\theta ^{})}{\mathrm{cosh}\eta \mathrm{cos}(\theta \theta ^{})}}.$$
(25)
Other variants of the expansions in Eqs. (9) and (24) follow from identities satisfied by the Legendre functions $`P`$ and $`Q`$. In particular, there is an interesting pair of relations involving index interchange given in Eqs. (8.2.7) and (8.2.8) of ref2 :

$`Q_{n\frac{1}{2}}^m(\mathrm{cosh}\eta )=(1)^m\sqrt{{\displaystyle \frac{\pi }{2\mathrm{sinh}\eta }}}\mathrm{\Gamma }(mn+{\displaystyle \frac{1}{2}})P_{m\frac{1}{2}}^n(\mathrm{coth}\eta ),`$ (26)
$`Q_{m\frac{1}{2}}^n(\mathrm{coth}\eta )=(1)^m\sqrt{{\displaystyle \frac{\pi \mathrm{sinh}\eta }{2}}}{\displaystyle \frac{\pi }{\mathrm{\Gamma }(mn+\frac{1}{2})}}P_{n\frac{1}{2}}^m(\mathrm{cosh}\eta ).`$ (27)


format: text

text:
# Electron confinement and optical enhancement in Si/SiO2 superlattices

## I Introduction

Informatics and telecommunications are the driving forces behind the development of new photonic devices.lampro  Such devices can be assembled from a wide variety of materials; among these, silicon is probably the most prevalent, although it suffers from having an indirect band gap. Direct band-gap materials such as GaAs, InP, and many other III-V and II-VI semiconductors, are efficient electron-photon energy converters. However they are costly for large-scale applications, and would be less interesting if full Si-based photonics were possible. These considerations have lead to considerable efforts towards the development of new types of Si-based materials having direct energy gap in the visible (viz. 1.5 to 3 eV),Zhang  suitable for photo-diodes and lasers, as well as solar-cell applications.

The simultaneous discovery by CanhamCanham  and Lehmann and GöseleLehmannGosele  of intense luminescence in porous silicon ($`\pi `$-Si) has opened up new horizons for Si-based materials. $`\pi `$-Si can be regarded as consisting of long and thin nanowires,ChuangCollinsSmith  forcing the electronic states to be confined within the finite dimension of the nanostructures. The shift of the luminescence spectrum towards the blue with decreasing nanowire sizeCanham  suggests that the confinement of the electronic states might be responsible for the visible luminescence. Unfortunately, most forms of $`\pi `$-Si are unstable, lasting only a few hours before a strong decrease of the luminescence efficiency is observed. This instability has sometimes been attributed to the volatile Si–H bonds at the surface.LehmannGosele  Annealing under nitrogen and oxygen has been suggested as a means of stabilizing the structures; however, the annealed material has weaker luminescence than the parent $`\pi `$-Si. Layered polysilanes (Si6nH6) also exhibit optical properties similar to $`\pi `$-Si. Uehara et alUehara  have performed first-principles calculations of layered polysilanes structures and showed that porous Si may be akin to thin Si6H6 sheets.

As a perhaps more promising avenue, confinement can also be achieved in superlattices (SLs) — or quantum wells — as well as quantum wires and quantum dots. Lu et al have reported enhanced luminescence in the visible part of the spectrum of Si/SiO2 SLs fabricated by molecular beam epitaxy (MBE),LuLockwoodBaribeau  and a blue shift was observed when the thickness of the Si layers was reduced from 6 to 2 nm. In contrast to $`\pi `$-Si, Si/SiO2 SLs present the advantage of being stable, with no significant decrease of the luminescence observed with time.

Motivated by these promising experimental observations, and to provide insight into the microscopic physics associated with the luminescence efficiency of SL-based devices, we undertook a detailed, first-principles investigation of the electronic and optical properties of Si/SiO2 SLs. Silicon and SiO2 are both already standard components of MOSFETs and other devices. SiO2 has an energy gap of $``$9 eV; Si has an indirect gap of $``$1.1 eV. Thus, the electrons are confined by the SiO2 barriers within the silicon quantum wells. Luminescence, however, is not determined solely by confinement since the nature of the energy gap is also important. The possibility of creating a direct energy gap in a SL can be discussed in simplified terms within the concept of Brillouin zone folding:gnutzmann  The periodicity of a SL, $`d`$, being larger than that of the constituent lattices, $`a_0`$, the Brillouin zone gets folded to dimension $`\pi /d`$ instead of $`\pi /a_0`$, which may bring the indirect minimum of the conduction band to the $`\mathrm{\Gamma }`$-point and thus produce an effectively direct energy-gap material. “Quasi-direct” energy gaps were so obtained in Si/Ge superlattices;pearsall  however, there was little enhancement in the optical properties.

In this work, we examine the luminescence properties of Si/SiO2 SLs using model structures and first-principles methodology. Transmission electron microscopy reveals that in Si/SiO2 SLs, both Si and SiO2 layers are amorphous.LuLockwoodBaribeau  Amorphous structures are however not easily amenable to quantum calculations. Here, we assume crystalline phases for both the Si and SiO2 slabs. This is further motivated by the need to study simple ideal systems as benchmarks, and in particular to determine how the difference in the interface affects the optical properties. We however note that crystalline-Si based SLs have very recently been fabricated.Lu\_private

Two simple model structures for Si/SiO2 SLs are considered: (i) The double-bonded model (DBM) of Herman and Batra,HermanBatra  shown in Fig. 1(a); here the Si and SiO2 layers are arranged so as to minimize the lattice mismatch and have relatively high symmetry, and the dangling bonds at the interfaces are saturated by the addition of a single double-bonded oxygen atom. (ii) A bridge-oxygen model (BOM), shown in Fig. 1(b), where the interface dangling bonds are saturated by adding an oxygen atom to bridge two Si bonds. Full details of the models are given in Section II.1.

The electronic properties of both the DBM and the BOM have been studied already within an empirical tight-binding framework by Tit and Dharma-wardana,TitDharma  who observed the bands along the high-symmetry lines of the BZ parallel to the $`k_z`$ axis to be essentially dispersionless, in agreement with the present calculations. However, the band gaps were found to be direct in both models, in contradiction with our results. The DBM has been studied from first principles by Punkkinen et al.Punkkinen  using a full-potential linear-muffin-tin orbital approach. However, no other interface-structure was examined and thus the effect of changing the interface could not be assessed; further, the optical properties were not calculated. Kageshima and Shiraishi have also used first-principles methods to examine different models for the Si–SiO2 interface in oxyde-covered Si slabs,Kageshima  but the SL geometry was not examined. Their calculations indicate that the interface plays an important role in the light-emitting properties of the system and that, in particular, interfacial Si-OH bonds are possibly related to the observed luminescence. Here we present results for the electronic structure, the dielectric function and the absorption coefficient of both the DBM and the BOM, and focus on a comparison between the two models.

## II Computational details

### II.1 Model structures

SLs can be fabricated by MBE;HermanSitter  a major advantage of this technique, as compared to others, is that it gives sharper interfaces. In the case of Si/SiO2 SLs, Lu et al.LuLockwoodBaribeau  report an interface thickness of about 0.5 nm, while the thickness of the SiO2 and Si layers vary between 2 and 6 nm. Here we assume that the interfaces are perfectly sharp.

The model Si/SiO2 SL structures are constructed by alternating layers of crystalline (diamond) silicon and crystalline SiO2 in the ideal $`\beta `$-cristobalite structure, viz. with a Si–O–Si bond angle of 180 (space group Fd$`\overline{3}`$m). The (experimental) lattice constants are 5.43 and 7.16 Å, respectively, with therefore a lattice mismatch of 32%. The mismatch can be reduced to less than 7% by rotating the Si unit cell by an angle of $`\pi /4`$ so that the diagonal of the Si-diamond structure fits, approximately, the cubic edge of the $`\beta `$-cristobalite unit cell. The superposition of the Si and SiO2 layers gives rise to dangling bonds on some of the interface Si atoms. These can be satisfied in different ways. In the DBM, this is done by adding a single O atom onto one of the Si atoms. The resulting unit cell has 23 atoms \[cf. Fig.1(a)\] and dimensions 5.43 Å in the $`xy`$ plane and 5.43 $`\times `$ (1 + $`\sqrt{2}`$) = 13.11 Å in the $`z`$ (growth) direction. In this model, the Si–Si and Si–O distances, 2.35 and 1.66 Å, respectively, are realistic, but the double-bonded oxygen Si=O is not found in naturally occuring silicates. In the BOM, Fig. 1(b), an oxygen atom saturates two dangling bonds provided by two different Si atoms; this model contains 21 atoms. Here the Si–O–Si angle is the usual 144 and the length of the Si–O bonds at the interface is 2.02 Å, i.e., somewhat longer than found in nature. No energy relaxation of the models have been carried out. However, as will be discussed in section III.3, we found gap states to be present in the BOM; in order to understand this, we investigated two other variants of the BOM, having Si–O–Si angles of 109 and 158, respectively.

The Brillouin zone for the SL, along with that for c-Si, is depicted in Figure 2. The principal symmetry axes are indicated: the $`Z\mathrm{\Gamma }`$, $`XR`$ and $`MA`$ directions correspond to possible growth directions where confinement is expected to occur, as discussed below.

### II.2 First-principles calculations

The electronic-structure calculations were carried out with the WIEN97 code,WIENref  which uses the all-electron, full-potential, linearized-augmented-plane-wave (LAPW) method,Singh  within the framework of density-functional theoryHohenbergKohn ; KohnSham  in the local-density approximation (LDA).Payne  Thus, the core states in the (spherical) atomic region are described using an atomic basis, while the delocalized states in the interstitial region are expanded in plane waves which are correctly formulated to account for the core region. The Kohn-Sham (KS) basis inside the atomic spheres is expressed as an angular-momentum expansion:

$$\psi (\stackrel{}{k}_n,\stackrel{}{r})=\underset{l,m}{}\left[A_{lm}^nu_l(r,E_l)+B_{lm}^n\frac{u_l}{E}|_{E=E_l}\right]Y_{lm}(\widehat{r}).$$

Here the energy $`E_l`$ is calculated and fixed at the first cycle of the self-consistent-field calculation. The KS equations are solved for a grid of $`\stackrel{}{k}`$ points as discussed below.

The integration of the radial secular equations was performed on a radial mesh containing 581 points. The radius of the spheres for silicon and oxygen atoms were set at 0.87 and 0.79 Å, respectively, and an energy cutoff of $`7.1`$ Ry was used for separating core from valence electrons. This generated 11 149 plane waves and the matrix size for the eigenvalue problem contained 5581 elements. These two parameters correspond to an angular momentum cutoff of 7.5 (see SinghSingh  for further details), a value which ensures proper energy convergence of the solution.

The DBM has nine non-equivalent atoms while the BOM has eight. The irreducible wedge of the Brillouin zone is obtained after four symmetry operations, which reduces considerably the computational load. The supercell in the $`z`$ direction being 3 times longer than in the $`x`$ or $`y`$ directions, only 181 $`\stackrel{}{k}`$ points were necessary to achieve convergence of the electronic energies. However, in order to ensure convergence of the wave functions needed to set up the optical matrix, the density of $`\stackrel{}{k}`$ points had to be doubled. Integration of the Brillouin zone was performed using the tetrahedron method and Broyden density mixing. The self-consistent calculations were performed in parallel mode (using four processors) on a Silicon Graphics Origin 2000 computer system. About 10 cycles were necessary to get convergence of the energy to about 10-4 Ry.

## III Results and discussion

### III.1 Band structures

The electronic band structures along the high symmetry axes of the tetragonal Brillouin zone (cf. Fig.2) for the two models are plotted in Fig. 3, together with the densities of states (DOS). For the BOM, unless otherwise noted, the interfacial Si-O-Si angle is the normal 144. The confinement in the $`z`$-direction can already be inferred from the essentially dispersionless character of the bands along the three symmetry axes which are parallel to the axis of the SL, namely $`XR`$, $`Z\mathrm{\Gamma }`$ and $`MA`$; this was also observed by Punkkinen et al.Punkkinen  Thus, optical transitions between bands lying in the growth direction are favored by the confinement.

The DOS of the two models differ strongly in the region of the gap, with an isolated conduction band clearly visible in the gap of the BOM which is absent in the DBM. The only difference between the two models lies in the way that the interfacial dangling bonds are fixed; the gap state in the BOM is thus connected to the bridge-bonded Si atom. The results, evidently, are sensitive to the way that the interface is formed, either in experiment, or in a theoretical model. This question will be discussed in more detail in section III.3.

The effect of the confinement in both models is clearly visible from the DOS. If we ignore the gap state in the BOM for the moment, the DOS at the Fermi level has a step-function-like threshold, as opposed to bulk-silicon which is knownCardona  to have a threshold depending on the photon energy $`\mathrm{}\omega `$ as $`(\mathrm{}\omega E_g)^{1/2}`$, where $`E_g`$ is the gap energy. Higher transition probabilities are thus expected in Si/SiO2 SLs.

Figure 4 zooms on the band structure of the two models near the Fermi level. The numerical values of the various possible transition energies are given in Table 1. For the DBM, the valence-band maximum (VBM) is at the $`\mathrm{\Gamma }`$ point but the conduction band minimum (CBM) is in between the $`\mathrm{\Gamma }`$ and $`X`$ points (labeled $`\mathrm{\Gamma }X/2`$ in Table 1), giving an indirect gap of 0.81 eV. For the BOM, the VBM is at $`X`$ while the CBM is at $`Z`$, giving an almost null (0.01 eV) indirect band gap. However, it should be remembered that the energy gaps are significantly underestimated in the LDA and other implementations of the DFT. Thus the bulk-Si LDA energy gap is 0.49 eV, about 0.6 eV below the true value of 1.1 eV. Better (but non-rigorous) estimates for the DBM and the BOM energy gaps would thus be 1.4 and 0.6 eV, respectively. The unusually low value for the energy gap of the BOM is an indication that the lowest conduction band is really a gap state. If this state \[“CB1” in Fig.4(b)\] is assumed to be an unphysical artifact of the model, and temporarily ignored, the indirect gap of 1.49 eV gets corrected to $``$2.1 eV. Thus, the energy gap would be in the near infrared or in the red, consistent with experiment.LuLockwoodBaribeau ; Mulloni ; Novikov

The dispersionless character of the bands in the growth direction provides evidence of the confinement effect, as noted above. However, confinement alone is not sufficient to explain the enhancement of the optical properties; comparison of the direct and indirect transitions for the two models can clarify this question. For the DBM, the direct transitions at $`\mathrm{\Gamma }`$ and $`\mathrm{\Gamma }X/2`$ have energies of 0.98 eV and 0.84 eV, respectively, while the $`\mathrm{\Gamma }\mathrm{\Gamma }X/2`$ indirect transition (i.e., the band gap) costs 0.81 eV. Thus the smallest direct–indirect gap (SDIG) in the LDA for the DBM is a mere 0.03 eV ($`0.840.81`$ eV). For the BOM (if the “CB1” gap state is ignored), the direct transitions at $`X`$ and $`Z`$ have energies of 2.19 and 1.87 eV, respectively, compared to 1.49 eV for the $`Z`$$`X`$ indirect transition, for a SDIG of only 0.35 eV. These numbers can be compared to the corresponding ones for bulk silicon (in the LDA): the direct transition at the $`\mathrm{\Gamma }`$ point has energy 2.52 eV while the direct transition at the CBM (near $`3/4\mathrm{\Gamma }X`$) costs 3.20 eV. The indirect band gap is 0.49 eV, so that the SDIG in this case is $``$2.0 eV. This is much larger than that for the DBM and the BOM. The latter, as a consequence, has much better optical properties than bulk Si. It is also clear that the DBM will have better optical properties than the BOM (see section III.4) since the SDIG of the DBM is smaller than that of the BOM.

### III.2 Comparison with tight-binding calculations

It is instructive to compare the present results to those obtained for the same model systems by Tit and Dharma-wardanaTitDharmaDBMBOM  within a tight-binding (TB) framework. In a TB description,Papaconstantopoulos  the wavefunctions are expressed as linear combinations of atomic orbitals and the exact many-body Hamiltonian is replaced by a parametrized Hamiltonian matrix whose eigenvalues and eigenvectors yield the energies and the wavefunctions of the corresponding electron levels. The size of the numerical problem is determined by the number of atomic orbitals chosen to describe the valence electrons, and about 4 or 5 orbitals per atom are typically used. This allows the study of much larger systems than is possible in first-principles calculations. It is therefore important to assess the validity of TB models in order to re-calibrate and improve the TB parametrization. Indeed, the parameters are typically fitted to experimental data or first-principles calculations for the pure crystalline phases. For mixed systems, e.g., a SL, the parameters are assumed to be transferable except for simple adjustments of energy reference or scaling to bond lengths, etc.

An important parameter entering TB calculations is the valence-band offset (VBO), which is the difference in energy between the VBM of the two materials constituting the SL.Williams  In their calculations, Tit and Dharma-wardana assumed a VBO of 3.75 eV, based on the experimental offset between Si and amorphous SiO2. However, this value need not be appropriate to the idealized DBM and BOM models.

In Fig. 5(a) we compare the LDA band structure for the DBM with the TB results of Tit and Dharma-wardana,TitDharmaDBMBOM  who used a VBO value of 3.75 eV; Fig. 5(c) is the corresponding plot for the BOM. The agreement with this particular value of the VBO is, at best, qualitative; this may indicate that the choice of VBO needs revision. In order to assess this, we have done new TB calculations using different VBO values. By adjusting as closely as possible the valence band structures to the TB ones, better values of the VBO were found to be 0.0 eV for the DBM and 1.0 eV for the BOM; the results are indicated in Fig. 5(b) and (d), where the band gaps from LDA are adjusted to the TB values to help comparison. For the BOM, the gap state “CB1” from LDA coincides well with the first conduction band obtained from TB, except in the $`X`$$`Z`$ region of the BZ. The first-principles “CB2”, in contrast, is somewhat higher than found from the TB calculations. Thus, while not perfect, the agreement with the first-principles calculations improves significantly with these new VBO values.

### III.3 Gap states in the BOM

In view of the presence of a gap state in the BOM — labeled “CB1” in Fig. 4 — and in order to assess the role of the interface on the electronic properties, other possibilities for the interfacial Si–O–Si angle were considered: (i) 109, corresponding to positioning the oxygen atom on a normal silicon site of the Si lattice; here, $`d_{\mathrm{Si}\mathrm{O}}=2.35`$ Å; (ii) 144 which corresponds to the experimental value of the Si–O–Si angle, used in the calculations discussed above, yielding $`d_{\mathrm{Si}\mathrm{O}}=2.02`$ Å; and finally (iii) 158, a value obtained by relaxing the oxygen atom position, and which corresponds to a (local) energy minimum; in this case, $`d_{\mathrm{Si}\mathrm{O}}=1.96`$ Å. The “normal” Si–O distance in silica is 1.61 Å, which cannot be accomodated by the crystalline silicon lattice in the BOM.

Figure 6 shows the band structures for the three cases; here we consider only the $`XRZ`$ direction, wherein lies the energy gap. Evidently, the precise value of the interfacial Si–O–Si bond angle, and corresponding Si–O bondlength, have a sizeable effect on the band structure. In fact, the 158–BOM is, within the LDA, a metal if CB1 is not assumed to be a gap state, as discussed earlier. Likewise, the 144–BOM is nearly metallic. These results indicate that the acceptable range of interfacial Si–O–Si angles, and the consequent longer than “normal” Si–O distances, are probably essential issues in explaining the luminescence in Si/SiO2 SLs.

### III.4 Optical properties

The Kohn-Sham calculations provide matrix elements and joint densities of states necessary for the calculation of the complex dielectric function $`\stackrel{}{ϵ}=\stackrel{}{ϵ}_r+i\stackrel{}{ϵ}_i`$. The absorption coefficient $`\alpha `$, which can then be deduced from $`\stackrel{}{ϵ}`$ as a function of the photon energy (see below), provides a detailed picture of one aspect of the optical properties of the material. In principle, the luminescence requires a knowledge of the excited states of the system, with electrons occupying the conduction band. Such states are not available from DFT calculations. Further, the electrons in the conduction band are associated with holes in the valence band. For low carrier concentrations, the screening is weak and hence exciton formation occurs.knox  This many-body effect is also required for a complete description of luminescence. Such calculations require, e.g., the solution of the Bethe-Salpeter equation for the electron-hole pair (see Chang et al., Ref. Chang, ) and are still prohibitive even for our model structures. Excitonic effects will thus be neglected here.

As already noted, the band gaps are underestimated by LDA-DFT and this has to be corrected if the absorption thresholds are to be realistic. In practice, this can be done to a reasonable approximation by rigidly shifting all the conduction bands to the appropriate energy. However, we remain within the DFT framework and calculate the absorption in the Fermi-golden-rule approximation. Also, we assume that the emission spectrum is similar to the absorption spectrum and neglect excitonic effects. In spite of these approximations, the main mechanisms of luminescence enhancement — the effects arising from the joint density of states and the matrix elements — would be correctly captured.

Using the program written by Abt et al.AbtDraxlKnoll  as part of the WIEN97 package,WIENref  the imaginary part of the dielectric function has been determined; it is given by:

$`ϵ_i^\alpha (\omega )=`$ $`\left({\displaystyle \frac{4\pi e^2}{m^2\omega ^2}}\right){\displaystyle \underset{v,c}{}}{\displaystyle \frac{2d\stackrel{}{k}^3}{(2\pi )^3}|<c\stackrel{}{k}|^\alpha |v\stackrel{}{k}>|^2}`$
$`<f_{v\stackrel{}{k}}(1f_{c\stackrel{}{k}})\delta (E_{c\stackrel{}{k}}E_{v\stackrel{}{k}}\mathrm{}\omega ),`$
where $`f_{v\stackrel{}{k}}`$ is the Fermi distribution and $`^\alpha `$ is the $`\alpha `$-component of the electron-radiation interaction Hamiltonian in the Coulomb gauge; it corresponds to the probability per unit volume for a transition of an electron in the valence band state $`|v\stackrel{}{k}>`$ to the conduction band state $`|c\stackrel{}{k}>`$ to occur. From Kramers-Kronig relations one then deduces the real part $`ϵ_r`$.Cardona  The dielectric function $`\stackrel{}{ϵ}`$ is defined as the square of the complex refractive index $`\stackrel{}{n}=n_r+in_i`$, so that

$$\alpha (E)=4\pi \frac{E}{hc}n_i=4\pi \frac{E}{hc}\left[\frac{(ϵ_r^2ϵ_i^r)^{1/2}ϵ_r}{2}\right]^{1/2},$$

with $`c`$ the speed of light in vacuum, $`h`$ Planck’s constant, and $`E`$ the photon energy.

Figure 7 shows the $`z`$-component of the absorption coefficients for the DBM and for the three different variations of the BOM. We also give, for comparison, the corresponding curves for silicon and for the ideal $`\beta `$-cristobalite SiO2 structures. The energy gaps are those directly from the LDA; correcting the gaps would only affect the position of the onset, but not the general aspect. The absorption coefficient for the two models are similar. This similarity is, to a large extent, a “zone folding” effect, as the two model structures have identical $`z`$-dimensions (the growth axis). However, closer inspection reveals that the DBM has better absorption properties than the three BOMs: the onset of absorption is sharper, and occurs earlier in the DBM than in all the BOM. This is consistent with the band structure analysis above that showed the SDIG to be smaller in the DBM than in the BOM ($``$ 0.03 vs $``$0.35 eV), implying a larger transition probability for the former than the latter.

In addition, the calculated SDIG for the BOM with an interfacial Si–O–Si angle of 109 is 0.18 eV — the band gap (or the indirect transition) being 1.53 eV while the direct transition between $`X`$ and $`Z`$ are 2.06 eV and 1.71 eV respectively — in between the SDIG for the DBM and the BOM with 144. This is in agreement with the absorption curves of Fig. 7, where the absorption onset is higher for the DBM, lower for the BOM with 109 and then even lower for the BOM with 144 or 158. The last two give similar absorption curves as can be seen from Fig. 7.

Thus, the four models (the DBM and the BOM with bond angles 109, 144 and 158), which differ only by their interface definition, give rather different optical properties. The onset of the absorption curve in the BOM with Si–O–Si angles of 144 and 158 do not differ that much from the one obtained from bulk-Si, while the DBM and the BOM with 109, viz. the BOM with higher Si–O bondlengths at the interface, give better optical properties. This is additional indication that the interfacial atomic structure has to be connected to the optical enhancement in the SLs, in agreement with Kageshima’s analysis.Kageshima  Likewise, it can be concluded from Fig. 7 that bulk Si has poor optical properties compared to the SLs, since the onset of absorption lags behind that for the two SLs and the SDIG is much larger.

## IV Concluding Remarks

We have used first-principles calculations to (i) study the electronic and optical properties of Si/SiO2 SL models and (ii) examine the applicability of tight-binding calculations for these systems.

Concerning the second point, our results indicate that, for both models, the VBO value of 3.75 eV used in TB calculations is excessive. In addition, the direct-transition nature of the gap from the TB calculations is incorrect. The present calculations indicate that the band gap in both models is indirect (albeit quasi-direct), in the red or infrared region of the spectrum. Confinement is further confirmed by the essentially dispersionless character of the electronic band structures in the growth direction (cf. Figs 3 and 4). Our calculations indicate, also, that the SLs have enhanced optical properties as compared to pure Si.

The influence on the optical properties of the Si–SiO2 interface has been assessed from the absorption coefficient. This quantity depends critically on the details of the interfacial structure. Realistic models are thus essential for determining the electronic and optical properties. We are presently studying improved models for this system.Tran

Acknowledgments – It is a pleasure to thank Gilles Abramovici, Ralf Meyer and Michel Côté for useful discussions. This work is supported by grants from the Natural Sciences and Engineering Research Council (NSERC) of Canada and the “Fonds pour la formation de chercheurs et l’aide à la recherche” (FCAR) of the Province of Québec. We are indebted to the “Réseau québécois de calcul de haute performance” (RQCHP) for generous allocations of computer resources.


text:
# Sub-wavelength lithography over extended areas

## I Introduction

Classically, to create an optical image, one has to modulate a wavefront of an electromagnetic wave in space. The minimum resolvable feature size of an imaged object corresponds, roughly speaking, to the minimum modulation period allowed, which turns to be of the order of the wavelength $`\lambda `$ of the light used. In fact, the Rayleigh criterion establishes that the best resolution that can be achieved classically is about $`\lambda /2`$, which is usually known as the diffraction limit.

When the quantum nature of light is considered, one is naturally confronted with the role that photon fluctuations play in setting fundamental performance limits for imaging systems. Even if all the technical noise sources are eliminated from the imaging system, the corpuscular nature of the photon induces fluctuations, or shot noise, that determine a seemingly fundamental spatial resolution, or standard quantum limit, of about $`\lambda /(2\sqrt{N})`$, where $`N`$ is the average number of photons. Sub-wavelength imaging have been used in number of applications Dragsten ; Flock ; Kamimura ; Jelles ; Denk ; Putman  and a careful analysis shows that indeed the shot-noise sets the resolution limit Denk ; Putman .

However, the quantum viewpoint allows for strategies that could significantly improve the spatial resolution beyond the standard quantum limit. A typical way of reducing photon-counting noise is by using multi-mode squeezed light Kolobov 1 ; Kolobov 2 ; Kolobov 3 ; Kolobov ; Fabre ; Fabre 2 . This possibility has been experimentally demonstrated in other precision measurement schemes Min  and allows one to attain an optimum spatial resolution proportional to $`\lambda /2N`$, usually known as the Heisenberg limit. These sub shot-noise imaging systems enable resolving, in principle, arbitrarily small details of an object in a diffraction-limited optical system.

However, writing images imposes even more stringent requirements than sub-wavelength resolution of spatial features because in image writing one wants to write small details with high contrast. A particular field where circumventing the classical resolution limit is becoming more and more important is optical lithography, which is the primary tool for writing electronic-circuit patterns. Current production technologies have tended to use light of shorter wavelengths to fabricate ever-smaller device features.

It has been known for some time that entangled photon pairs can be used to achieve Heisenberg limited resolution of time Bollinger ; Huelga ; Buzek  and phase Rarity ; Margolus ; Fonseca ; Trifonov , but only very recently it has been proposed the use of entanglement to increase high contrast image resolution indefinitely Boto ; Kok ; Bjork . The reason that these entangled quantum states show increased resolution can be traced back to the fact that they allow the modulation period to be as small as $`\lambda /(2N)`$, and thus they approach Heisenberg-limited resolution. The process can be envisioned as the photons clustering into a $`N`$-photon quasi-particle with a linear momentum $`N`$ times large as that of a single photon, and therefore with a shorter de Broglie wavelength Jacobson . It is the de Broglie wavelength that ultimately determines the interference resolution. This has been appreciated for a long time in atomic, molecular and solid-state physics, but has only recently been noticed for electromagnetic waves.

In an earlier paper we discussed the use of reciprocal binomial states in sub-wavelength resolution lithography Bjork . Our method works for even number of photons $`2N`$ and it is especially germane to determine the exposure sequence to generate any pixellated pattern on a $`(N+1)\times (N+1)`$ grid, occupying a square with a half wavelength long side. An advantage with the method is that only one particular entangled state needs to be generated: all other necessary states can be produced from the first by means of, e.g., a small bank of phase plates with a prescribed birefringence. Unfortunately, it is not possible to generate larger patterns since the deposition methods proposed hitherto all are periodic with a period of half a wavelength Boto ; Kok ; Bjork . Restricting the exposure source of the lithography to four modes, with pairwise opposite wavevectors, Fig. 1, one can only increase the size of the pattern by some factor by sacrificing the pattern resolution by exactly the same factor. To be able to adjust the size of the deposited pattern independently of the resolution one must use more modes. In this paper we report a systematic multi-mode extension to the method we proposed earlier Bjork .

## II Sub-wavelength lithography

### II.1 One dimension, two modes

Our goal is to establish how to create arbitrary two-dimensional patterns on a squared substrate of side $`L`$. Suppose we have two counter-propagating beams in a direction we shall denote $`X`$, see Fig. 1. The beams propagate at angles $`\pm \theta `$ to the normal of the substrate. This substrate is coated with lithographic resist (in the following we will refer to the resist as the film) and situated in the region where the beams overlap. In general, provided the coherence lengths of the wave packets are much longer than the side of the film, we do not need to take into account the mode shapes and we will assume that they are plane over the side of the film.

Following Ref. Boto , the lithographic film absorption is modeled by an $`M`$-photon absorption process. Thus, the absorption process can be modeled by the operator $`\widehat{e}^{}{}_{}{}^{M}\widehat{e}_{}^{M}`$ given by

$$\widehat{e}^{}{}_{}{}^{M}\widehat{e}_{}^{M}\left(\frac{1}{\sqrt{W}}\underset{i}{}\widehat{a}_i^{}\right)^M\left(\frac{1}{\sqrt{W}}\underset{i}{}\widehat{a}_i\right)^M,$$
(1)
where $`\widehat{a}_i`$ is the annihilation operator of mode $`i`$, and $`W`$ is the number of excited modes impinging on the film. In this way, higher-order interference effects are naturally brought out. The photosensitive “grains” in the film must be much smaller than the shortest de Broglie wavelength encountered in the exposure process. Therefore, from the point of view of a “grain”, the photon packets in the respective modes will be indistinguishable in spite of their different linear momenta, as manifested by (1).

Let us now discuss the interference in one dimension between a pair of modes, labeled $`1,1`$, propagating in the plane defined by the $`X`$ axis and the film normal. We restrict ourselves to consider states that are eigenstates $`N_1`$ of the total photon number in the modes $`\pm 1`$. According to our assumptions the deposition rate in the substrate $`\mathrm{\Delta }_M`$ is proportional to the expectation value of the operator $`\widehat{e}^{}{}_{}{}^{M}\widehat{e}_{}^{M}`$, where $`\widehat{e}=(\widehat{a}_1+\widehat{a}_1)/\sqrt{2}`$. Let us further assume that the two beams impinge at the angles $`\theta _{\pm 1}=\pm \pi /2`$. The beams will hence strike the film surface at grazing incidence. Furthermore, we shall assume that the modes are prepared in a two-mode reciprocal binomial state of the general form

$$|\psi ^{(N_i)}=\frac{1}{\sqrt{𝒩_i}}\underset{n=0}{\overset{N_i}{}}\sqrt{n!(N_in)!}|n_i|N_in_i,$$
(2)
where $`N_i`$ is the total photon number of the two modes and $`𝒩_i=_{n=0}^{N_i}n!(N_in)!`$ is a normalization factor.

Let the $`X`$ coordinate normalized to the optical wavelength $`\lambda `$ be denoted $`x`$. Since the two modes $`1`$ and $`1`$ impinge over the film in anti-parallel directions, the accumulated phase of mode $`1`$ (propagating in the positive $`X`$ direction) at a distance $`\lambda x`$ from the left edge of the film will be

$$\widehat{U}_1=\mathrm{exp}(ik\lambda x\widehat{a}_1^{}\widehat{a}_1)=\mathrm{exp}(i2\pi x\widehat{a}_1^{}\widehat{a}_1),$$
(3)
where $`k=2\pi /\lambda `$, while mode $`1`$ will have accumulated the phase

$$\widehat{U}_1=\mathrm{exp}[ik\lambda (1x)\widehat{a}_1^{}\widehat{a}_1]=\mathrm{exp}[i2\pi (1x)\widehat{a}_1^{}\widehat{a}_1]$$
(4)
at the same location. Using these free-space unitary propagation operators, we find that at the location $`x`$, the state (2) for modes $`\pm 1`$ is transformed into

$$|\psi _x^{(N_1)}=\frac{1}{\sqrt{𝒩_1}}\underset{n=0}{\overset{N_1}{}}e^{i2\pi x(2nN_1)}\sqrt{n!(N_1n)!}|n_1|N_1n_1.$$
(5)
We can now translate the substrate a distance $`\lambda /[4(N_1+1)]`$ to the left, and at the same time phase-shift mode $`1`$ by $`2\pi \mathrm{}_{1x}/(N_1+1)`$ $`(\mathrm{}_{1x}=1,2,\mathrm{},N_1+1)`$ relative to mode $`1`$. The corresponding state will be labeled $`|\psi _x^{(N_1,\mathrm{}_{1x})}`$, where

$`|\psi _x^{(N_1,\mathrm{}_{1x})}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{𝒩_1}}}{\displaystyle \underset{n=0}{\overset{N_1}{}}}e^{i\pi (2x\frac{\mathrm{}_{1x}1/2}{N_1+1})(2nN_1)}`$ (6)
$`\times `$ $`\sqrt{n!(N_1n)!}|n_1|N_1n_1.`$
As shown in Bjork , this state will deposit a “one-dimensional pixel”, that is, the deposition rate $`\mathrm{\Delta }_M`$ will have a single pronounced peak $`\lambda /2(N_1+1)`$ wide, occupying the interval on the $`X`$-axis between $`\lambda (\mathrm{}_{1x}1)/2(N_1+1)`$ and $`\lambda \mathrm{}_{1x}/2(N_1+1)`$.

To make a qualitative comparison between the sub-wavelength resolution lithographic method proposed in Refs. Boto ; Kok  and our method, we have calculated the deposition pattern when the target pattern is a rectangular trench. In Refs. Boto ; Kok  such a trench, $`\lambda /4`$ wide, was used as a trial target function for a 10-photon state. We have done the same, but since a 10-photon state will define a 11-pixel pattern (in one dimension) the natural target trench function in our case is an integer number of pixel width wide. The pixel width for a 10-photon state is $`\lambda /22`$. In Fig. 2 we have calculated the deposition rate for a four pixel wide trench, that is, a trench $`2\lambda /110.18\lambda `$ wide. In order to make this pattern we can e.g. expose the pixels sequentially employing the states indicated in the figure caption. Although not shown in the figure, remember that this two-mode deposition rate is periodic with the period $`\lambda /2`$.

We see that the result of our method is almost the same, both in the respect of edge sharpness and in exposure penalty, to those obtained by the method proposed by Boto et al. Boto ; Kok . As we shall show below, neither the edge sharpness nor the exposure penalty need to be sacrificed when the lithographic pattern is extended over areas larger than a half a wavelength in each dimension. A fundamental difference between the methods is that the pattern producing state is pure in the proposal of Boto et al., while our proposal is based on mixed states (or a sequence of pure states if each pixel is deposited separately).

### II.2 One dimension, four modes

To overcome the limiting $`\lambda /2`$ periodicity of the deposition rate we introduce another pair of modes $`2`$ and $`2`$ impinging along the $`X`$ directions at angles $`\theta _{\pm 2}=\pm \mathrm{arcsin}[(N_2+1)^1]`$ from the film normal. Hence, they have only the wavevector components $`\pm 2\pi /[\lambda (N_2+1)]`$ in the $`X`$ direction. We assume that this pair of modes are prepared in a reciprocal binomial state, too. Consequently, their state at location $`x`$ is given by

$$|\varphi _x^{(N_2)}=\frac{1}{\sqrt{𝒩_2}}\underset{n=0}{\overset{N_2}{}}e^{i2\pi x\frac{2nN_2}{N_2+1}}\sqrt{n!(N_2n)!}|n_2|N_2n_2.$$
(7)
(We use the symbol $`\varphi `$ in the ket above to indicate that the modes corresponding to the state do not impinge at grazing incidence over the film surface.)

Now suppose the $`N_1+N_2`$ photon, product state $`|\psi _x^{(N_1)}|\varphi _x^{(N_2)}`$ is prepared. Calculating the pattern deposition rate $`\mathrm{\Delta }_M`$, where now $`\widehat{e}=(\widehat{a}_2+\widehat{a}_1+\widehat{a}_1+\widehat{a}_2)/2`$ and $`M=N_1+N_2`$, we find that

$`\mathrm{\Delta }_M`$ $``$ $`\left|{\displaystyle \underset{m=0}{\overset{N_1}{}}}e^{i2\pi x(2mN_1)}{\displaystyle \underset{n=0}{\overset{N_2}{}}}e^{i2\pi x\frac{2nN_2}{N_2+1}}\right|^2`$ (8)
$``$ $`{\displaystyle \frac{1}{[(N_1+1)(N_2+1)]^2}}{\displaystyle \frac{\mathrm{sin}^2[2(N_1+1)\pi x]}{\mathrm{sin}^2[2\pi x/(N_2+1)]}}.`$
The deposition rate $`\mathrm{\Delta }_M`$ has a highest oscillation period in $`x`$ of $`1/[2(N_1+1)]`$ and an overall periodicity of $`(N_2+1)/2`$, corresponding to the physical lengths $`\lambda /[2(N_1+1)]`$ and $`\lambda (N_2+1)/2`$, respectively. A plot of (8) for the case $`N_1=N_2=3`$ is shown in Fig. 3. Note that the deposition function spatial resolution is $`\lambda /8`$ and its periodicity is $`2\lambda `$.

When we translate the substrate a distance $`\lambda /[4(N_1+1)]`$ to the left, and at the same time phase-shift mode $`2`$ by $`2\pi \mathrm{}_{2x}/(N_2+1)+2\pi \mathrm{}_{1x}/[(N_1+1)(N_2+1)]`$ $`(\mathrm{}_{2x}=1,2,\mathrm{},N_2+1)`$ relative to mode $`2`$, the state $`|\varphi _x^{(N_2)}`$ in modes $`\pm 2`$ is transformed to

$`|\varphi _x^{(N_2,\mathrm{}_{1x},\mathrm{}_{2x})}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{𝒩_2}}}{\displaystyle \underset{n=0}{\overset{N_2}{}}}e^{i\pi (2x\mathrm{}_{2x}\frac{\mathrm{}_{1x}1/2}{N_1+1})\frac{2nN_2}{N_2+1}}`$ (9)
$`\times `$ $`\sqrt{n!(N_2n)!}|n_2|N_2n_2.`$
Using (6) and (9) we see that the four-mode state $`|\psi _x^{(N_1)}|\varphi _x^{(N_2)}`$ will consequently be transformed into the state $`|\psi _x^{(N_1,\mathrm{}_{1x})}|\varphi _x^{(N_2,\mathrm{}_{1x},\mathrm{}_{2x})}`$ after the translation and respective relative phase shifts. The deposition rate for this state can readily be calculated to be

$`\mathrm{\Delta }_M(\mathrm{}_{1x},\mathrm{}_{2x})`$ $``$ $`{\displaystyle \frac{1}{[(N_1+1)(N_2+1)]^2}}`$ (10)
$`\times `$ $`{\displaystyle \frac{\mathrm{sin}^2\{[2(N_1+1)x\mathrm{}_{1x}+1/2]\pi \}}{\mathrm{sin}^2\left[\left(2x\mathrm{}_{2x}\frac{\mathrm{}_{1x}1/2}{N_1+1}\right)\frac{\pi }{N_2+1}\right]}}.`$
If we divide the part of the $`X`$ axis between the origin and the point $`x=(N_2+1)/2`$ into $`(N_1+1)(N_2+1)`$ pieces, each $`1/[2(N_1+1)]`$ long, each interval will represent a “one-dimensional pixel”. With a specific choice of $`\mathrm{}_{1x}`$ and $`\mathrm{}_{2x}`$ we can deposit (or expose) pixel number $`\mathrm{}_{1x}+(N+1)\mathrm{}_{2x}`$ \[numbered from left to right and the number taken modulo $`(N_1+1)(N_2+1)`$\] with a negligible exposure penalty (that is, negligible unwanted exposure of nominally unexposed pixels). This can be clearly seen in Fig. 4, where we have assumed that $`N_1=N_2=3`$, $`\mathrm{}_{1x}=2`$, and $`\mathrm{}_{2x}=1`$, leading to the exposure of pixel number 6. The relative phase-shifts between the modes, labeled by the numbers $`\mathrm{}_{1x}`$ and $`\mathrm{}_{2x}`$, can be accomplished via a bank of appropriately chosen birefringent plates, provided that the modes $`\pm 1`$, and $`\pm 2`$, respectively, are originally prepared in spatially and temporally degenerate modes, but with orthogonal polarizations as discussed in Ref. Bjork .

The deposition-rate function (10) has two special properties that are worth pointing out. The first is that the deposition rate will be identically zero at the center of all unexposed pixels regardless of what other pixels are exposed. Hence, a nominally unexposed pixel surrounded by exposed pixels will remain unexposed at the pixel center. This is a very appealing feature of the proposed method since the exposure penalty hardly depends at all on the particular pixel pattern one intends to expose. The second nice feature is that the sum of the deposition-rate functions for all pixels add up to unity; i.e.

$$\underset{\mathrm{}_{1x}=1}{\overset{N_1+1}{}}\underset{\mathrm{}_{2x}=1}{\overset{N_2+1}{}}\mathrm{\Delta }_M(\mathrm{}_{1x},\mathrm{}_{2x})1,$$
(11)
for all values of $`x`$. This, in turn, means that we never risk overexposure, even if we expose two or more adjacent pixels. In fact, if a row, or column, of adjacent pixels are exposed, the resulting deposition function ridge will hardly have any modulation Bjork . The identity (11) also means that if one wants to make the negative image of some pixel pattern one can construct identically the negative image deposition rate by exposing all previously unexposed pixels, and vice versa.

Let us now discuss the geometrical scaling properties of the deposition rate. By decreasing the modes’ wavevector components in the film plane by a fixed factor, both the minimum feature size resolution and the fundamental period of deposition rate will increase by the same factor. If we, e.g., let modes $`\pm 1`$ impinge at angles $`\pm \mathrm{arcsin}(1/2)`$ from the film normal and modes $`\pm 2`$ impinge at the angles $`\pm \mathrm{arcsin}([2(N_2+1)]^1)`$, then the minimum feature size resolution (i.e., pixel size) becomes $`\lambda /(N_1+1)`$ and the period of deposition rate becomes $`(N_2+1)\lambda `$. However, the wavevector component parallel to the film is not only governed by the modes’ propagation angles but are also governed by the de Broglie wavelength of the impinging states. Therefore, the pixel and pattern sizes are intimately connected to how we prepare the states. If the (one-dimensional) film is modeled as a $`M=2N`$-photon absorber, the choice to partition the $`2N`$ photons equally between the two pairs of modes $`\pm 1`$, and $`\pm 2`$, as assumed in Fig. 4, is by no means necessary. Instead we can, e. g., prepare modes $`1`$ and $`1`$ in a two-mode $`(N1)`$-photon state $`|\psi _x^{(N1)}`$ and the modes $`2`$ and $`2`$ in a two-mode $`(N+1)`$-photon state $`|\varphi _x^{(N+1)}`$. The appropriate relative phase shifts are $`2\pi \mathrm{}_{1x}/N`$, where $`\mathrm{}_{1x}=1,2,\mathrm{},N`$ and $`2\pi \mathrm{}_{2x}/(N+2)+2\pi \mathrm{}_{1x}/(N^2+2N)`$ where $`\mathrm{}_{2x}=1,2,\mathrm{},N+2`$, respectively. In this case, the minimum lithographic feature size resolution becomes $`\lambda /(2N)`$, the number of individually depositable pixels become $`N(N+2)`$, and the fundamental period of the deposition rate becomes $`\lambda (N+2)/2`$. An illustration of an ensuing deposition rate function is given in Fig. 5. Continuing this re-partition one can either increase the fundamental period of the deposition rate at the expense of increasing the minimum resolution by increasing the photon number in modes $`\pm 2`$ at the expense of the photon number in modes $`\pm 1`$, or vice versa. The attainable minimum size resolution and deposition rate period are shown in Table 1.

### II.3 Two dimensions, eight modes

One can now extend the lithographic exposure procedure to two dimensions by simply introducing two additional pairs of modes $`\pm 1^{}`$ and $`\pm 2^{}`$, impinging towards the film at corresponding angles to modes $`\pm 1`$ and $`\pm 2`$, but in the $`Y`$ direction, perpendicular to $`X`$ Bjork . If the eight modes are prepared in the initial state $`|\psi _x^{(N,\mathrm{}_{1x})}|\varphi _x^{(N,\mathrm{}_{1x},\mathrm{}_{2x})}|\psi _y^{(N,\mathrm{}_{1y})}|\varphi _y^{(N,\mathrm{}_{1y},\mathrm{}_{2y})}`$, the deposition rate is given by the product of the corresponding deposition rates in the $`X`$ and in the $`Y`$ directions. Of course, if the number of photons in each of the two-mode states $`|\psi _x^{(N,\mathrm{}_{1x})}`$, …, $`|\varphi _y^{(N,\mathrm{}_{1y},\mathrm{}_{2y})}`$ is $`N`$, then the film must have a non-negligible $`M=4N`$-photon absorption cross section. If so, the assumed state will expose the pixel $`(\mathrm{}_{1x}+(N+1)\mathrm{}_{2x},\mathrm{}_{1y}+(N+1)\mathrm{}_{2y})`$ and leave the remaining pixels essentially unexposed. In order to expose a pattern, such as a line of adjacent pixels, one must prepare a statistical mixture of the pixels’ associated states or one can expose the pixels sequentially. As shown in Ref. Bjork , diagonal lines composed of exposed pixels will have an unacceptably large deposition rate fluctuation along the diagonal center-line. However, the minima can be “filled in” by depositing intermediate pixels (with their centers at the intersection points between four adjacent regular pixels). Again, the states corresponding to these intermediate pixels can be prepared by appropriate relative phase shifts between the mode pairs Bjork . For the rest of this paper we shall only study the deposition rates in one dimension, bearing in mind that with our method the two-dimensional deposition rate function is simply the product of two one-dimensional functions.

## III Generalized multi-mode quantum lithography

It is clear that the procedure to increase the deposition rate period is not limited by considering only two pairs of modes with opposite wavevectors in each dimension. We can continue this procedure by introducing a third pair of modes, labeled $`3`$ and $`3`$, impinging towards the film at the angles $`\theta _{\pm 3}=\pm \mathrm{arcsin}[1/(N+1)^2]`$ from the film normal. If the number of photons $`M`$ contributing to the film absorption process is divisible by 3 (or 6, in two dimensions), so that $`M=3N`$, and this photon number is partitioned equally between the three modes, then one will be able to deposit any pixellated patterns with the minimum feature size resolution of $`\lambda /[2(N+1)]`$ over a length of $`L=\lambda (N+1)^2/2`$.

However, in order to cover the maximum area for a given number of photons $`M`$ and resolvable feature size $`\lambda /[2(N+1)]`$, where $`1N<M`$, the following product state should be prepared:

$$|\psi _x^{(N,\mathrm{}_{1x})}|\varphi _x^{(1,\mathrm{}_{1x},\mathrm{}_{2x})}\mathrm{}|\varphi _x^{(1,\mathrm{}_{1x},\mathrm{}_{2x},\mathrm{},\mathrm{}_{(MN)x})},$$
(12)
where $`\mathrm{}_{1x}=1,2,\mathrm{},N+1`$, $`\mathrm{}_{2x},\mathrm{},\mathrm{}_{(MN)x}=1,2`$, and modes $`\pm 1`$ impinge at grazing incidence while modes $`\pm i`$, $`i=2,\mathrm{},MN`$, impinge at the angles $`\theta _{\pm i}=\pm \mathrm{arcsin}(2^{(i1)})`$. The rationale for preparing this state is that the state $`|\psi _x^{(N,\mathrm{}_{1x})}`$ will determine the feature size resolution and let us deposit any one of $`N+1`$ pixels each with a size of $`\lambda /[2(N+1)]`$. With the remaining $`MN`$ photons, each of the one photon states $`|\varphi _x^{(1,\mathrm{}_{1x},\mathrm{}_{2x})},\mathrm{},|\varphi _x^{(1,\mathrm{}_{1x},\mathrm{}_{2x},\mathrm{},\mathrm{}_{(MN)x})}`$ will allow us to double the fundamental period of the deposition rate function. If we compare this to a case where the $`MN`$ photons are partitioned between a smaller number of more highly excited states, it is clear that the state (12) gives a longer fundamental deposition rate period since $`2^{MN}MN+1`$ for all relevant numbers $`M`$ and $`N`$. With the initial state (12) one will be able to deposit any one of $`2^{MN}(N+1)`$ pixels in one dimension, where each pixel is $`\lambda /[2(N+1)]`$ wide. The fundamental period of the pixel pattern will be $`L=2^{MN1}\lambda `$. With $`2M`$ photons one will be able to make a two-dimensional pattern with this resolution and periodicity in both dimensions. This is the major result in this paper. In Fig. 6 we have plotted the one-dimensional deposition rate from the interference between eight modes in a statistical mixture of the states

$$|\psi _x^{(4,1)}|\varphi _x^{(1,1,1)}|\varphi _x^{(1,1,1,1)}|\varphi _x^{(1,1,1,1,0}$$
(13)
and

$$|\psi _x^{(4,3)}|\varphi _x^{(1,3,1)}|\varphi _x^{(1,3,1,1)}|\varphi _x^{(1,3,1,1,0},$$
(14)
where modes $`\pm 1`$ impinge at $`\theta _{\pm 1}=\pm \pi /2`$ and the remaining three pairs of modes impinge at the angles $`\theta _{\pm 2}=\pm \mathrm{arcsin}(1/2)`$, $`\theta _{\pm 3}=\mathrm{arcsin}(1/4)`$, and $`\theta _{\pm 4}=\pm \mathrm{arcsin}(1/8)`$. The total photon number in both states is six, the same number assumed in Fig. 4 and Fig. 5. The ensuing 32-pixel pattern is the largest one-dimensional pattern one can make with six photons provided that the minimum feature size is fixed to $`\lambda /8`$. The price for the large number of depositable pixels is the difficulty one will have generating this eight-mode state. If the states’ relative phase-shifts are generated by birefringent phase plates, one will need two plates for the first pair of modes, and three, four and five plates for the remaining three pairs of modes, respectively. We can in principle continue this process ad infinitum, but for a fixed minimum feature size, this requires $`M`$, the number of photons contributing to the exposure process, to increase. Since the absorption cross section quite generally decreases rapidly with increasing $`M`$, there will be a practical limit to such an extension.

## IV Imperfections due to losses and competing multi-photon absorption processes

Above we have discussed how sub-wavelength imaging can work under ideal conditions. However, in order for the proposed method to be of practical use it is necessary that it is robust against imperfections. Below we shall discuss three mechanisms that will deteriorate the image forming ability of entangled states: linear losses, competing multi-photon absorption processes, and exposure noise due to light quantization.

First we will discuss the effect of losses occuring between the state generator and the film. As long as the film strictly absorbs the same number of photons as the generated multi-mode state contains, losses will not affect the lithographic resolution, it will only lower the deposition rate by a fixed amount. This is rather obvious, because if one or more photons are lost from an $`M`$-photon state, no $`M`$-photon absorption process can be triggered by the state. However, if there exists a competing $`(M1)`$-photon absorption process in the film, the film may be exposed even after a photon is lost. In this case the modified deposition rate will be the same whether the photon is lost before impinging on the film or if only $`M1`$ out of $`M`$ photons impinging on the film are absorbed. Losses prior to the film will, however, shift the relative proportion between $`(M1)`$\- and $`M`$-photon absorption processes in favor of the former by decreasing the probability that the state impinging on the film contains $`M`$ photons. Therefore, losses in the optical system prior to the film should be kept as low as possible. Fortunately, losses will only gradually increase the required exposure dose and shift the probability of absorption toward absorption processes involving a smaller number of photons.

Next, we shall examine how the deposition rate of a $`M`$-photon state is affected by absorption processes involving less than $`M`$ photons. Two physical effects will deteriorate the resolution and the exposure penalty in this case. One is that if the state contains $`M`$ photons, but only $`M1`$ of them are absorbed, there are as many final states as there are modes. If we look at the simplest case, an impinging two-mode state $`|\psi _x^{(M)}`$, the possible final states are $`|1,0`$ and $`|0,1`$. Since these states are distinguishable, the absorption probability amplitudes leading to one of these final states will not interfere with the amplitudes leading to the other. Only the initial state $`|M,0`$ ($`|0,M`$) will evolve into the state $`|1,0`$ ($`|0,1`$) with certainty upon absorption of $`M1`$ photons. (All other number-difference states $`|Mn,n`$, $`n0,M`$ can evolve either into $`|1,0`$ or $`|0,1`$). This means that the two extreme number-difference states in the expansion of $`|\varphi _x^{(M)}`$ cannot interfere in a $`M1`$ photon absorption process. Therefore, the Fourier component with the highest spatial frequency will be absent in the ensuing deposition rate. Hence, the spatial resolution will decrease monotonically with decreasing order of the absorption process. A second effect will be that the destructive interference between the absorption probabilities outside the designated pixel will be incomplete, leading to an increased exposure penalty.

In Fig. 7 we have drawn the deposition rates of the state $`|\psi _x^{(4,3)}`$ due to 4-, 3-, 2-, and 1-photon absorption processes. The curves have been normalized such that the deposition rates all have a maximum of unity to facilitate comparisons. It can be seen that the width of the deposition rate peak increases, and so does the exposure penalty as the mismatch between the state and the absorption process increases. It is, of course, possible to make a better deposition rate function for, e. g., a 3-photon absorption film by exposing the film by a state of the kind $`|\psi _x^{(3)}`$ instead of the state $`|\psi _x^{(4,3)}`$. On the positive side, it is seen that if the film allows both for a 3-photon and a 4-photon absorption process, our method will yield a somewhat poorer result than if the 3-photon cross-section were identically zero. Still, the ensuing deposition rate is moderately deteriorated as compared to the ideal case.

In Fig. 8 we have drawn similar curves for a six-mode, four photon state that exposes pixel 7 (from the left) of 12 pixels. Each pixel is $`\lambda /6`$ wide, and the fundamental period of the deposition rate is $`2\lambda `$. Again, 3-, 2-, and 1-photon absorption processes will deteriorate the deposition rate, and in this case by a larger amount than for the two-mode state. The physical reason is that with a larger number of modes, there can also be a larger number of final states that will separate the interference paths into a larger number of distinguishable groups. This will primarily affect the (destructive) interference between the different absorption paths outside the deposition rate peak. However, in the six-mode case too, the deterioration is gradual. Hence, even if the 3-photon cross section is not identically zero, the consequences are not catastrophic.

The final effect we wish to discuss regards the fact that the calculated deposition rate is an ensemble average. In an experiment, the actual deposition rate may look rather different than its expectation value. This effect has not been discussed in any of the previous papers on entangled-state sub-wavelength lithography Boto ; Kok ; Bjork . If a state of the type given by (12) impinges on the film, the probability of state absorption is with all likelihood low. In addition, as have been discussed previously, each pixel defined on the film must contain many photo-sensitive grains since, for the method to make sense, the grains must be smaller (preferably much smaller) than the pixel size. Therefore, in order to expose a pixel, many exposure shots (per pixel) are needed. Since the absorption process is stochastic, the ensuing exposure will also be stochastic. The consequences of such an effect was studied in the context of the opposite process, namely image recognition, by Rose already in the 1940s Rose ; Rose 2  (an up-to-date review of this early work was recently written by Burgess Burgess ). In a first order approximation, the exposure of each pixel can be modeled by a Poisson distribution with a mean determined by the relation between the absorption probability of a grain, the deposition rate at the pixel, and the exposure dose. It is clear that in order to have pixels with a relative variation in exposure of, say less than 10%, the mean number of exposed grains must be larger than 100. This in turn implies that the number of states that need to impinge on the film to expose this particular pixel must be much larger than 100. As a consequence, it is desirable that the state generator emits states with a high repetition frequency in order to expose the film swiftly.

## V Conclusions

We have developed a multi-mode extension of the sub-wavelength lithographic method based on multi-photon absorption proposed in Ref. Bjork . We have shown the optimal way (in terms of the number of photons used in the process) to generate a pixellated pattern of given (sub-wavelength) resolution and pattern size. A salient feature of our proposal is that only one particular multi-mode photon number eigenstate needs to be prepared. Any pixellated pattern can be generated from this state by applying differential phase-shifts between the modes. For a given lithographic feature resolution, the price to be paid for a larger pattern is that the complexity of the state used to expose the pattern increases. For each doubling of the pattern area, two more modes in a one-photon state must be used. This also requires the film absorption process to increase by one in terms of photons absorbed in the process. Therefore it seems unlikely that it will be possible to make very large patterns. We have also studied the effects of imperfections, in terms of losses, competing lower order multi-photon processes, and deposition-rate fluctuations due to the quantization of the light used for the exposure. We have shown that neither effect results in catastrophic consequences, demonstrating that the proposed method is somewhat robust against imperfections.

###### Acknowledgements.

This work was supported by the Swedish Research Council for Engineering Sciences (TFR), the Swedish Foundation for Strategic Research (SSF) and L M Ericssons stiftelse för främjande av elektroteknisk forskning.