text:
# Spectra and eigenvectors of scale-free networks

## Abstract

We study the spectra and eigenvectors of the adjacency matrices of scale-free networks when bi-directional interaction is allowed, so that the adjacency matrix is real and symmetric. The spectral density shows an exponential decay around the center, followed by power-law long tails at both spectrum edges. The largest eigenvalue $`\lambda _1`$ depends on system size $`N`$ as $`\lambda _1N^{1/4}`$ for large $`N`$, and the corresponding eigenfunction is strongly localized at the hub, the vertex with largest degree. The component of the normalized eigenfunction at the hub is of order unity. We also find that the mass gap scales as $`N^{0.68}`$.

PACS numbers: 05.10.-a, 05-40.a, 05-50.+q, 87.18.Sn

Complex systems such as social, biological, and economic systems consist of many constituents such as individuals, substrates, and companies, respectively. Each constituent reacts and adapts to the pattern created in the system through diverse interactions. Cooperative phenomena between constituents in such complex systems may be described in terms of random graphs, consisting of vertices and edges, where vertices (edges) represent constituents (their interactions). The study of complex systems in terms of random graph was initiated by Erdös and Rényi (ER) . In the ER model, the number of vertices is fixed, and edges connecting from one vertex to another occur with probability $`p`$. Then there exists a probability threshold $`p_c`$, above which the system is percolated. Recently Watts and Strogatz (WS) introduced the small-world network  to consider local clustering, while the number of vertices is also fixed. The WS model offers the first indication that real networks could be more complex than predicted by the ER model.

Recently Barabási and Albert (BA)  introduced an evolving model where the number of vertices increases linearly with time rather than fixed, and a newly introduced vertex is connected to $`m`$ already existing vertices, following the so-called preferential attachment rule that the vertices with more edges are preferentially selected for the connection to the new vertex. The number of edges $`k`$ incident upon a vertex is the degree of the vertex. Then the degree distribution $`P(k)`$ of vertices, equivalent to the connectivity distribution, follows a power-law $`P(k)k^3`$ for the BA model, while for the ER and WS models, it follows a Poisson distribution. The BA model is interesting in a sense that a lot of complex interactive networks such as the world wide web , the actor network , the citation network of scientific papers and the author collaboration network, metabolic network  and food web  in biological systems all exhibit the power-law in the degree distribution, implying that a characteristic degree is absent in such systems. Thus the BA model is called a scale-free network. In the BA network, one may assign a direction on each edge pointing from the younger vertex to the older one. However, when the direction on each edge is ignored, allowing bi-directional interactions such as download and upload communications in Internet, the BA network may be regarded as a simple model for the Internet topology.

While it is known that the BA network follows the power-law in its degree distribution, further structural properties are not so well known. When $`m=1`$, the BA network forms a tree structure without forming any loop, but for $`m>1`$, loops are formed, and network topology becomes much complicated. So it would be interesting to investigate the spectrum of the BA network, because generally the spectrum of a random graph and corresponding eigenvectors are closely related to topological features of the random graph . In this paper, we study the spectrum and the corresponding eigenvectors of the adjacency matrix of the BA network, comparing spectral properties with structural features. When a BA network is composed of $`N`$ vertices, the adjacency matrix $`𝐀`$ consists of $`N\times N`$ elements, $`\{a_{i,j}\}`$ ($`i,j=1,\mathrm{},N`$), defined as $`a_{i,j}=1`$ if vertices $`i`$ and $`j`$ are connected via an edge, and $`a_{i,j}=0`$ otherwise. $`k_i=_ja_{i,j}`$ is the degree of the vertex $`i`$. In the BA model, the vertex with the largest degree is singled out and is called the hub, denoted by $`h`$ hereafter, in this work. From previous studies , we know that $`k_h`$, the degree of the hub, scales as $`k_hmN^{1/2}`$. The vertices are ordered in their ages; the vertex $`i`$ is the one born at time $`i`$.

Since we consider the bi-directional case (sometimes called the undirected case), the adjacency matrix $`𝐀`$ is real and symmetric, so that all eigenvalues are real and the largest eigenvalue is not degenerate. We obtained the spectrum of the BA network via exact diagonalization for $`N`$ up to $`5,000`$ and for the first few largest eigenvalues, via the Lanczos method  for $`N`$ up to $`400,000`$. Throughout this paper, numerical simulations were carried out for $`m=2`$, the simplest case including the loop structure.

$`Eigenvalues`$ We consider the distribution of the eigenvalues, shown in Fig.1. Analytic formula for the spectrum is not known yet, but one can easily see that the spectrum does not fit to the semi-circular equation derived by Wigner appropriate to the ER random graph. The density of eigenvalues $`\rho (\lambda )`$ in the middle part of the spectrum is likely to fit to the formula, $`\rho (\lambda )\mathrm{exp}(|\lambda |/a)`$ where $`a1.25`$ (see the left inset of Fig.1), while the density further out follows the power-law, $`\rho (\lambda )|\lambda |^4`$, (see the right inset of Fig.1). Since the power-law decays much slower than the exponential one, the spectrum shows long tails at both spectrum edges (see Fig.1). Such behavior has also been observed in the financial time series .

The size dependence of the largest eigenvalue $`\lambda _1`$ is of interest and we found numerically that $`\lambda _1`$ increases as $`N^{1/4}`$ for large systems (see Fig.2). To consider the relation between the largest eigenvalue $`\lambda _1`$ and the structure of the BA network, we consider the followings. Let $`𝐱`$ and $`𝐲`$ be $`N\times 1`$ column vectors, related through

$$𝐲(n)=𝐀^n𝐱.$$
(1)
Let $`v_{i,l}`$ be the $`i`$-th component of the normalized eigenvector corresponding to the $`l`$-th eigenvalue of $`𝐀`$, $`\lambda _l`$, ($`\lambda _1>\lambda _2\mathrm{}\lambda _N`$). Then Eq.(1) may be rewritten as

$$y_i(n)=\underset{l}{}\underset{j}{}\lambda _l^nv_{i,l}v_{j,l}x_j,$$
(2)
When we set $`x_j=1`$ for some $`j`$ and $`x_r=0`$ for $`rj`$, $`y_i(n)`$ becomes the number of possible ways of $`n`$-step walks starting from the vertex $`j`$ and terminating at the vertex $`i`$. This will be denoted as $`y_{ji}(n)`$ hereafter. In particular, when $`i=j`$, $`y_{jj}(n)`$ is the number of possible ways to return to the starting vertex $`j`$ after $`n`$-step walks. When $`n`$ is sufficiently large, $`y_{ji}(n)`$ can be represented in terms of the largest eigenvalue $`\lambda _1`$ alone as

$$y_{ji}(n)\lambda _1^nv_{i,1}v_{j,1}.$$
(3)
In particular, when $`i=j`$,

$$y_{jj}(n)\lambda _1^nv_{j,1}^2.$$
(4)
Even though the relation Eq.(4) holds for any vertex $`j`$, we consider it particularly at the hub. Now let us consider the number of ways of $`n+2`$ steps returning to the starting vertex $`h`$. It can be split into two parts,

$$y_{hh}(n+2)=y_{hh}(n)y_{hh}(2)+\underset{jh}{}y_{hj}(n)y_{jh}(2).$$
(5)
According to Eq.(4),

$$\lambda _1^2\frac{y_{hh}(n+2)}{y_{hh}(n)}=y_{hh}(2)+\frac{\underset{jh}{}y_{hj}(n)y_{jh}(2)}{y_{hh}(n)},$$
(6)
where $`y_{hh}(2)`$ corresponds to the degree of the hub, $`y_{hh}(2)=m\sqrt{N}`$, while $`y_{jh}(2)`$ is found numerically to behave as $`N^{0.05}/j^{0.43}`$, weakly depending on $`N`$ (see Fig.3). The second term of the right hand side of Eq.(6) can be written using Eqs.(3) and (4) as $`_{jh}\left(\frac{v_{j,1}}{v_{h,1}}\right)y_{jh}(2)`$. We will show later that $`v_{h,1}1/2`$, while $`v_{j,1}(Nj)^{1/4}`$ for $`j`$ being away from $`h`$, leading to the result that the second term becomes $`𝒪(N^{0.1})`$. Consequently, $`\lambda _1`$ is contributed dominantly by the first term in the right hand side of Eq.(6), leading to $`m^{1/2}N^{1/4}`$.

The result of $`\lambda _1N^{1/4}`$ can also be understood through two toy models. First, one vertex is located at the center, and all the other $`N1`$ vertices are linked only to the center vertex (Fig.4a). This structure, called the radial structure, is an extreme case of “winner-takes-all”. The largest eigenvalue of this structure is $`\lambda _1=\sqrt{N1}`$, while the largest degree is $`N1`$. Second, we consider a two-level Cayley tree structure (Fig.4b). When the coordination number is chosen as $`k=\sqrt{N}`$, the total number of vertices becomes $`N+1`$, and the largest degree $`\sqrt{N}`$, as in the case of the BA network. The largest eigenvalue for this toy model is found to be $`\lambda _1=\sqrt{2k1}N^{1/4}`$. So, in both models, $`\lambda _1\sqrt{k_h}`$.

The second largest eigenvalue in its absolute magnitude is located on the negative side of the spectrum. The absolute magnitude of this eigenvalue $`|\lambda _N|`$ is also found to scale as $`|\lambda _N|N^{1/4}`$ (see Fig.2). The difference between the first two largest eigenvalues in their absolute magnitude, $`\lambda _1|\lambda _N|`$, is found to scale as $`N^{0.43}`$ (see the inset of Fig.2). Thus, the mass gap, defined as $`\mathrm{log}(\lambda _1/|\lambda _N|)`$ , scales as $`N^z`$ with $`z0.68`$.

Next, we consider the size dependence of the moments of the spectrum. First, since the matrix $`𝐀`$ is traceless, the first moment becomes zero, that is, $`_1=_i^N\lambda _i=0`$. Second, for each vertex, the number of ways to return to the starting vertex by a 2-step walk is the same as the degree of that vertex, so that $`_2=_i^N\lambda _i^2=2mN`$. Third, since the path to return to the starting vertex by a 3-step walk forms a triangle, the total number of triangles $`𝒯_N`$ in the system can be obtained through the third moment, $`_3=_i^N\lambda _i^3=6𝒯_N`$. We found numerically that $`_3N^{0.40}`$ (see Fig.5). Note that while the first moment vanishes, the third moment does not, implying that the spectrum is not completely symmetric with respect to $`\lambda =0`$.

$`Eigenvector`$ Since the adjacency matrix $`𝐀`$ is real and symmetric, every component of the eigenvector corresponding to the largest eigenvalue is positive. Since the largest eigenvalue and the corresponding eigenfunction is important, as in the case of Eq.(3), we focus on the eigenvector $`\{v_{j,1}\}`$ for the largest eigenvalue. We consider the square of each component of the eigenvector $`\{v_{j,1}^2\}`$ instead of $`\{v_{j,1}\}`$ itself because of the normalization, $`_jv_{j,1}^2=1`$, and compare it with the normalized degree of each vertex $`\{k_j/_lk_l\}`$. As seen in Fig.6, the two quantities are correlated in such a way that up-down behavior occurs in a similar fashion.

The components $`\{v_{j,1}^2\}`$ ($`j=1,\mathrm{},N`$) is strongly localized at the hub. This is in contrast to the ER case where the corresponding eigenfunction is extended over all vertices (see the inset of Fig.6). We found that $`v_{h,1}^2`$ at the hub increases with increasing $`N`$ for small $`N`$, but converges to a constant close to $`1/2`$ for large $`N`$ (see Fig.7). The value $`1/2`$ can be obtained from the radial structure analytically. Since the radial structure is the extreme case that the connectivity is mostly localized at the hub, the value $`1/2`$ can become the upper bound of $`v_{h,1}^2`$. The result that $`v_{h,1}^2`$ approaches a constant as $`N\mathrm{}`$ is interesting, because $`v_{h,1}^2`$ behaves differently from the normalized degree at the hub, decreasing as $`N^{1/2}`$ (see Fig.7). Since the number of returning $`n`$-step walks from a vertex $`j`$ is $`y_{jj}(n)v_{j,1}^2`$ from Eq.(4), $`v_{j,1}^2`$ can be considered as a measure of contribution from the vertex $`j`$ to transport processes in a network. The fact that $`v_{h,1}^2𝒪(1)`$ can be interpreted as that the hub plays much more dominant role in transports on the BA network, compared with the contribution measured by the normalized degree at the hub, which is $`𝒪(1/\sqrt{N})`$. Thus, the presence of the hub in the scale-free networks diversifies pathways, and enhances the efficiency of transport, while the network systems are vulnerable to the attack on the hub . On the other hand, for other components of the eigenfunction, $`\{v_{j,1}^2\}`$, ($`j=1,\mathrm{},N`$, but $`jh`$), it was found numerically that $`v_{j,1}^21/4\sqrt{Nj}`$ (see Fig.8).

Each component of eigenvectors can be obtained through the diagonalization of the adjacency matrix $`𝐀`$. However, the eigenfunction for $`\lambda _1`$ can be obtained using the relation, Eq.(4). The knowledge of $`y_{jj}(n)`$ and $`\lambda _1`$ enables one to obtain each component of the eigenfunction through $`v_{j,1}^2y_{jj}(n)/\lambda _1^n`$ for sufficiently large $`n`$. However, when $`n`$ is not large enough, $`y_{jj}(n)`$ for vertex $`j`$ is affected by local topology around $`j`$, so that the value of $`v_{j,1}^2`$ obtained by this method would be different from that by the exact diagonalization. Thus, we can define a characteristic number of steps $`n_c`$, such that for $`n>n_c`$, Eq.(4) holds, while for $`n<n_c`$, it breaks down. To find $`n_c`$, we introduce the quantity,

$$\delta _n\left|\underset{j}{}\left(\frac{y_{jj}(n)}{\lambda _1^n}v_{j,1}^2\right)\right|,$$
(7)
which is found to decay exponentially as $`\delta _n\mathrm{exp}(n/n_c)`$ (see the inset of Fig.9). We found numerically that $`n_c`$ shows a crossover across a characteristic system size $`N_c`$ such that for $`N<N_c`$, $`n_cN^{0.35}`$, while for $`N>N_c`$, $`n_cN^{0.50}`$ (see Fig.9). On the other hand, Eq.(7) may be rewritten as $`\delta _n=|_{l>1}(\lambda _l/\lambda _1)^n|.`$ Then, for sufficiently large $`n`$, $`\delta _n`$ is contributed dominantly by the term $`(|\lambda _N|/\lambda _1)^n`$ alone. Combined with $`\delta _n\mathrm{exp}(n/n_c)`$, one obtains that $`n_c=\lambda _1/(\lambda _1|\lambda _N|)`$, equivalent to the inverse of the mass gap. Thus, $`n_cN^z`$. However, the numerical value of the exponent $`z`$ obtained in this way, $`z0.50`$, deviates from the value $`z0.68`$ obtained by the direct measurement of $`\lambda _1|\lambda _N|`$. This can be attributed to the fact that $`\delta _n`$ includes contributions from other eigenvalues.

$`Shortest`$ $`paths`$ Since the role of the hub in transport is much more dominant, compared with the contribution by the normalized degree, we also study topological feature of shortest paths between two vertices . The transport from one position to another is mainly carried along the shortest path(s) between them, and is contributed dominantly by the largest eigenvalue of the adjacency matrix $`𝐀`$. We define a set, composed of the vertices on the shortest path(s) from one position to another. Then, since there are $`N(N1)/2`$ pairs of vertices, the same number of sets exist in the system. Among them, we are interested in how many number of different sets a certain vertex $`j`$ belongs to. This number is called the involving number $`P_j`$, while the normalized involving number is defined as $`p_jP_j/_lP_l`$ for each vertex $`j`$. Fig.6 also shows $`p_j`$ versus $`j`$. It behaves similarly to the up-down behavior of $`v_{j,1}^2`$ and $`k_j`$. In particular, the involving number at the hub $`P_h`$ is found numerically to scale as $`P_hN^2`$, while the total involving number summed over all vertices is found numerically to scale as $`_{j=1}^NP_jN^2\mathrm{log}N`$ (see Fig.10) . So the normalized involving number at the hub $`p_hP_h/_jP_j`$ behaves as $`p_h1/\mathrm{log}N`$, decaying much more slowly, compared with the normalized degree decreasing as $`N^{0.5}`$. This weak dependence on $`N`$ is comparable to the result that $`v_{h,1}𝒪(1)`$ for large $`N`$. Therefore, the contribution of the hub to shortest paths is much larger than that of the naive estimate based on the normalized degree at the hub, $`𝒪(1/\sqrt{N})`$.

In conclusion, we have considered the spectrum and eigenvectors of the adjacency matrix of the BA network, when bi-directional interaction is allowed. The density of the eigenvalues decays exponentially for small $`|\lambda |`$s, followed by power-law tails at both spectrum edges. This is different from the Wigner’s formula appropriate to random graphs. We found that the largest two eigenvalues $`\lambda _1`$ and $`|\lambda _N|`$ depend on system size $`N`$ as $`N^{1/4}`$ for large $`N`$, and the mass gap scales as $`N^{0.68}`$. The eigenfunction corresponding to the largest eigenvalue $`\lambda _1`$ is strongly localized at the vertex with the largest degree, called the hub. The component of the normalized eigenfunction for $`\lambda _1`$ at the hub is independent of $`N`$, implying that the role of the hub in transport on the scale-free network becomes much important for larger system, and its contribution becomes much more dominant than expected according to the normalized degree at the hub, which scales as $`N^{1/2}`$. Therefore, it is very efficient in communication networks to construct central vertices, through which most of information traffic pass.

We would like to thank S.Y. Park for providing us with the Lanczos algorithm codes. This work is supported by grants No.2000-2-11200-002-3 from the BRP program of the KOSEF. While writing this manuscript we learned of a work by Farkas $`et`$ $`al.`$ , which overlaps some of our results.


format: text

text:
# Untitled Document

NEW IMPROVEMENTS FOR MIE SCATTERING CALCULATIONS

V. E. Cachorro

Departamento de Física Aplicada I

Valladolid University, 47071 Valladolid, SPAIN

L. L. Salcedo

Departamento de Física Moderna

Granada University, 18071 Granada, SPAIN

ABSTRACT

New improvements to compute Mie scattering quantities are presented. They are based on a detailed analysis of the various sources of error in Mie computations and on mathematical justifications. The algorithm developed on these improvements proves to be reliable and efficient, without size ($`x=2\pi R/\lambda `$) nor refractive index ($`m=m_R\mathrm{i}m_I`$) limitations, and the user has a choice to fix in advance the desired precision in the results. It also includes a new and efficient method to initiate the downward recurrences of Bessel functions.

1. INTRODUCTION

The Mie theory of light scattering by a homogeneous sphere is used for many problems of atmospheric optics and also in other fields in Physics. The application of Mie theory still needs modern computers for numerical calculations of the many functions and coefficients involved. The primary difficulty is in the precise evaluation of expansion coefficients $`a_n`$ and $`b_n`$. This is further aggravated as $`x`$ gets large, and when the calculation of size distribution is needed. An optimization of computer time for reliable computation is clearly of necessity.

The formulas for Mie scattering are well known1,2. Here we follow the notation of Bohren and Huffman3. The scattering and extinction efficiency factors are given by

$$\begin{array}{cc}\hfill Q_s& =\frac{2}{x^2}\underset{n=1}{\overset{N}{}}(2n+1)\left(|a_n|^2+|b_n|^2\right)\hfill \\ \hfill Q_e& =\frac{2}{x^2}\underset{n=1}{\overset{N}{}}(2n+1)\mathrm{Re}(a_n+b_n)\hfill \end{array}$$
$`(1)`$
where $`x=2\pi R/\lambda `$ is the size parameter of the problem, $`R`$ being the radius of the sphere, $`\lambda `$ the wavelength of the light and $`N`$ a large enough number. The Mie scattering coefficients $`a_n`$ and $`b_n`$ are functions of $`x`$ and the relative refractive index $`m=m_R\mathrm{i}m_I`$, with $`m_R1,m_I0`$.

$$\begin{array}{cc}\hfill a_n& =\frac{x\psi _n(x)\psi _n^{}(y)y\psi _n^{}(x)\psi _n(y)}{x\zeta _n(x)\psi _n^{}(y)y\zeta _n^{}(x)\psi _n(y)}\hfill \\ \hfill b_n& =\frac{y\psi _n(x)\psi _n^{}(y)x\psi _n^{}(x)\psi _n(y)}{y\zeta _n(x)\psi _n^{}(y)x\zeta _n^{}(x)\psi _n(y)}\hfill \end{array}$$
$`(2)`$
where $`y=mx`$ and $`\psi _n(z),\zeta _n(z)`$ are the Riccati-Bessel functions related to the spherical Bessel functions $`j_n(z)`$ and $`y_n(z)`$:

$$\begin{array}{cc}\hfill \psi _n(z)& =zj_n(z)\hfill \\ \hfill \zeta _n(z)& =zj_n(z)\mathrm{i}zy_n(z)\hfill \end{array}$$
$`(3)`$
These functions are known in closed form (Ref. 4, p. 437) but it is more convenient to use the recurrence relation

$$\begin{array}{cc}\hfill X_{n+1}(z)& =F_n(z)X_n(z)X_{n1}(z),\hfill \\ \hfill F_n(z)& =(2n+1)/z.\hfill \end{array}$$
$`(4)`$
where $`X`$ is any of the functions in eqn. (3).

Presently, there are many versions of Mie scattering computer codes (Dave5,6, Blattner7, Grehan and Gouesbet8,9, Wiscombe10,11, Goedecke et al.12, Miller13) and authors who had been doing Mie calculations (Kattawar and Plass14, Deirmendjian15, Quenzel and Müller16, Bohren and Huffman3). These are reflected in performing our work.

Some essential points should be addressed by any Mie scattering algorithm:

1) How to determine the number $`N`$ for truncating a Mie series.

2) Whether the Riccati-Bessel functions will be computed by upward recursion or by downward recursion.

3) If downward recursion is used, how to initialize it.

4) How to structure the algorithm in an efficient way.

Answers to all the above questions constitute the objective of this paper. We focus particularly on analyzing the numerical error sources and show that our Mie algorithm permits users to prescribe a precision $`ϵ`$ beforehand, to effect an efficient, reliable Mie coefficients calculation. Needless to say, the precisely evaluated Mie coefficients $`a_n,b_n`$ are required for calculating the angular scattering amplitudes1,2,3,5,6,10.

2. CONVERGENCE PROPERTIES OF THE MIE SERIES

In this section we shall estimate the error introduced in some typical quantity such as the efficiency factors, by keeping a finite number $`N`$ of partial waves in the Mie series. We shall also find a criterion for choosing the value of $`N`$. In this section the quantities $`a_n,b_n`$ themselves are assumed to be computed exactly.

In order to investigate the convergence properties of the scattering coefficients $`a_n,b_n`$ we shall make use of very well known properties of the spherical Bessel functions (e.g. ref. 4, p. 438 and ff.). Let us recall some properties which are relevant for us:

i)

$$\underset{n\mathrm{}}{lim}\psi _n(z)=0,\underset{n\mathrm{}}{lim}\zeta _n(z)=\mathrm{}.$$
$`(5)`$
ii) For $`z=x`$ real, $`\psi _n(x)`$ and $`\zeta _n(x)`$ have two distinct regimes as functions of $`n`$:

a) oscillating regime for $`n<x`$. $`\psi _n(x)`$ and $`\zeta _n(x)`$ keep changing their sign regularly, and $`|\psi _n(x)|`$ and $`|\zeta _n(x)|`$ are bounded by slowly changing functions of $`n`$.

b) exponential regime for $`n>x`$. $`\psi _n(x)`$ becomes exponentially decreasing and $`|\zeta _n(x)|`$ becomes exponentially increasing.

In view of these considerations one concludes from eqn. (2), that all the partial waves $`n<x`$ ($`x`$ being the size parameter from now on) will contribute to the Mie series and convergence will appear only after $`n`$ enters in the exponential regime. This is so because $`\psi _n(x),\psi _n^{}(x)`$ go very quickly to zero in the numerator and $`\zeta _n(x),\zeta _n^{}(x)`$ go to infinity in the denominator. On the other hand $`\psi _n(y),\psi _n^{}(y)`$ appear both in numerator and denominator and therefore seem to play no role in the convergence. We can emphasize this fact by writing

$$\begin{array}{cc}\hfill a_n& =\frac{\psi _n(x)}{\zeta _n(x)}[a]_n=\frac{\psi _n(x)}{\zeta _n(x)}\frac{n(y/xx/y)+xA_n(y)yA_n(x)}{n(y/xx/y)+xA_n(y)yB_n(x)}\hfill \\ \hfill b_n& =\frac{\psi _n(x)}{\zeta _n(x)}[b]_n=\frac{\psi _n(x)}{\zeta _n(x)}\frac{yA_n(y)xA_n(x)}{yA_n(y)xB_n(x)}\hfill \end{array}$$
$`(6)`$
where we have extracted the factor $`\psi _n(x)/\zeta _n(x)`$ responsible for the convergence of $`a_n`$ and $`b_n`$ and also we have reexpressed the ratios $`\psi _n^{}(z)/\psi _n(z)`$ and $`\zeta _n^{}(x)/\zeta _n(x)`$ in terms of (ref. 4, p. 439)

$$A_n(z)=\frac{\psi _{n1}(z)}{\psi _n(z)},B_n(x)=\frac{\zeta _{n1}(x)}{\zeta _n(x)}$$
$`(7)`$
Let us state more clearly our assumption: we shall assume that the quantities $`[a]_n,[b]_n`$ are bounded by slowly varying functions of $`n`$ in the exponential regime $`n>x`$. The validity of this assumption will be analyzed in a later section.

If $`[a]_n`$ and $`[b]_n`$ are well behaved for large $`n`$, we can approximate them by their asymptotic values in order to discuss the convergence of $`a_n`$ and $`b_n`$. In order to take advantage of this approximation we can use the asymptotic expansion of the Bessel functions for large orders (ref. 4, p. 365),

$$A_n(z)F_n(z),B_n(x)F_n^1(x)\mathrm{as}n\mathrm{}$$
$`(8)`$
where the next term in the expansion has a higher power of $`1/n`$. We obtain

$$[a]_n\frac{1m^2}{1+m^2}+O\left(\frac{1}{n}\right),[b]_nO\left(\frac{1}{n}\right)$$
$`(9)`$
In practice, for $`xnN`$$`[a]_n`$ and $`[b]_n`$ are both of the order of unity, (unless $`m`$ is nearly 1, in which case $`[a]_n,[b]_n0`$). On the other hand, recalling that $`m_R1`$, it can be proved that $`|\frac{1m^2}{1+m^2}|<2`$, therefore a good enough estimate is

$$[a]_n,[b]_n1$$
$`(10)`$
Using this and the asymptotic values (8), it can be shown that the truncation error in $`Q_e`$ is bounded by

$$\delta Q_e|a_N|$$
$`(11)`$
The proof is presented in Appendix I where it is shown that the series $`_{n=N+1}^{\mathrm{}}|a_n|`$ converges faster than some geometric series. Let us note that what actually appears in $`Q_e`$ is Re$`a_n`$, not $`|a_n|`$, therefore the bound (11) will usually be conservative. This is especially true for small $`m_I`$ because in this case Re$`a_n|a_n|^2`$ (i.e. $`Q_eQ_s`$) and $`|a_n|^2|a_n|`$ for $`n>N`$.

Let us now find a criterion for choosing the number $`N`$ of partial waves that should be taken into account. For this purpose let $`ϵ`$ be the error allowed in the calculation, and let us take $`\delta Q_e`$ as a typical quantity in the problem. Then $`N`$ should be taken so that

$$\delta Q_eϵ$$
$`(12)`$
Taking the quantity $`Q_e`$ has the advantage of being simple and also that $`\delta Q_s\delta Q_e`$, because $`|a_n|^2<\mathrm{Re}a_n`$ (i.e. $`Q_sQ_e`$ for each partial wave). Other interesting quantities, such as the scattering amplitudes, have similar convergence properties as $`Q_e`$ and $`Q_s`$.

Putting together the bound (11), the criterion (12) and the estimate (10) we find the following prescription

$$\left|\frac{\psi _N(x)}{\zeta _N(x)}\right|ϵ$$
$`(13)`$
In order to find something more convenient let us make use of the Wronskian identity (ref. 4, p.439)

$$\psi _n(x)\zeta _{n1}(x)\psi _{n1}(x)\zeta _n(x)=\mathrm{i},$$
$`(14)`$
and the asymptotic values of $`A_n(x)`$ and $`B_n(x)`$. In this way we obtain (within approximations keeping the order of magnitude)

$$\psi _n(x)\zeta _n(x)\mathrm{i}F_n^1(x),$$
$`(15)`$
This allows us to remove $`\psi _n(x)`$ from (13) and finally we obtain the prescription for $`N`$

$$|\mathrm{Im}\zeta _N(x)|\sqrt{\frac{1}{ϵ}},$$
$`(16)`$
which has been written in a form convenient for being checked while $`\zeta _n(x)`$ is being computed by upward recurrence. In getting (16) we have neglected a factor $`F_N(x)`$ from (15) because by doing so $`N`$ may increase at most by one unit (recall that $`\zeta _n(x)/\zeta _{n1}(x)F_n(x)`$). Also we have used that Re$`\zeta _n(x)=\psi _n(x)`$ is negligible as compared to Im$`\zeta _n(x)`$ in the exponential region.

It is remarkable that the value of $`N`$ obtained from (16) for $`ϵ=10^8`$ is virtually identical to the standard prescription $`N=x+cx^{1/3}+1`$, with $`c=4.3`$. It is shown in Appendix II that it must be so using asymptotic expansions for $`\zeta _n(x)`$, and also how to modify $`c`$ if some other precision $`ϵ`$ is desired. To know $`N(x)`$ in advance is necessary if the computer code is to be vectorized10,11.

3. NUMERICAL ERROR AND UPWARD RECURRENCE

In this section we shall discuss the propagation of numerical error through the calculation.

It is known that the determination of $`\psi _n(z)`$ by upward recursion is intrinsically unstable (see e.g. ref. 5). Let us clarify this point.**We thank one of the referees for providing us with a simpler proof of this statement. For the sake of simplicity let us assume that the numerical error is coming from the initial values

$$\stackrel{~}{\psi }_0(z)=\psi _0(z)+ϵ_0,\stackrel{~}{\psi }_1(z)=\psi _1(z)+ϵ_1$$
$`(17)`$
but the recursion itself is free of roundoff error, i.e.

$$\stackrel{~}{\psi }_{n+1}(z)=F_n(z)\stackrel{~}{\psi }_n(z)\stackrel{~}{\psi }_{n1}(z)$$
$`(18)`$
$`ϵ_0,ϵ_1`$ being small numbers depending on the precision of the computer, and $`\stackrel{~}{\psi }_n(z)`$ being the numerical sequence that is actually obtained instead of the exact one, $`\psi _n(z)`$. Subtracting the exact recursion for $`\psi _n(z)`$ from (18) we find

$$\delta \psi _{n+1}(z)=F_n(z)\delta \psi _n(z)\delta \psi _{n1}(z)$$
$`(19)`$
where $`\delta \psi _n(z)=\stackrel{~}{\psi }_n(z)\psi _n(z)`$ is the error in our numerical sequence. Any sequence satisfying the recurrence relation (4) is a linear combination of $`\psi _n(z)`$ and $`\zeta _n(z)`$, therefore

$$\delta \psi _n(z)=\eta \psi _n(z)+\eta ^{}\zeta _n(z)$$
$`(20)`$
The small numbers $`\eta ,\eta ^{}`$ are directly related to $`ϵ_0,ϵ_1`$ through eqn. (17), namely

$$\begin{array}{cc}\hfill \eta & =\mathrm{i}(ϵ_0\zeta _1ϵ_1\zeta _0)\hfill \\ \hfill \eta ^{}& =\mathrm{i}(ϵ_0\psi _1ϵ_1\psi _0)\hfill \end{array}$$
$`(21)`$
Recalling now that $`\zeta _n(z)`$ diverges for large $`n`$ we conclude that the absolute error in $`\stackrel{~}{\psi }_n(z)`$ will eventually blow up. More generally, if the recursion itself is not exact due to computer roundoff error, $`\stackrel{~}{\psi }_n(z)`$ is rather given by

$$\delta \psi _n(z)=\eta _n\psi _n(z)+\eta _n^{}\zeta _n(z)$$
$`(22)`$
where $`\eta _n,\eta _n^{}`$ are of the order of the roundoff error or the initial values error, whichever the largest. In any case the conclusion is still that $`\delta \psi _n(z)`$ is small for small $`n`$ (or while $`n`$ is in the oscillating regime for $`z`$ nearly real), but blows up when $`n`$ enters in the exponentially increasing regime of $`\zeta _n(z)`$. Since $`\psi _n(z)`$ itself goes to zero in the exponential regime, $`\stackrel{~}{\psi }_n(z)`$ has less and less correct figures at each step.

We can extract some corollaries from the previous discussion:

1) The upward recursion is always unstable for computing $`\psi _n(z)`$ for large $`n`$, depending on $`z`$. The error $`\delta \psi _n(z)`$ grows as $`|\zeta _n(z)|`$. On the other hand the upward recursion is perfectly stable for computing $`\zeta _n(x)`$ for any value of $`n`$. This is because $`\delta \zeta _n(x)`$ still grows as $`|\zeta _n(x)|`$, therefore the relative error in $`\zeta _n(x)`$ is kept small. Note however that the relative error in the quantity Re$`\zeta _n(x)=\psi _n(x)`$ is not at all small.

2) A downward recursion is stable for computing $`\psi _n(z)`$, because $`|\zeta _n(z)|`$ is either slowly changing (in the oscillating regime) or quickly decreasing with decreasing $`n`$ (in the exponential regime). This allows for taking even very rough estimates for the initial values of $`\psi _n(z)`$ in the downward recursion and the ratio $`\stackrel{~}{\psi }_{n1}(z)/\stackrel{~}{\psi }_n(z)`$ will still quickly approach the exact value $`A_n(z)`$. On the other hand, a downward recursion is not appropriate for computing $`\zeta _n(x)`$ or the ratio $`B_n(x)`$ if it starts in the exponential regime.

Now let us study the influence of the numerical error on the $`a_n,b_n`$ coefficients, and hence on $`Q_e`$ if an upward recursion is used to compute $`\psi _n(x)`$. In this analysis $`\zeta _n(x)`$ and $`B_n(x)`$ are assumed to be exact due to previous considerations. On the other hand $`A_n(y)`$ is also assumed to be exact. The effect of using approximate values of $`\psi _n(y)`$ will be considered later. We can make the discussion for $`a_n`$. Similar conclusions will hold for $`b_n`$. Eqn. (6) can be rewritten as

$$a_n=\frac{\psi _n(x)}{\zeta _n(x)}f(A_n(x)),$$
$`(23)`$
where only the $`A_n(x)`$ dependence is shown explicitly as it is the only relevant one for error analysis. The relative error in $`a_n`$ will be given by

$$\frac{\delta a_n}{a_n}\frac{\delta \psi _n}{\psi _n}+\frac{f^{}}{f}\frac{\delta A_n}{A_n}.$$
$`(24)`$
Recalling the definition (7), the relative error in $`A_n`$ can be estimated to be of the same order of magnitude as that of $`\psi _n`$, and taking into account that $`f`$ is a smooth function of the order of unity (cf. eqn. (10)), one gets the estimate

$$\delta a_na_n\frac{\delta \psi _n}{\psi _n}a_n\eta ^{}\frac{\zeta _n}{\psi _n}=\eta ^{}f\eta ^{}.$$
$`(25)`$
where use has been made of eqn. (22) and $`\eta ^{}`$ is some typical value of $`\eta _n^{}`$.

This means that the absolute error in $`a_n`$ or $`b_n`$, remains roughly constant throughout the computation. Of course eqn. (24) holds only for small $`\delta \psi _n`$, but this is guaranteed as $`N`$ is of the order of $`x`$ and so the recurrence does not go deep inside the exponential region. The important consequence of eqn. (25) is that the upward recursion can be used to obtain $`\psi _n(x)`$ because the error introduced is of the order of the roundoff error (see however the comment at the end of Section 6). Let us note that this fact is consistent with available algorithms for doing Mie calculations, where $`\psi _n(x)`$ and $`\zeta _n(x)`$ are always computed by upward recursion (e.g. refs. 5,11).

Let us consider now the effect of the numerical error coming form $`\psi _n(y)`$. We have argued before that an upward recursion would not be appropriate for computing $`\psi _n(z)`$ in general, however we have just shown that it can be used in the case of $`\psi _n(x)`$. The reason for this was that the relative error in $`\psi _n(x)`$ grew as $`\zeta _n(x)/\psi _n(x)`$ but the quantities $`a_n`$ and $`b_n`$ themselves converged to zero as $`\psi _n(x)/\zeta _n(x)`$. Both factors cancel rendering $`\delta a_n`$ and $`\delta b_n`$ bounded. We cannot apply a similar argument to $`\delta \psi _n(y)`$ and therefore an upward recursion is not reliable to compute $`\psi _n(y)`$ for arbitrary $`y`$. We can consider two limiting cases

a) $`m_I=0`$. In this case $`y`$ is real and greater than $`x`$, thus the instability in $`\psi _n(y)`$ starts only after that in $`\psi _n(x)`$, therefore the upward recursion can be used.

b) Large $`m_I`$. From the initial values4

$$\begin{array}{cc}\hfill \psi _0(z)& =\mathrm{sin}(z),\psi _1(z)=\frac{1}{z}\mathrm{sin}(z)\mathrm{cos}(z)\hfill \\ \hfill \zeta _0(z)& =\mathrm{i}\mathrm{exp}(\mathrm{i}z),\zeta _1(z)=\left(\frac{\mathrm{i}}{z}1\right)\mathrm{exp}(\mathrm{i}z)\hfill \end{array}$$
$`(26)`$
one can see that $`\psi _n\mathrm{exp}(m_Ix),\zeta _n\mathrm{exp}(m_Ix)`$, for small $`n`$, thus $`\psi _n`$ is much larger than $`\zeta _n`$. On the other hand $`ϵ_{0,1}`$ are related to the computer precision, typically $`ϵ_{0,1}r\psi _{0,1}`$ with $`r10^{16}`$ in double precision. Upon substitution in (21) we find that $`\eta `$ is small but $`\eta ^{}r\mathrm{exp}(2m_Ix)`$ which is not necessarily small. The relative error in $`\stackrel{~}{\psi }_n(z)`$ goes as

$$\frac{\delta \psi _n(z)}{\psi _n(z)}r\left|\frac{\psi _0(z)}{\zeta _0(z)}\right|\left|\frac{\zeta _n(z)}{\psi _n(z)}\right|$$
$`(27)`$
For small $`n`$ the relative error is small, of the order of $`r`$, however for $`n|z|`$, where $`\psi _n`$ and $`\zeta _n`$ are of the order of unity, the relative error is $`r|\psi _0/\zeta _0|r\mathrm{exp}(2m_Ix)`$ which is large for large $`m_I`$. Therefore the upward recursion is not stable in this case.

To summarize, the upward recursion to compute $`\psi _n(y)`$ can be used if $`m_I`$ is small enough but becomes unstable for large $`m_I`$. We have not analyzed in any detail in which cases the upward recursion for $`\psi _n(y)`$ is reliable, therefore we shall only consider downward recurrences for this quantity. See however refs. 10,11 for an extensive analysis of this problem through computer experiments. Noting that all we need is the ratio $`A_n(y)`$, for $`1nN`$, we can use the downward recursion

$$A_n(y)=F_n(y)\frac{1}{A_{n+1}(y)}.$$
$`(28)`$
Computing the initial value $`A_N(y)`$ requires some algorithm such as that of Lentz17 or the one we present in the next section. Let us estimate now the precision required in $`A_N(y)`$ in order not to introduce an error in $`Q_e`$ larger than the prescribed precision $`ϵ`$. By arguments similar to those used for $`\psi _n(x)`$, we have

$$\frac{\delta a_n}{a_n}\frac{\delta A_n(y)}{A_n(y)}$$
$`(29)`$
where $`\delta a_n`$ is the error introduced by $`\delta A_n(y)`$. Given that the downward recursion is stable we can assume that

$$\left|\frac{\delta a_n}{a_n}\right|\left|\frac{\delta A_N(y)}{A_N(y)}\right|\mathrm{for}nN$$
$`(30)`$
Using this relationship one gets for the numerical error in $`Q_e`$

$$\delta Q_e\frac{1}{x^2}\underset{n=1}{\overset{N}{}}(2n+1)\delta a_nQ_e\left|\frac{\delta A_N(y)}{A_N(y)}\right|.$$
$`(31)`$
Therefore the numerical error from $`A_n(y)`$ will be under control by imposing

$$\left|\frac{\delta A_N(y)}{A_N(y)}\right|\frac{ϵ}{Q_e}.$$
$`(32)`$
Let us note that this criterion will be conservative in general. An exception would be the case of $`y`$ being real and bigger than $`N`$. In this case the recurrence (28) has no healing properties (for it already starts in the oscillatory regime) and hence the equal sign is reached in (30).

4. INITIALIZATION OF THE DOWNWARD RECURRENCE

In this section we present a new method to compute $`A_N(z)`$, of similar efficiency to that due to Lentz17 (actually ours needs one multiplication less at each step). This method has the advantage of being able to implement a precision condition as that in eqn. (32), hence controlling the required precision in $`A_n(y)`$.

Let $`X_n(z)`$ and $`Y_n(z)`$ be two sequences satisfying the recurrence (4) for some value of $`z`$ (the dependence on $`z`$ is irrelevant here). Then they will satisfy the Wronskian identity

$$C=X_nY_{n+1}X_{n+1}Y_n$$
$`(33)`$
where $`C`$ is independent of $`n`$. We can rewrite it as a difference equation

$$C=Y_nY_{n+1}\left\{\left(\frac{X}{Y}\right)_n\left(\frac{X}{Y}\right)_{n+1}\right\},$$
$`(34)`$
and solve it in $`X_n`$

$$X_n=DY_n+CY_n\underset{k=n}{\overset{\mathrm{}}{}}(Y_kY_{k+1})^1,$$
$`(35)`$
$`D`$ being a constant. To write (35) we have assumed that $`Y_n`$ is a sequence going to infinity for large $`n`$, which is true for almost any solution of the recurrence (4). If we take $`Y_n`$ as a fixed sequence and regard $`C,D`$ as free parameters, then $`X_n`$ is the most general solution of the recurrence relation (4). In particular for $`D=0`$$`X_n`$ goes to zero as $`n`$ goes to infinity, as a consequence it must be proportional to $`\psi _n`$,

$$\psi _n(z)=C(z)Y_n(z)\underset{k=n}{\overset{\mathrm{}}{}}(Y_k(z)Y_{k+1}(z))^1$$
$`(36)`$
The constant $`C`$ cancels after computing the ratio $`A_n(z)`$

$$A_n(z)=Y_n^1\left\{Y_{n1}+Y_n^1\left[\underset{k=n}{\overset{\mathrm{}}{}}(Y_kY_{k+1})^1\right]^1\right\}.$$
$`(37)`$
Finally, a simpler formula can be obtained for $`A_N(z)`$ by choosing as starting values for the sequence $`Y_n`$

$$Y_{N1}=0,Y_N=1$$
$`(38)`$
$$A_N(z)=\left[\underset{k=N}{\overset{\mathrm{}}{}}(Y_k(z)Y_{k+1}(z))^1\right]^1.$$
$`(39)`$
About the convergence of the series in (39), we note that it is very fast when $`Y_k`$ enters in its exponential regime. Note that for real $`y`$ the convergence begins only after $`ky`$. A similar conclusion was reached by other authors11 in Lentz’s method which basically follows the same principle as ours and so has similar convergence properties.

The sequence in eqn. (39) must be truncated at some value $`k=M`$ in such a way as to fulfill the requirement (32). This can be easily done by noting that the error introduced in $`A_N^1(y)`$ is of the order of the last term taken into account (this follows from $`|Y_k/Y_{k1}||F_k|>2`$ for large $`k`$),

$$\delta A_N^1\left(Y_MY_{M+1}\right)^1.$$
$`(40)`$
On the other hand we should require

$$|\delta A_N^1(y)|\left|A_N^1\frac{\delta A_N}{A_N}\right|\left|\frac{1}{F_N(y)}\frac{ϵ}{Q_e}\right|$$
$`(41)`$
where we have made use of eqns. (8) and (32). Recall now that for $`x1`$, $`F_N`$ and $`Q_e`$ are of the order of unity whereas for $`x1`$ the product of $`F_NQ_e`$ is still of the order of unity, therefore the final criterion to truncate (39) is

$$\left|\left(Y_M(y)Y_{M+1}(y)\right)^1\right|ϵ.$$
$`(42)`$
To finish this section we shall show how to avoid ill-conditioning in (39), which will appear if $`Y_k`$ gets too near to zero for some value of $`k`$. To do this we can use the recurrence relation (4) to write

$$\frac{1}{Y_{k1}Y_k}+\frac{1}{Y_kY_{k+1}}=\frac{Y_{k1}+Y_{k+1}}{Y_{k1}Y_kY_{k+1}}=\frac{F_k}{Y_{k1}Y_{k+1}},$$
$`(43)`$
which is well behaved even for $`Y_k=0`$.

5. COMPUTATIONAL ALGORITHM

Using the previous ideas, we have developed a computational algorithm which we shall briefly describe now. The input is $`x,m`$ and $`ϵ`$ and the main output are the coefficients $`a_n`$ and $`b_n`$, and $`N`$. To start with, analytic expressions for $`\zeta _0(x)`$ and $`\zeta _1(x)`$ are taken to initiate an upward recurrence for $`\zeta _n(x)`$. This quantity is kept in a (complex) array variable. The recurrence stops when the condition (16) is fulfilled, providing the value of $`N`$. The quantities $`\psi _n(x)`$ are automatically obtained as the real part of $`\zeta _n(x)`$. As a second step, $`A_N(y)`$ is computed using eqns. (38), (39) and (42). Here we note that from a computational point of view an equivalent form of (42) is more convenient, which consist in doing the check for the absolute values of the real and imaginary parts. This is much faster than computing the modulus of a complex number.

Then a downward recurrence is performed for $`A_n(y)`$, eqn. (28), until $`n=1`$. Simultaneously, $`a_n`$ and $`b_n`$ are computed using $`\zeta _n(x)`$ and $`A_n(y)`$. The quantities $`Q_s`$ and $`Q_e`$ can then be computed. We have not developed any especial algorithm for computing the scattering amplitudes $`S_1`$ and $`S_2`$. To do this efficiently see for instance ref. 11.

The criteria developed above are intended to be robust, hence they are rather conservative. As a consequence the error in $`Q_e`$ is smaller than the prescribed precision $`ϵ`$. This is especially true for small values of $`x`$, whereas for $`x1`$, about two more figures than expected are obtained. We point out also that $`Q_s`$ is always obtained as accurately as $`Q_e`$ or more. This fact was expected because the criteria were stated for $`|a_n|`$ while $`Q_s`$ goes as $`|a_n|^2`$ which converges faster.

6. RESONANT TERMS IN THE MIE SERIES

Let us recall that after eqn. (7) we stated a smoothness assumption for the quantities $`[a]_n,[b]_n`$, namely that they are nearly constant in the $`x`$ exponential regime and do not play any role in the convergence of the Mie series, which is only controlled by the ratio $`\psi _n(x)/\zeta _n(x)`$. In particular this assumption implied that the highest partial wave with a relevant contribution is independent of $`m`$ (cf. eqn. (16)). In other words, $`N`$ is a function of $`x`$ only. This result is also supported numerically, (see for instance refs. 10,11). Therefore it was a surprise for us to discover that strictly speaking such a statement must be false. Moreover, for any choice of $`N`$ as a function of $`x`$ only, and for any prescribed value of n, $`n>N`$, one can always pick a value of $`m`$ (in fact infinitely many of them) in such a way that the $`n`$-th term in the Mie series is not negligible, for instance one can make $`a_n=1`$. The consequence of this that in order to guarantee that the numerical value of $`Q_e`$ is correct within some prescribed precision, $`N`$ should depend on $`m`$ as well as on $`x`$.

In order to clarify the point let us consider the worst case, which is also the simplest, namely $`m_I=0`$, i.e. $`y`$ real. This is the only case in which $`|a_n|`$ or $`|b_n|`$ can reach the value 1. The point can be made for $`a_n`$: recalling that for $`z`$ real Re$`\zeta _n(z)=\psi _n(z)`$, eqn. (2) can be rewritten as

$$\begin{array}{cc}\hfill a_n& =\frac{\mathrm{Re}D_n}{D_n}\hfill \\ \hfill D_n& =x\zeta _n(x)\psi _n^{}(y)y\zeta _n^{}(x)\psi _n(y)\hfill \end{array}$$
$`(44)`$
where $`D_n`$ is a complex quantity. Obviously $`a_n=1`$ if and only if

$$\mathrm{Im}D_n=0.$$
$`(45)`$
Let us regard $`x`$ and $`n`$ as given and look for solutions of (45) in the variable $`y`$. The equation can be rewritten as

$$\frac{1}{y}\frac{\psi _n^{}(y)}{\psi _n(y)}=\frac{1}{x}\frac{\mathrm{Im}\zeta _n^{}(x)}{\mathrm{Im}\zeta _n(x)}$$
$`(46)`$
In the interval $`y>n`$, $`\psi _n(y)`$ is a real oscillating function of $`y`$ with infinitely many zeroes. Between two zeroes of $`\psi _n(y)`$, the l.h.s. of eqn. (46) takes every real value, therefore there are infinitely many solutions to our equation for any values of $`x`$ and $`n`$, no matter how large is $`n`$ as compared to $`x`$. For these values of $`x,m`$, and $`n`$, $`a_n`$ will not at all be negligible.

Let us now show that these resonances do not occur for unrealistic values of $`m`$. Typically (and asymptotically for large $`y`$) the distance between two consecutive zeroes of $`\psi _n(y)`$ is of the order of $`\pi `$, therefore for given $`x`$ and $`n`$ the lowest resonant value of $`m`$ will occur near the interval $`(\frac{n}{x},\frac{n+\pi }{x})`$ approximately. For large $`x`$ this happens for $`m`$ near to unity, and all the other resonant values will follow at a distance of about $`\pi /x`$ from each other.

From a rigorous point of view these findings would invalidate the estimates (10) and their consequences. They would also invalidate any algorithm in which $`N`$ depends on $`x`$ only, namely every existent algorithm known to us. In fact the only practical way to make sure that the resonant partial waves have been accounted for would be to take $`N`$ greater than $`y`$ in order to guarantee that $`\psi _n(y)`$ has no zeroes for $`n>N`$.

Nevertheless it is clear that in practice the existent algorithms to do Mie scattering calculations work. To account for this fact we should consider not only the existence of resonant partial waves but also their width. Let us show that for sensible choices of $`N`$ (as a function of $`x`$) and for $`n>N`$ the resonances are so narrow that they will not normally show up. Let $`y_0`$ be one the values of $`y`$ such that $`a_n=1`$. A look to eqn. (44) shows that for generic $`y`$, Re$`D_n`$ goes as $`\psi _n(x)`$ whereas $`D_n`$ goes as $`\zeta _n(x)`$, therefore $`a_n`$ is very small. However for the especial value $`y_0`$ there is a cancellation between two huge numbers in Im$`D_n`$, leaving $`a_n`$ of the order of unity. The range of values of $`y`$ for which a partial cancellation takes place is related to the slope of $`D_n`$ in $`y=y_0`$, namely

$$\mathrm{\Gamma }\left|\frac{D_n}{D_n^{}}\right|_{y=y_0}=\left|\frac{\mathrm{Re}D_n}{D_n^{}}\right|_{y=y_0}\left|\frac{\psi _n(x)}{\zeta _n(x)}\right|.$$
$`(47)`$
Where $`D_n^{}=\mathrm{d}D_n/\mathrm{d}y`$. In other words, if $`N`$ is large enough only by a very careful choice of $`m`$ or $`x`$ can one find one these resonant contributions. More precisely, recalling eqn. (13), we can see that $`m`$ or $`x`$ should be fine tuned at least with a precision $`ϵ`$ in order to pick a resonant term for some $`n>N`$. On the other hand, except for these rare cases, $`a_n,b_n`$ are indeed small and of the order of $`\psi _n(x)/\zeta _n(x)`$, therefore our analysis applies. If $`m`$ is allowed to be complex, a more involved analysis would be needed, but we expect that the conclusion would not differ.

Let us finally note another consequence of the resonant terms on the calculation, even when they are taken into account. For one of these terms the quantity $`f`$ in eqn. (23) is no longer of the order of unity, on the contrary it is rather large, and the last step in eqn. (25) cannot be taken. This means that a resonant term amplifies the error due to $`\psi _n(x)`$. The cure is simply to compute $`\psi _n(x)`$ by downward recursion for $`x<n<N`$. This has in fact been observed in selected quantities such as the backscattering efficiency for suitable values of $`x`$ and $`m`$ (Ref. 5).

7. CONCLUSIONS

In this paper we have addressed several points relevant to Mie scattering calculations. To be specific:

a) We have estimated the error introduced in the calculation by truncating the Mie series, thereby finding a prescription for choosing $`N`$. We have found that in the generic case $`N`$ depends on $`x`$ only.

b) The possible instabilities in the recursions used to compute $`\psi _n`$ and $`\zeta _n`$ have been analyzed. We have found that upward recursion is always unstable for computing $`\psi _n(z)`$ if $`n`$ is large enough. However it can be used to compute $`\psi _n(x)`$ in Mie calculations. As a matter of fact $`\psi _n(x)`$ is computed in this way in nowadays available algorithms. We have also found that upward recursion can be used for $`\psi _n(y)`$ if $`m_I`$ is small enough, but no criterion is given for how small $`m_I`$ should be.

c) A criterion has been established for the allowed error in $`\psi _{n1}(y)/\psi _n(y)`$.

d) A new method to compute $`\psi _{n1}(y)/\psi _n(y)`$ is presented which is efficient and allows for controlling the error and removing ill-conditioning.

e) It has been shown the existence of resonant terms in the Mie series which can also appear for $`n>N`$. Strictly speaking the existence of these terms invalidates any algorithm in which $`N`$ is a function of $`x`$ only. However we have also shown that those resonant terms are extremely rare, namely they appear with a probability of the order of $`ϵ`$.

A specific algorithm is also described. It is meant to be robust and efficient for a wide range of size parameters and refractive indices. With this algorithm we have written the computer program LVEC-MIE18, which is available both in single and double precision contacting V.E. Cachorro.

APPENDIX I

Let us justify the bound (11). To do so we shall study the convergence rate of the terms left out in the series, $`n>N`$. In this region we can make use of the estimate (10),

$$\begin{array}{cc}\hfill \delta Q_e& =\frac{2}{x^2}\underset{n=N+1}{\overset{\mathrm{}}{}}(2n+1)\mathrm{Re}(a_n+b_n)\hfill \\ & \frac{2}{x^2}\underset{n=N+1}{\overset{\mathrm{}}{}}(2n+1)\left(|a_n|+|b_n|\right)\hfill \\ & \frac{8}{x^2}\underset{n=N+1}{\overset{\mathrm{}}{}}n\left|\frac{\psi _n(x)}{\zeta _n(x)}\right|\hfill \\ & =\frac{8}{x^2}\underset{n=N+1}{\overset{\mathrm{}}{}}n\left|\frac{B_n(x)}{A_n(x)}\frac{B_{n1}(x)}{A_{n1}(x)}\mathrm{}\frac{B_{N+1}(x)}{A_{N+1}(x)}\right|\left|\frac{\psi _N(x)}{\zeta _N(x)}\right|.\hfill \end{array}$$
$`(\mathrm{I}.1)`$
Now making use of (8) and recalling that $`F_n(x)`$ is a monotonically increasing function of $`n`$, we obtain

$$\begin{array}{cc}\hfill \delta Q_e& \frac{8}{x^2}\underset{n=N+1}{\overset{\mathrm{}}{}}n\frac{1}{F_n^2(x)}\frac{1}{F_{n1}^2(x)}\mathrm{}\frac{1}{F_{N+1}^2(x)}\left|\frac{\psi _N(x)}{\zeta _N(x)}\right|\hfill \\ & \frac{8}{x^2}\underset{n=N+1}{\overset{\mathrm{}}{}}n\left(F_N(x)\right)^{2(Nn)}\left|\frac{\psi _N(x)}{\zeta _N(x)}\right|\hfill \\ & =\frac{8}{x^2}\left(\frac{N}{F_N^2(x)1}+\frac{F_N^2(x)}{\left(F_N^2(x)1\right)^2}\right)|a_n|\hfill \end{array}$$
$`(\mathrm{I}.2)`$
For small $`x`$, $`N=1`$ and $`F_N(x)`$ is large, hence

$$\delta Q_e2|a_n|$$
$`(\mathrm{I}.3)`$
on the other hand, for large $`x`$, $`Nx`$ and $`F_N(x)2`$,

$$\delta Q_e\frac{8}{3}\frac{1}{x}|a_n|.$$
$`(\mathrm{I}.4)`$
In both cases eqn. (11) is valid (up to factors of the order of unity).

APPENDIX II

In order to know in advance the value of $`N`$ that will be obtained from the prescription (16) for given $`x`$ and $`ϵ`$, let us recall that Im$`\zeta _n(x)=\sqrt{\pi x/2}Y_{N+\frac{1}{2}}(x)`$$`Y_\nu (z)`$ being the Bessel function of the second kind. Let $`\nu `$ and $`c`$ be defined by

$$\begin{array}{cc}\hfill N& =\nu \frac{1}{2}\hfill \\ \hfill x& =\nu c\nu ^{1/3}.\hfill \end{array}$$
$`(\mathrm{II}.1)`$
Note that for large $`\nu `$, eqn. (II.1) can be inverted to give $`Nx+cx^{1/3}`$. Now we can make use of the leading order term in the asymptotic expansion of $`Y_\nu `$ for large $`\nu `$ and fixed $`c`$, ref. 4, p. 367:

$$\mathrm{Im}\zeta _N(x)\sqrt{\pi }\left(\frac{\nu }{2}\right)^{1/6}\mathrm{Bi}\left(2^{1/3}c\right),$$
$`(\mathrm{II}.2)`$
where Bi$`(z)`$ is the Airy function of the second kind, ref. 4, p. 446. This function is given by

$$\mathrm{Bi}(z)=z^{1/4}f(z)\mathrm{exp}(\frac{2}{3}z^{3/2}),$$
$`(\mathrm{II}.3)`$
where $`f(z)`$ is nearly constant for $`z>1`$ with $`f(z)1/\sqrt{\pi }`$, ref. 4, p. 449. Thus

$$\left|\mathrm{Im}\zeta _N(x)\right|\left(\frac{\nu }{2\sqrt{2}}\right)^{1/6}c^{1/4}\mathrm{exp}\left(\frac{1}{3}(2c)^{3/2}\right).$$
$`(\mathrm{II}.4)`$
The right hand side of (II.4) has a very strong dependence on $`c`$ whereas it depends very smoothly on $`\nu `$. Actually $`(\nu /2\sqrt{2})^{1/6}`$ is of the order of unity for $`\nu =1`$ up to $`10^5`$. Therefore using eqn. (16), $`c`$ will be determined by $`ϵ`$. We find that $`c=4.3`$ corresponds to $`ϵ=10^8`$. Other values are $`c=4.0,ϵ=10^7`$, and $`c=5.0,ϵ=10^{10}`$, computed for $`\nu =100`$ in (II.4).

REFERENCES

1. H. C. van de Hulst, Light Scattering by Small Particles, John Wiley, N. Y. 1957.

2. M. Kerker, The Scattering of Light and Other Electromagnetic Radiation, Academic Press. N. Y., 1969.

3. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles, Wiley Interscience, N. Y. 1983.

4. M. Abramowitz and I. A. Stegun ed., Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables, Dover Pub. Inc., N. Y., 1965.

5. J. V. Dave, Subroutines for Computing the Parameters of Electromagnetic Radiation Scattered by a Sphere, Report No. 320-3237, IBM Scientific Center, Palo Alto, California, USA, 1968.

6. J. V. Dave, Scattering of Electromagnetic Radiation by Large Absorbing Spheres, IBM J. Res. Develop., Vol.13, 1302-1313, 1969.

7. W. Blattner, Utilization Instruction for Operation of the Mie Programs on the CDC-6600 Computer at AFCRL, Radiation Center Associates, Ft. Worth, Texas, Res. Note RRA-N7240, 1972.

8. G. Grehan and G. Gouesbet, The Computer Program SUPERMIDI for Mie Theory Calculation, without Practical Size nor Refractive Index Limitations, Internal Report TTI/GG/79/03/20, Laboratoire de Génie Chemique Analytique, U. de Rouen, 76130 Mt-St-Aignan (France), 1979. Also Private communication.

9. G. Grehan and G. Gouesbet, Mie theory calculations: new progress, with emphasis on particle sizing, Appl. Opt., Vol. 18, 3489-3493, 1979.

10. W. J. Wiscombe, Mie scattering calculations: Advances in technique and fast vector speed computer codes. NCAR Technical Note NCAR/TN-140+STR (National Center for Atmospheric Research) Boulder, Colorado, 80307, 1979, and private communication.

11. W. J. Wiscombe, Improved Mie Scattering Algorithms, Appl. Opt., Vol. 19, 1505-1509, 1980.

12. G. H. Goedecke, A. Miller and R. C. Shirkey, Simple Scattering Code Agausx, in Atmospheric Aerosols: Their Formation, Optical Properties and Effects. Ed. A. Deepak, Spectrum Press, Hampton, Virginia, 1982.

13. A. Miller, Comments on Mie Calculations, Am. J. Phys., Vol. 54, 297-297, 1986. Also private communication.

14. G. W. Kattawar and G. N. Plass, Electromagnetic Scattering from Absorbing

Spheres, Appl. opt., vol. 6, 1377, 1967.

15. D. Deirmendjian, Electromagnetic Scattering on Spherical Polydispersion, Elsevier, N. Y. 1969.

16. H. Quenzel and H. Müller, Optical properties of single particles diagrams of intensity, extinction scattering and absorption efficiencies, Wissenschaftliche Mitteilung, n. 34. Metereologisches Institut, Universität München, 1978.

17. W. J. Lentz, Generating Bessel Functions in Mie Scattering Calculations using Continued Fractions, Appl. Opt., vol. 15, 668-671, 1976.

18. V. E. Cachorro, L. L. Salcedo and J. L. Casanova, Programa LVEC-MIE para el cálculo de las magnitudes de la teoría de esparcimiento de Mie, Anales de Física, vol. 85, Serie B, 198-211, 1989.


text:
# Galactic Bulge Pixel Lensing Events

## 1 Introduction

The previous and the current microlensing surveys towards the Galactic bulge and the Magellanic Clouds (Alcock et al., 1993; Aubourg et al., 1993; Udalski et al., 1993; Alard & Guibert, 1997; Abe et al., 1997) have been and are being carried out by using the conventional photometric method based on the extraction of the individual source stars’ point spread function (PSF). With this method, the lensing-induced light variations can be measured only for events with resolved source stars, and thus events occurred to unresolved faint sources cannot be detected. Pixel lensing, on the other hand, detects and measures the variation of a source star flux by subtracting successive images of the same field from a normalized reference image (Tomaney & Crotts, 1996; Alard & Lupton, 1998; Alard, 1999), enabling one to measure the light variations even for unresolved source stars. The pixel lensing technique was originally proposed to search for lensing events towards the unresolved star field of M31 (Crotts, 1992; Baillion et al., 1993). Recently, pixel lensing searches have been conducted towards the resolved Galactic bulge field (Alcock et al., 1999a, b, 2000), but they are intended to make use of the additional sources lying below the detection threshold, and so to increase the sensitivity of the existing experiments.

Gould & DePoy (1998) proposed a radically different type of pixel lensing survey towards the Galactic bulge using a small aperture ($`65`$ mm) camera with a large pixel-size ($`10^{\prime \prime }`$) detector and deliberately degraded optics achieving $`30^{\prime \prime }`$ PSF. With these instruments, the observations are analogous to the normal pixel lensing observations of the bulge of M31, but are carried out under conditions where the detected events can be followed up in detail. Since the proposed PSF set by the optics will not vary in time with atmospheric conditions, the data reduction process will be very simple. Compared to classic lensing surveys based on PSF photometry, the proposed pixel lensing survey has an advantage of the broad coverage of most or all of the ($`50\mathrm{deg}^2`$) southern bulge field of interest with the large field of view ($`36\mathrm{deg}^2`$ with a 2K$`\times `$2K CCD camera). In addition, the instrumentation required to carry out the survey will be inexpensive. Gould & DePoy (1998) described the possibility of the Galactic bulge pixel lensing survey in general, but detailed estimation of the event rate and investigation about the characteristics of the detectable events have not been done.

In this paper, we estimate the event rate of the proposed Galactic pixel lensing survey expected under various detection criteria and investigate the characteristics of the detectable events.

## 2 Event Rate Estimation

### 2.1 Event Production

To estimate the event rate, we perform simulations of Galactic bulge pixel lensing events by producing a large number of artificial light curves. The light curve of a pixel lensing event is represented by

$$F=F_0(A1),$$
(1)
$$A=\frac{u^2+2}{u\sqrt{u^2+4}};u=\left[\beta ^2+\left(\frac{t}{t_\mathrm{E}}\right)^2\right]^{1/2},$$
(2)
where $`F_0`$ is the baseline flux of the source star, $`A`$ is the amplification, $`u`$ is the lens-source separation (normalized by the angular Einstein ring radius $`\theta _\mathrm{E}`$), $`t_\mathrm{E}`$ is the Einstein ring radius crossing time (Einstein timescale), and $`\beta `$ is the impact parameter (also normalized by $`\theta _\mathrm{E}`$). The Einstein timescale is related to the physical parameters of the lens system by

$$t_\mathrm{E}=\frac{r_\mathrm{E}}{v};r_\mathrm{E}=D_{\mathrm{ol}}\theta _\mathrm{E}=\left[\frac{4Gm}{c^2}\frac{D_{\mathrm{ol}}(D_{\mathrm{os}}D_{\mathrm{ol}})}{D_{\mathrm{os}}}\right]^{1/2},$$
(3)
where $`r_\mathrm{E}`$ is the Einstein ring radius, $`m`$ is the lens mass, $`v`$ is the lens-source transverse speed, and $`D_{\mathrm{ol}}`$ and $`D_{\mathrm{os}}`$ represent the distances to the lens and source from the observer, respectively. Therefore, to produce a pixel lensing event light curve, it is required to define the lensing parameters of $`F_0`$, $`\beta `$, $`m`$, $`D_{\mathrm{os}}`$, $`D_{\mathrm{ol}}`$, and $`v`$.

We produce light curves by selecting the lensing parameters so that the event rate, $`\mathrm{\Gamma }`$, becomes

| $`\mathrm{\Gamma }`$ | $`{\displaystyle 𝑑D_{\mathrm{os}}\rho (D_{\mathrm{os}})𝑑D_{\mathrm{ol}}\rho (D_{\mathrm{ol}})𝑑v_y𝑑v_zvf(v_y,v_z)}`$ |
| --- | --- |
|  | $`{\displaystyle 𝑑m\mathrm{\Phi }(m)r_\mathrm{E}𝑑L\mathrm{\Phi }(L)𝑑\beta f(\beta )},`$ |

(4)
where $`\rho (D_{\mathrm{ol}})`$ and $`\rho (D_{\mathrm{os}})`$ are the density distributions of the lens and source locations, $`(v_y,v_z)`$ are the two components of the transverse velocity, $`f(v_y,v_z)`$ is their distribution, $`\mathrm{\Phi }(m)`$ is the lens mass function, $`\mathrm{\Phi }(L)`$ is the absolute luminosity function of source stars, and $`f(\beta )`$ is the distribution of impact parameters. In eq. (4), the factors $`r_\mathrm{E}`$ and $`v=(v_y^2+v_z^2)^{1/2}`$ are included because events with larger cross-sections (i.e. $`r_\mathrm{E}`$) and higher transverse speeds have higher chances of being lensed.

For the source and lens density distributions, we adopt a barred bulge model of Dwek et al. (1995) and a double-exponential disk model of Bahcall (1986). The Dwek et al. bulge model is represented by

$$\rho =\rho _0\mathrm{exp}\left[0.5r_s^2\right]M_{}/\mathrm{pc}^3,$$
(5)
where $`r^4=[(x^{}/x_0)^2+(y^{}/y_0)^2]^2+(z^{}/z_0)^4`$, $`(x^{},y^{},z^{})`$ are the three axes of the triaxial bulge ($`x^{}`$ is the longest and $`z^{}`$ is the shortest), and $`(x_0,y_0,z_0)=(1.58,0.62,0.43)\mathrm{kpc}`$ are the scale lengths. The normalization constant $`\rho _0`$ is set so that the total bulge mass becomes $`1.8\times 10^{10}M_{}`$ (Han & Gould, 1995). The Bahcall disk model is represented by

$$\rho =0.06\mathrm{exp}\left[\left(\frac{RR_0}{h_R}+\frac{z}{h_z}\right)\right]M_{}/\mathrm{pc}^3,$$
(6)
where $`h_R=3.5\mathrm{kpc}`$ and $`h_z=325\mathrm{pc}`$ are the radial and vertical scale heights of the disk and $`R_0=8\mathrm{kpc}`$ is the galactocentric distance of the sun.

For the transverse velocity distribution, we adopt the model of Han & Gould (1995). In this model, the velocity distributions for both disk and bulge components have a Gaussian form, i.e.

$$f(v_i)\mathrm{exp}\left[\frac{(v_i\overline{v}_i)^2}{2\sigma _i^2}\right];iy,z,$$
(7)
where the means and the standard deviations of the individual velocity components are

$$(\overline{v}_y,\sigma _y)=\{\begin{array}{cc}(220.0,30.0)\mathrm{km}\mathrm{s}^1\hfill & \text{(disk)}\hfill \\ (0.0,82.5)\mathrm{km}\mathrm{s}^1\hfill & \text{(bulge),}\hfill \end{array}$$
(8)
$$(\overline{v}_z,\sigma _z)=\{\begin{array}{cc}(0.0,20.0)\mathrm{km}\mathrm{s}^1\hfill & \text{(disk)}\hfill \\ (0.0,66.3)\mathrm{km}\mathrm{s}^1\hfill & \text{(bulge).}\hfill \end{array}$$
(9)
Unlike the models of the matter density and velocity distributions that are constrained by observations, our knowledge about the lens mass function is very poor. Therefore, we start with several proposed models of the lens mass function and compute the resulting $`t_\mathrm{E}`$ distributions, $`d\mathrm{\Gamma }(t_\mathrm{E})/dt_\mathrm{E}`$, combined with the adopted models of the matter density and velocity distributions. Then the mass function model producing $`t_\mathrm{E}`$ distribution that best matches the observed one is used for our simulation. We test three mass function models. These models are a main-sequence mass function of Kroupa, Tout & Gilmore (1993) (with three powers of $`\alpha =4.5`$ for $`m>1.0M_{}`$, $`\alpha =2.5`$ for $`0.5M_{}<m1.0M_{}`$, and $`\alpha =1.2`$ for $`0.08M_{}<m0.5M_{}`$), a single power-law distribution of Han & Gould (1996) (with $`\alpha =2.3`$ and a mass cutoff of $`m_{\mathrm{cut}}=0.04M_{}`$), and another power-law distribution of Zhao, Spergel, & Rich (1995) (with $`\alpha =2.0`$ and $`m_{\mathrm{cut}}=0.08M_{}`$). In Figure 1, we present the distributions of Einstein timescales resulting from the model mass functions and they are compared with that of the latest $`t_\mathrm{E}`$ data of the 97 events detected by the MACHO group (Alcock et al., 1999b, 2000). One finds that the $`t_\mathrm{E}`$ distribution resulting from the Zhao et al. (1995) model is consistent with the observed one, while those resulting from the Han & Gould (1996) and Kroupa et al. (1993) models produce short timescale events that are more than and less than the observed one, respectively. Therefore, we adopt the Zhao et al. (1995) model for our simulation.

For the source brightness distribution, we adopt the $`I`$-band absolute luminosity function determined by Kroupa (1995). We assume that photometry of the pixel lensing survey is carried out in $`I`$ band because more light can be obtained due to its large band width and small amount of extinction compared to other optical bands. Once the location of a source star and its absolute brightness, $`M_I`$, are selected, the apparent source star brightness is obtained by calculating the corresponding distance modulus, $`\mu `$, and extinction, $`A_I`$, i.e. $`I=M_I+\mu (D_{\mathrm{os}})+A_I(D_{\mathrm{os}})`$. Due to the patchy and irregular characteristics of the Galactic interstellar absorbing matter distribution, it is difficult to model its detailed distribution. Hence, we approximate that the interstellar matter is distributed similar to the Galactic disk matter distribution. With this approximation, the source flux is decreased from its de-reddened value by a factor $`\mathrm{exp}\left[\kappa \mathrm{\Sigma }(D_{\mathrm{os}})\right]`$, where $`\mathrm{\Sigma }(D_{\mathrm{os}})`$ represents the column density of the absorbing matter located between the observer and source star and $`\kappa `$ is a normalization constant. Stanek (1996) determined that the extinction in $`V`$ band towards the Baade’s Window (BW) varies $`A_V1.26`$$`2.79`$ with the most frequent value of $`A_V1.7`$. With the relation $`A_V/E(VI)2.49`$ also provided by him, this corresponds to $`A_I0.76`$$`1.67`$ with $`A_I1.02`$. Therefore, we set the value $`\kappa `$ so that the average amount of extinction towards the BW becomes $`A_I=1.02`$ at a distance of $`D_{\mathrm{os}}=8\mathrm{kpc}`$. We discuss the effect of variable extinction with position on the sky in the following subsection.

Under the definition of a lensing event as ‘a close lens-source encounter within the Einstein ring of a lens’, the impact parameters are randomly chosen, i.e. $`f(\beta )\mathrm{constant}`$, in the range $`0\beta 1`$.

### 2.2 Efficiency Determination

Once events are produced, we then determine the efficiency of pixel lensing event detection by imposing the following detection criteria.

First, we require detectable events should have high signal-to-noise ratios during their peak amplifications. Since the uncertainty of the light variation measurement of a pixel lensing event is almost totally dominated by the background flux, the signal-to-noise ratio is computed by

$$S/NF_0(A1)\left(\frac{t_{\mathrm{exp}}}{F_0A+B}\right)^{1/2},$$
(10)
where $`t_{\mathrm{exp}}=2`$ min is the mean exposure time and $`B`$ is the average background flux. Determination of $`B`$ is described below.

Second, since event detection is also restricted by the duration of events, $`t_{\mathrm{dur}}`$, we additionally require detectable events should last long enough for solid confirmation. We define the duration as ‘the length of the time period during which the signal-to-noise ratio exceeds a threshold value of $`(S/N)_{\mathrm{th}}`$’ and it is computed by

$$t_{\mathrm{dur}}=2t_\mathrm{E}\left(\beta _{\mathrm{th}}^2\beta ^2\right)^{1/2},$$
(11)
where $`\beta _{\mathrm{th}}`$ represents the maximum allowed impact parameter for an event to satisfy the threshold signal-to-noise ratio (Han, 1999). To compute $`\beta _{\mathrm{th}}`$, we first compute the threshold amplification $`A_{\mathrm{th}}`$ corresponding to $`(S/N)_{\mathrm{th}}`$ by numerically solving eq. (10). With the obtained value of $`A_{\mathrm{th}}`$, then, the threshold impact parameter is obtained by

$$\beta _{\mathrm{th}}=\left[2(1A_{\mathrm{th}}^2)^{1/2}2\right]^{1/2}.$$
(12)
With these detection criteria, the detection efficiency is determined by obtaining the ratio of the number of events satisfying the above detection criteria out of the total number of tested events.

The proposed pixel lensing survey will be performed on a large area of the sky. Then, the event rate per unit angular area, $`d\mathrm{\Gamma }/d\mathrm{\Omega }`$, will be different for different regions of the sky because the stellar distribution is not uniform due to the variation both in geometric distance from the Galactic center and extinction. In addition, the spatial variation of the optical depth, $`\tau `$, and the disk/bulge event ratio with the position on the sky also makes the event rate depend on the observed field. To consider this effect, we construct maps of the surface brightness, $`\mu _I(l,b)`$, optical depth, $`\tau (l,b)`$, and the disk event fraction, $`f_{\mathrm{\Gamma },\mathrm{disk}}(l,b)`$. The surface brightness map is constructed by using the matter density distribution and extinction models described in § 2.1. The map is normalized so that the de-reddened surface brightness of the BW field contributed by bulge stars becomes $`17.6\mathrm{mag}\mathrm{arcsec}^2`$ (Terndrup, 1998) and the amount of the extinction towards the same field becomes $`A_I=1.02`$ mag. With the additional flux from disk stars, which is $`25\%`$ of the flux from bulge stars, and the sky background flux, which contributes $`20\%`$ of the total stellar background flux, then the average background flux within the PSF (with an angular area of $`\mathrm{\Omega }_{\mathrm{PSF}}0.75\mathrm{arcmin}^2`$) of the BW field becomes $`800\mathrm{e}^1\mathrm{s}^1`$. The optical depth represents the average probability for a single source to be gravitationally amplified by a factor $`A1.34`$ at a given moment and it is computed by

| $`\tau =`$ | $`{\displaystyle _0^{\mathrm{}}}𝑑D_{\mathrm{os}}\rho (D_{\mathrm{os}}){\displaystyle _0^{D_{\mathrm{os}}}}𝑑D_{\mathrm{ol}}\rho (D_{\mathrm{ol}})\pi r_\mathrm{E}^2`$ |
| --- | --- |
|  | $`\left[{\displaystyle _0^{\mathrm{}}}𝑑D_{\mathrm{os}}\rho (D_{\mathrm{os}})\right]^1.`$ |

(13)
The disk event fraction is approximated by $`f_{\mathrm{\Gamma },\mathrm{disk}}\eta f_{\tau ,\mathrm{disk}}`$, where $`f_{\tau ,\mathrm{disk}}`$ represents the ratio of disk event contribution to the optical depth and $`\eta 0.7`$ is included because disk lensing events tend to have longer $`t_\mathrm{E}`$ (Han & Gould, 1995) and thus $`f_{\mathrm{\Gamma },\mathrm{disk}}`$ is smaller than $`f_{\tau ,\mathrm{disk}}`$. In Figure 2, we present the maps of $`\mu _I(l,b)`$, $`\tau (l,b)`$, and $`f_{\mathrm{\Gamma },\mathrm{disk}}(l,b)`$.

In Figure 3, we present the baseline source star brightness distribution for pixel lensing events (effective source luminosity function, $`\mathrm{\Phi }_{\mathrm{eff}}`$). The presented $`\mathrm{\Phi }_{\mathrm{eff}}`$ is for the source stars of the events detectable towards the BW under the criteria of $`(S/N)_{\mathrm{th}}=10`$ and $`t_{\mathrm{dur},\mathrm{th}}=6`$ hr. It is normalized so that the surface flux contributed by all source stars matches the background stellar flux. In the figure, we also present the effective luminosity functions for events with high amplifications: light shaded region for events with $`A20`$ and dark shaded region for events with $`A200`$ (extreme microlensing events, EMEs).

### 2.3 Total Event Rate

With the obtained effective luminosity function, the total pixel lensing event rate is then calculated by

| $`\mathrm{\Gamma }_{\mathrm{tot}}={\displaystyle _{\mathrm{\Omega }_{\mathrm{tot}}}}{\displaystyle \frac{d\mathrm{\Gamma }}{d\mathrm{\Omega }}}𝑑\mathrm{\Omega };`$ |
| --- |
| $`{\displaystyle \frac{d\mathrm{\Gamma }}{d\mathrm{\Omega }}}={\displaystyle \frac{\pi }{2}}t_{\mathrm{obs}}N_{\mathrm{eff}}\tau _\mathrm{p}{\displaystyle \frac{1}{t_\mathrm{E}}},`$ |

(14)
where $`N_{\mathrm{eff}}=\mathrm{\Phi }_{\mathrm{eff}}(I)𝑑I`$ is the total number of effective source stars of a field, $`1/t_\mathrm{E}`$ is the mean value of the inverse Einstein timescales of the events detected towards the field, and $`\tau _\mathrm{p}`$ is the pixel lensing optical depth. Compared to the source stars of events detectable from classic lensing surveys (classic lensing events), those of pixel lensing events will be systematically brighter, and thus closer. This can be seen in Figure 4, where we plot the distributions of the source locations for classic, $`f_\mathrm{c}(D_{\mathrm{os}})`$, and pixel lensing events, $`f_\mathrm{p}(D_{\mathrm{os}})`$. As a result, the pixel lensing optical depth will be smaller than the value determined by the classic lensing survey, $`\tau _\mathrm{c}`$. We, therefore, apply $`\tau _\mathrm{p}`$ decreased from $`\tau _\mathrm{c}`$ by a factor

| $`{\displaystyle \frac{\tau _\mathrm{p}}{\tau _\mathrm{c}}}=`$ | $`{\displaystyle \frac{[d\tau (D_{\mathrm{os}})/dD_{\mathrm{os}}]f_\mathrm{p}(D_{\mathrm{os}})𝑑D_{\mathrm{os}}}{f_\mathrm{p}(D_{\mathrm{os}})𝑑D_{\mathrm{os}}}}`$ |
| --- | --- |
|  | $`\left[{\displaystyle \frac{[d\tau (D_{\mathrm{os}})/dD_{\mathrm{os}}]f_\mathrm{c}(D_{\mathrm{os}})𝑑D_{\mathrm{os}}}{f_\mathrm{c}(D_{\mathrm{os}})𝑑D_{\mathrm{os}}}}\right]^1.`$ |

(15)
We set the mean value of the classic lensing optical depth towards the BW to be $`\tau _\mathrm{c}=2.4\times 10^6`$ by adopting the determination of Alcock et al. (2000). In the lower right panel of Fig. 2, we present the distribution $`d\mathrm{\Gamma }/d\mathrm{\Omega }`$ for the events detectable under the criteria of $`(S/N)_{\mathrm{th}}=10`$ and $`t_{\mathrm{dur},\mathrm{th}}=6`$ hr. We assume that the total observation time per year is $`t_{\mathrm{obs}}=150`$ days and the total angular area of the monitored field is $`\mathrm{\Omega }_{\mathrm{tot}}=50\mathrm{deg}^2`$.

## 3 Results

In Table 1, we list the determined Galactic pixel lensing event rates expected under various criteria, which are characterized by $`(S/N)_{\mathrm{th}}`$ and $`t_{\mathrm{dur},\mathrm{th}}`$. The important results from the determined detection rates of pixel lensing events and their characteristics are summarized as follows.

First, the total event rate varies significantly in the range $`6\mathrm{yr}^1`$$`120\mathrm{yr}^1`$ depending strongly on the imposed detection criteria. As expected, the event rate increases as $`(S/N)_{\mathrm{th}}`$ and $`t_{\mathrm{dur},\mathrm{th}}`$ decreases. The strong dependency of the detection rate on $`t_{\mathrm{dur},\mathrm{th}}`$ implies that majority of events will have short durations. Therefore, to maximize the event rate it will be essential to identify pixel lensing events by diligently inspecting light variations with signals as low as possible and to conduct followup observations for the identified events as promptly as possible. If events are identified for all signals with $`S/N7`$ and followup observations are carried out within an hour (corresponding to $`t_{\mathrm{dur},\mathrm{th}}=2`$ hr) after the identification, we find that the event rate can reach up to $`\mathrm{\Gamma }_{\mathrm{tot}}120\mathrm{yr}^1`$.

Second, the majority of pixel lensing events will have high amplifications. For example, we estimate that $`90\%`$ of the pixel lensing events detectable under the criteria of $`(S/N)_{\mathrm{th}}=10`$ and $`t_{\mathrm{dur},\mathrm{th}}=6`$ hr will have amplifications $`A20`$ and $`40\%`$ will be EMEs. These fractions increase as the imposed $`t_{\mathrm{dur},\mathrm{th}}`$ becomes shorter. As will be discussed in § 4, events with $`A20`$ are excellent followup monitoring targets for high efficiency extra-solar planet detections (Griest & Safizadeh, 1998; Han & Kim, 2001). In addition, followup observation of EMEs enables one to uniquely determine the masses, distances, and transverse speeds of the individual lenses (Gould, 1997; Sumi & Honma, 2000).

Third, compared to classic lensing events, the events detectable by the proposed pixel lensing experiment will generally involve brighter source stars, as predicted by Gould & DePoy (1998). Under the detection criteria of $`(S/N)_{\mathrm{th}}=10`$ and $`t_{\mathrm{dur},\mathrm{th}}=6`$ hr, for example, we find that the baseline brightness distribution of the pixel lensing event source stars peaks at $`I17`$ (see Fig. 3), which is $`2`$ mag brighter than that of the classic lensing events \[cf. Fig. 7 of Alcock et al. (2000)\]. The source star brightness becomes even brighter as one applies stricter detection criteria, i.e. higher $`(S/N)_{\mathrm{th}}`$ and longer $`t_{\mathrm{dur},\mathrm{th}}`$. Since pixel lensing events tend to have high amplifications, the source during its peak amplification will become even brighter. In Figure 5, we present the distribution of peak brightnesses of source stars for events expected to be detected under the conditions of $`(S/N)_{\mathrm{th}}=10`$ and $`t_{\mathrm{dur},\mathrm{th}}=6`$ hr. One finds that the peak brightnesses of nearly all events are brighter than $`I14`$.

## 4 Discussion

Due to the difference in the characteristics of pixel lensing events from classic ones, the proposed pixel lensing experiment has both advantages and disadvantages. The greatest disadvantage is that the pixel lensing survey should be carried out in parallel with followup observations to obtain meaningful results. However, several groups (Albrow, et al., 1998; Rhie et al., 1999) are already conducting such followup observations, and thus the pixel lensing observations can be made in parallel with these existing observations. Once this problem is resolved, the proposed pixel lensing survey will be able to provide the following important scientific rewards.

First, one can determine the lensing parameters of events with improved precision and accuracy. Precise determination of the lensing parameters is achieved because high precision photometry is possible for pixel lensing events because they generally involve brighter source stars and have higher amplifications. The accuracy, on the other hand, is improved because the effect of blended light from unresolved nearby faint stars is minimized with the increased signal from the lensed source. With the reduced uncertainties in the determined lensing parameters, one can better constrain the nature of Galactic lenses.

Second, one can obtain various useful extra information both about the lens and source star by efficiently detecting lensing light curve anomalies. A lensing light curve deviates from the standard form due to various reasons, e.g. the extended source effect (Schneider, Ehlers, & Falco, 1992; Witt & Mao, 2000), parallax effect (Refsdal, 1966; Gould, 1994), binary lens effect (Schneider & Weiss, 1986; Witt, 1990), planetary lens effect (Mao & Paczyński, 1991; Gould & Loeb, 1992; Bolatto & Falco, 1994), and binary source effect (Griest & Hu, 1992; Han & Gould, 1997). However, the amount of these deviations is usually very small, and thus detecting them requires high-precision observations. Since better photometry is possible for pixel lensing events, the efficiency of detecting light curve anomalies will be increased.

Third, a significant fraction of high amplification pixel lensing events will be ideal monitoring targets of high-efficiency planet detections. Griest & Safizadeh (1998) pointed out that by intensively monitoring high amplification events, planets located in the lensing zone ($`0.6\theta _\mathrm{E}`$$`1.6\theta _\mathrm{E}`$) can be detected with a nearly 100% efficiency. Since majority of Galactic pixel lensing events will have high amplifications, they will be ideal targets for intensive followup monitoring for high-efficiency planet detections. Covone et al. (2000) pointed out that monitoring the pixel lensing events detectable towards the M31 (Crotts & Tomaney, 1996; Ansari et al., 1997; Crotts et al., 2000) will also provide a channel for high efficiency planets detections via microlensing. We note, however, that since followup observations cannot be done for M31 pixel lensing events, although one might detect planet-induced anomalies in the light curves, it will be very difficult to determine the planet parameters (i.e. the mass ratio and separation) just based on the peak part of the light curve. By contrast, since one can perform high-precision followup observations for Galactic bulge pixel lensing events, one can accurately determine the planet parameters.

Another importance of the pixel lensing survey is that one can uniquely determine the physical parameters of individual lenses for EMEs, which will comprise a significant fraction of all pixel lensing events. Gould (1997) showed that followup observations of an EME enable one to measure both the lens proper motion and parallax that together would yield individual mass, distance, and transverse-speed of the lensing object. The proper motion is determined by analyzing the light curve anomalies induced by the extended source effect which occur when the lens approaches very close to the source star. The parallax is determined by observing the difference in the peak parts of the light curves as seen from two Earth observers separated by $`R_{}`$.

## 5 Conclusion

We estimate the event rate of the pixel lensing survey towards the Galactic bulge proposed by Gould & DePoy (1998) expected under various detection criteria and investigate the characteristics of the detectable events. From this investigation, we find that the event rate varies significantly in the range $`6\mathrm{yr}^1`$$`120\mathrm{yr}^1`$ depending strongly on the imposed detection criteria. This implies that to maximize event detections it will be essential to identify events by diligently inspecting light variations and to promptly conduct followup observations for the identified events. Compared to classic lensing events, the events detectable by the proposed pixel lensing survey will generally involve brighter source stars and have higher amplifications. For the pixel lensing events detectable under the criteria of $`(S/N)_{\mathrm{th}}=10`$ and $`t_{\mathrm{dur},\mathrm{th}}=6`$ hr, for example, we find that the baseline source brightness will be $`2`$ mag brighter than those of classic lensing events and most of them ($`90\%`$) will have high amplifications with $`A20`$ and a significant fraction of them ($`40\%`$) will be EMEs. The high amplification events will be especially important because they are excellent targets for effcient extra-solar planet detections. In addition, one can uniquely determine the mass, distance, and transverse speed of individual lensing objects for EMEs.

We would like to thank A. Gould for carefully reading the manuscript and making useful comments. We also would like to express our gratitude to the anonymous referee for making useful comments that contributed to significantly improve the paper. This work was supported by a grant (2000-015-DP0449) from the Korea Research Foundation (KRF).


text:
# Epitaxial MgB2 superconducting thin films with a transition temperature of 39 Kelvin

\[

\]

The recent discovery of the binary metallic MgB2 superconductor  having a remarkably high transition temperature (Tc) of 39 K has attracted great scientific interest . With its metallic charge carrier density and the strongly linked nature of the inter-grains in a polycrystalline form , this material is expected to be a very promising candidate for superconducting device  as well as large-scale applications. However, the fabrication of thin films of this material has not been reported yet. Here, we report the growth of high-quality epitaxial MgB2 thin films on (1102) Al2O3 and (001) SrTiO3 substrates by using a pulsed laser deposition technique. The thin films grown on Al2O3 substrates show a Tc of 39 K with a sharp transition width of $``$ 0.7 K. The critical current density in zero field is $``$ 8 x 106 A/cm2 at 5 K, suggesting that this compound has great potential for industrial applications. Also, Al2O3 is a very promising substrate in terms of superconducting device applications because of its high thermal conductivity and small dielectric constant. For the films deposited on Al2O3, X-ray diffraction patterns indicate a c-axis-oriented crystal structure perpendicular to the substrate surface whereas the films deposited on SrTiO3 show {101} plane as the preferred orientation.

The growth technique used in this study was reported elsewhere. MgB2 thin films were deposited on Al2O3 (AO) and SrTiO3 (STO) substrates by using pulsed laser deposition. The laser energy density was 20 - 30 J/cm2 at a laser flux of 600 mJ/pulse. The resistivity measurements were carried out using the dc four-probe method. The dc magnetic properties were measured with a Quantum Design MPMS superconducting quantum interference device magnetometer. The structures were analyzed using a x-ray diffractometer (XRD).

The typical temperature dependence of the resistivity of MgB2 grown on AO is shown in Fig. 1. The inset is a magnified view near the Tc region. The resistivity of the thin film begins to enter the superconducting transition at 39 K and goes to zero resistance at 37.6 K. A very sharp transition, with a width of $``$ 0.7 K from 90% to 10% of the normal state resistivity, is evident in the inset. This is comparable to most reported values for high-quality bulk samples . The normal-state resistivity at 290 K was $``$ 4.7 $`\mu `$cm, indicating an intermetallic nature with a relatively high charge carrier density . This resistivity is smaller than those for polycrystalline MgB2 wire  and for bulk samples synthesized under high pressure . Most of our films fabricated under the same conditions showed a similar superconducting transition around 39 K.

Fig. 2 shows the zero-field-cooled (ZFC) and the field-cooled (FC) dc magnetization (M) curves of a MgB2 thin film in a 10 Oe field applied parallel to the c-axis. The irreversibility temperature detected at 37.5 K coincides with the zero-resistance temperature obtained from resistivity measurement. The ZFC curve shows a rather broad diamagnetic transition compared to the resistivity data. To estimate the critical current density (Jc), we measured the M-H loop for the same sample as a function of temperature, as shown in the inset of Fig. 2. At 5 K, the Jc calculated using the Bean model was $``$ 8 x 106 A/cm2 at zero field and $``$106 A/cm2 at 1 T, which are 20 times higher than the values obtained for a MgB2 wire  by using transport measurements and slightly smaller than the values for Hg-based superconducting thin films . It should be noted that we used the sample size rather than the grain size in calculating Jc.

Structural analysis was carried out by XRD, as shown in Fig. 3. The a- and the c-axis lattice constants determined from the (101) and the (00l) peaks are observed to be 0.310 and 0.352, respectively. Interestingly, the XRD patterns indicate that the films deposited on (a) AO are epitaxially aligned with the c axis while the films on (b) STO is well aligned with the (101) direction normal to the substrate planes. These results suggest that we may be able to control the orientation of MgB2 thin film simply by using different substrates if the optimum growth condition is explored.

We also find that the AO substrates are chemically very stable during the heat treatment at high temperatures in Mg vapor. For the STO substrates, however, the chemical reaction between the substrate and Mg was observed and appeared as minor impurity phases, as shown in Fig. 3 (a). Considering that STO has a cubic structure with a lattice constant of 0.389 nm and R-plane AO has a0 = 0.476 nm and b$`{}_{}{}^{}{}_{0}{}^{}`$ = (c$`{}_{}{}^{2}{}_{0}{}^{}`$ \+ 3a$`{}_{}{}^{2}{}_{0}{}^{}`$)1/2 = 1.538 nm, it is quite striking that MgB2 thin films grow with preferred orientations on AO and STO substrates even though the lattice-matching relationship between the MgB2 and the substrates is not well satisfied. Further microscopic investigation of the MgB2 crystallographic relations on AO and STO substrates is required.

###### Acknowledgements.

This work is supported by the Ministry of Science and Technology of Korea through the Creative Research Initiative Program.