## Abstract

The ability to accurately and efficiently characterize multiple scattering of waves of different nature attracts substantial interest in physics. The advent of photonic crystals has created additional impetus in this direction. An efficient approach in the study of multiple scattering originates from the Rayleigh method, which often requires the summation of conditionally converging series. Here summation formulae have been derived for conditionally convergent Schlömilch type series *Z*_{n}(*z*) stands for any of the following cylindrical functions of integer order: Bessel functions *J*_{n}(*z*), Neumann functions *Y* _{n}(*z*) or Hankel functions of the first kind

## 1. Introduction

The primary purpose of this paper is to present summation formulae for conditionally convergent Schlömilch type series
*Z*_{n}(*z*) stands for the following cylindrical functions of integer order *n*: Bessel functions *J*_{n}(*z*), Neumann functions *Y* _{n}(*z*) or Hankel functions of the first kind *x* and *D* are real numbers such that 0<|*x*|<*D* and *J*_{n}(*z*) of an even order *n* is given in [2]. Summation equations for general types of series (1.1) have not been derived until now.

Meanwhile, the summation of series (1.1) has both considerable theoretical interest and important applications in physics. For example, the series (1.1) arises naturally in two-dimensional scattering problems of plane waves on diffraction gratings composed of multiple inclusions per unit cell [3–5]. In this context, *θ*_{0} is the angle of a plane wave incident on a grating measured from the normal to the grating. The parameter *D* is given by *D*=*kd*=2*πd*/*λ*, while *x*=*kc*_{lq}, where |*c*_{lq}| represents the distance between two inclusions in a grating unit cell. Here *d* is the grating period and *λ* is the radiation wavelength. In physical applications, the Schlömilch series are often called the lattice sums and they characterize the contribution to the wave's scattering amplitude at the vicinity of an inclusion *l* from an inclusion *q* and all its periodic replicas.

Recently, we have witnessed the emergence of the field of photonic [6,7] and phononic [8] crystals. Such structures can be viewed as a stack of diffraction gratings. Therefore, the ability to effectively and accurately characterize grating structures has become even more important in the broader context of applications in photonic and phononic crystals. A typical approach of modelling such structures requires the consideration of a unit cell with multiple inclusions [9,10]. By introducing intricate defects into such consecutive gratings, it is possible to model wave transport through rather complicated interconnected nano-scale networks of waveguides and resonators embedded in such crystals.

An analogous series to (1.1) with *x*=0 appears in plane wave scattering problems on diffraction gratings with a single inclusion per unit cell
*s* means that the term with *s*=0 is absent from the series. The summation equation of series (1.2) at normal incidence *θ*_{0}=0 was first given by von Ignatovsky [11,12]. By using intricate algebraic transformations von Ignatovsky expressed the series in terms of an absolutely converging series of elementary functions and a finite sum consisting of Bernoulli numbers *B*_{m}. Von Ignatovsky's result has been generalized for off-axis incidence *θ*_{0}≠0 in an elegant approach developed by Twersky [13]. It is shown there that, for the off-axis case, the series (1.2) can be expressed as a finite sum of Bernoulli polynomials *B*_{m}(*x*) and an absolutely converging series of elementary functions. Twersky used Euler's summation formula in the form given by Nörlund [14] to deduce his final results [13]. Summation formulae of the zero order series (*n*=0) (1.2) were included in major reference books for a while [2,15]; the summation formulae for more general series (1.2) are yet to be included.

In [3], Graf's addition theorem [16] has been used to convert the initial series (1.1) into the form
*m*≤100 terms are required to achieve convergence in (1.3). Therefore, use of (1.3) requires calculation of *n*≈100. While equation (1.3) is quite efficient for the summation of the series (1.1) for *D*. The deficiency associated with equation (1.3) is that the sum in (1.3) contained the product of increasing *J*_{m}(|*x*|) factors, which leads to the loss of accuracy for larger *D*.

The ability to model structures with larger unit cells which can accommodate more inclusions has broad and important applications. These applications include, for example, the characterization of wave transport through complex waveguide structures embedded in photonic crystals, consideration of the effects of manufacturing imperfections on operational properties of devices embedded in such crystals or the study of the fundamental problem of Anderson localization of waves in disordered media [17–19]. As it was once well stated—*more is different* [20].

Therefore, it is important and interesting both analytically and from an application perspective to derive Ignatovsky–Twersky like formulae that can be used to directly find values of *D*≫100. In this paper, we adopt the approach developed by Twersky [13] and derive such summation formulae for the series (1.1). It is shown here that the Schlömilch series (1.1) involving Hankel functions or Neumann functions can be expressed as a finite sum of Lerch transcendent functions *Φ*(*z*,*s*,*v*) and an absolutely converging elementary series. For the particular case of normal incidence *θ*_{0}=0, the Schlömilch series (1.1) of Hankel functions or Neumann functions can be expressed as a finite sum involving polylogarithm functions *Li*_{s}(*z*) and an absolutely converging elementary series. The summation formulae for the Schlömilch series (1.1) involving Bessel functions have also been derived and are expressed in terms of a finite sum of elementary functions only. For the reported result in [2] for the even order series,

First, in §2 the key relations are derived, which are used subsequently in §3 to find the summation formulae for the series

## 2. Derivation of the basic relations

In this section, the basic relations are derived, which are used subsequently to transform series (1.1) into a form which allows efficient and accurate summation of the series *π*. Therefore, it is sufficient to derive the summation formulae for only non-negative orders *n*≥0 and the negative orders *n*<0 can be deduced using equation (2.1). Relation (2.1) also holds for the series (1.2) for the same reason.

By direct inspection of series (1.1), it can be established that *x* with −*x* in (2.2), it is easy to establish the relation
*x*, while the values for *x* to be positive.

First, we consider the series *x*- and *y*-axes. Therefore,
*y*>0 and using Sommerfeld integral representation of the Hankel function, we can write
*t*|>1. After substitution of (2.7) into (2.4) and using either the Poisson summation formula
*y*>0. Here the angles *θ*_{p} given by the relation

represent the diffraction orders of plane waves in grating terminology. The representation of (2.4) for *y*<0 is obtained from (2.10) by replacing *θ*_{p} by *π*−*θ*_{p} and it takes the form
*p*_{−}≤*p*≤*p*_{+}, where *m*_{±}. The rest of the orders *p*≥*p*_{+}+1>*m*_{+} and *p*≤−*p*_{−}−1<−*m*_{−} are evanescent orders satisfying *p*≥*p*_{+}+1 we have *θ*_{p}=*π*/2−*i*|*η*_{p}| and for *p*≤−*p*_{−}−1 we have *θ*_{p}=−*π*/2+*i*|*η*_{p}|, where *m*_{±}=0,1,2,… correspond to the Rayleigh wavelength values in which the diffracted wave propagates along the surface of a grating.

For *n*=0, both of the relations (2.10) and (2.12) reduce to
*y*=0, which coincides with the initial series (1.1) for zero order *n* is positive, *n*>0 unless stated otherwise.

## 3. Summation formula for S n H ( x )

In this section, we derive the summation formulae for the series

### (a) The derivation of the key relation

We use the integral representation of Hankel's function

Next, we split the first term in (3.5) into two terms
*D*/2*π*. Then we add and subtract the term _{+} and introduce the *F*_{n}(*x*) term
*ε*=0 in these terms and split the integrals over the propagating orders and evanescent orders and, after changing the variable of integration to *θ*_{p}, we deduce
*θ*_{p}=∓*π*/2±*it*. This leads to the following representation of (3.11)
*n* separately.

For even *n*, we combine the last two integrals over the infinite region together and the first two integrals over a finite region together in (3.12) then expand the factor *e*^{−inθp} into its trigonometric parts, which leads to
*H*^{(1)}(|*x*|) in (3.13) can be dropped, given it is equal to unity regardless of the sign of *x*. Now using the Bessel function representation ([2], p. 79)
*J*_{2n}(*x*) of the Hankel

For odd vales of *n*, we combine the integrals over the infinite region together in (3.12) and we combine the first two integrals over the finite region together. Then we expand the factor *e*^{−i(2n+1)θp} into its trigonometric parts transforming (3.12) into
*J*_{2n+1}(*x*) ([2], p. 79)
*J*_{2n+1}(|*x*|) of the Hankel function cancels given that integral (3.17) always has the opposite sign to the last term in (3.16).

It appears that the integrals over the infinite intervals in (3.15) and (3.18) were first considered by Coates in [21] and these integrals are also subsequently mentioned by Watson for more general cases of non-integer order *n* ([1], p. 313). It turns out that the Coates integrals appearing in equations (3.15) and (3.18) and integrals over the finite region in (3.15) and (3.18) can be calculated in closed form. The details of this calculation are given in appendix A and in the electronic supplementary material. In the next section, we provide details for the calculation of the *F*_{n}(*x*) term using the limiting process based on the Nörlund summation formula [14].

### (b) Calculation of the *F*_{n}(*x*) term

The first two series in (3.10) converge without requiring the application of a convergence limiting process. By contrast, the *F*_{n}(*x*) term is calculated by using the Nörlund summation formula [14] expressed in the form
*B*_{n}(*x*) is a Bernoulli polynomial. In the subsequent calculations, we use the following form of relation (3.19)
*e*^{−εz} in (3.19), given *p*. We consider even and odd *n* cases separately. First, we derive the expression for *F*_{2n}(*x*).

#### (i) Derivation of the expression for *F*_{2n}(*x*)

For even values of order *n*, *F*_{2n}(*x*) takes the form
*Φ*(*z*,*s*,*α*) defined by [2,29]
*t*|<2*π*, equation (3.24) can be written as
*x*|/*D*<1. The final expression for *F*_{2n}(*x*) can be deduced from (3.27) after straightforward algebraic transformations
*x* cannot take zero value in the initial sum *F*_{2n}(*x*) can be calculated for *x*=0 . Using the limit [2]

where

#### (ii) Derivation of the expression for *F*_{2n+1}(*x*)

The calculation for the odd order terms, *F*_{2n+1}(*x*), is quite similar to the even order case. For odd values of *n* relation (3.9) takes the form
*n* case. The final expression for *F*_{2n+1}(*x*) takes the form
*n* case, expression (3.37) reduces to the corresponding term of Twersky's result [13] for *x*=0.

#### (iii) Final forms of the series S n H ( x )

The final expression for the series *n*, we substitute expression (3.28) for *F*_{2n}(*x*) and the values for integrals (A 4) and (A 9) into (3.15). Some pleasing cancellations occur. The double sums of integral (A 4) cancel the double sums of the Coates integral (A 9), while the Neumann function of the Coates integral (A 9) *Y* _{2n}(|*x*|) cancels the corresponding term in (3.15). The second term of *F*_{2n}(*x*) (3.28) cancels the last term of integral (A 4) and the expression for *F*_{2n+1}(*x*) (3.37) and the corresponding integrals (A 8) and (A 10) into (3.18). The double sums of integral (A 8) cancel the double sums of the Coates integral (A 10) and the Neumann function of the closed-form expression of the Coates integral (A 10) cancels the last term of (3.18). The first two terms of *F*_{2n+1}(*x*) (3.37) cancel the single sum of integral (A 8) and the expression for *θ*_{0}≠0. Hence the series (3.39) can be rewritten as
*x*=*D*/2, the Lerch transcendent functions in equations (3.38) and (3.41) reduce to Euler polynomials [2] given the connection *E*_{m}(*x*)=2*Φ*(−1,−*m*,*x*). Therefore, for this case also the Schlömilch series involving Hankel functions are represented in terms of elementary functions only. Given the rich properties of the Lerch transcendent function there are other interesting argument values where

It is worthwhile mentioning that the Lerch transcendent function in (3.38) and (3.41) is also called the Lerch zeta function, denoted by *L*(*x*/*D*,*s*,*α*) provided the following relation *Φ*(*e*^{i2πx/D},*s*,*α*)=*L*(*x*/*D*,*s*,*α*). In §6, the accuracy, efficiency and the numerical verification of the derived summation formulae (3.38) and (3.41) will be demonstrated.

## 4. The summation formulae for the S n J ( x ) and S n Y ( x ) series

In this section, we provide details of the derivation of the summation formulae for the

The Bessel series for *t*|<1. Now we substitute (4.2) into (4.1) then interchange the integration and summation, and, using either the Poisson summation formula (2.8) or the identity (2.9), the expression for *n* order cases. For the even case, we deduce
*T*_{n}(*x*) using the following connection

The expressions for the *n* is replaced by zero throughout.

## 5. Series S n Z ( x ) for *θ*_{0}=0

In this section, we consider the Schlömilch series (1.1) for *θ*_{0}=0. These series arise in the case of a plane wave incidence normal to the surface of a grating. As the incidence angle vanishes *θ*_{0}=0, for the normal incidence the initial Schlömilch series takes the form

The formulae for the normal incidence can be readily deduced from the corresponding general formulae for off-axis incidence by taking the limit *p*_{−}=*p*_{+}=[Δ]=[*m*_{±}], where *p*_{±} are the largest integer smaller than *m*_{±} . For the case of the Schlömilch series involving Hankel functions *Li*_{s}(*z*) [31]
*Li*_{s}(*z*) is defined by

Relation (5.5) can be used to derive the summation equation for the series *n*=0 in (5.5) and dropping the last term to obtain
*θ*_{0}=0 involving Bessel functions *θ*_{0}=0 involving the Neumann series take the forms

## 6. Numerical verification

Apart from the considerable theoretical interest of the presented results, the obtained expressions for the Schlömilch series are important for applications also. In this section, we provide details of the numerical verifications for the derived analytical expressions (3.38) and (3.41) for the Schlömilch series

In applications, the numerical values for *D*≫1. At the same time, larger values of *D* are required either for the modelling of wave propagation in photonic or phononic crystals with embedded intricate structures [9,10] or for investigating the fundamental phenomenon of Anderson localization of classical waves [17–19].

The analytical expressions (3.38) and (3.41) are ideally suited for numerical calculations. Indeed, the obtained expressions involve the absolutely converging series of elementary functions and Lerch transcendent functions which can be calculated with high precision. The rate of convergence of the elementary series *p*^{−(2n+1)}, where *p* is the summation index of the series. Therefore, the slowest convergence rate for the even order series occurs for *p*^{−3}. The rate of convergence for odd order series *p*^{−(2n+2)}. So the slowest convergence rate for odd order series occurs for *p*^{−2}. If required, the convergence rate can be accelerated further using Kummer's transformation [16,27,28]. Note that the convergence rate increases substantially with the increase in the series order *n*.

First, the consistency between the obtained equations (3.38), (3.41) and the relation (1.3) has been verified. In table 1, we present the comparison data for the values *D*=80, *θ*_{0}=*π*/6 and *x*=79. There is an excellent agreement between the computed values. The values calculated based on (3.38) and (3.41) are correct with all presented figures. The values based on the relation (1.3) have been calculated with 160 terms, which requires calculation of series (1.2) for up to 80 orders of *N*≈30 inclusions per unit cell for typical wavelengths. The obtained expressions for the Schlömilch series (3.38) and (3.41) enable the modelling of structures with substantially larger unit cells, thus allowing a substantial increase in the number of inclusions per unit cell. In figure 1, we plot the partial sums of the series *D*=1000, *θ*_{0}=*π*/6 and *x*=250. The oscillating curves in figure 1*a*,*b* are the real and imaginary parts of the partial sums of the series

The summation equations for the series *θ*_{0}=0. The plots of the partial sums of *θ*_{0}=*π*/6 and they are not provided here. The obtained summation formulae (3.38) and (3.41) satisfy analytically the quasi-periodicity conditions (2.2) and (2.3). These conditions have also been verified numerically.

Given the obtained series (3.38) and (3.41) converge absolutely, the number of significant figures that can be achieved is only limited by the standard function accuracy and limits set by a computer. The presented calculations have been done using Wolfram Mathematica [32] with double precision accuracy. The computational time is of the same order as for the calculation of the series (1.2) based on the Twersky formulae [13].

Having verified numerically the accuracy and efficiency of the obtained summation formulae, it is interesting to plot *x*. In figure 5, we plot *x* for *θ*_{0}=*π*/6, *D*=80 in the domain −*D*<*x*<*D*. The real part of *a*, whereas the imaginary part is plotted in figure 5*b*. Both the real and imaginary parts of *x*. The real part of *x*=0, whereas the imaginary part diverges. Both the real and imaginary parts of *x*=±*D*.

In figure 6, we plot the same relationship as in figure 5 but for the series *D* increases the oscillatory behaviour of *x* and we do not provide their plots here.

## 7. Conclusion

Accurate and efficient summation equations for the Schlömilch type series (1.1) have been derived. This is achieved by using the Euler summation formula in the form given by Nörlund [14]. The obtained analytical formulae (3.38) and (3.41) for the Schlömilch series involving the Hankel functions or Neumann functions are expressed in terms of an absolutely convergent series of elementary functions and a finite sum of Lerch transcendent functions. The obtained equations have been verified numerically and their accuracy and efficiency have been demonstrated. The summation equations for the Schlömilch series involving the Bessel functions are expressed in terms of a finite sum of elementary functions only (equations (4.4) and (4.5)).

The summation equations for the Schlömilch series at the normal incidence *θ*_{0}=0 have also been derived (equations (5.10)–(5.13)). For the Schlömilch series involving the Hankel functions or Neumann functions, the summation formulae are expressed in terms of an absolutely convergent series of elementary functions and a finite sum of polylogarithm functions. The closed-form expressions for the Coates integrals [21], required for the derivation of the summation equations, have been obtained (see equations (A 9), (A 10) and the electronic supplementary material). These analytical expressions for the Coates integrals have been verified numerically.

The obtained results, apart from their considerable analytical importance, also allow the modelling of scattering problems with substantially larger unit cells. Therefore, the presented results permit accurate and efficient modelling of extended devices embedded in photonic crystals, or characterization of the transmission properties of phononic crystals. The effects of manufacturing imperfections on their functional properties can also be investigated. The obtained results are also important for the accurate characterization of the phenomenon of Anderson localization of classical waves [17–19], where the ability to consider substantially larger unit cells is essential.

Therefore, the derived analytical formulae presented here complement the results by von Ignatovsky [11,12] and Twersky [13] and extend the available analytical toolbox for modelling a variety of important and exciting fundamental and applied problems in physics.

## Data accessibility

The data used here have been produced using the obtained equations.

## Competing interests

I have no competing interests.

## Funding

I received no funding for this study.

## Acknowledgements

The computational resources used in this work were provided by Intersect Ltd. The unequivocal support and encouragement from my family and the full understanding of my wife (most of the time) is gratefully acknowledged.

## Appendix A. Closed form of integrals over the finite interval in equations (3.15) and (3.18)

Here we provide details for the calculation of integrals over the finite region appearing in (3.15) and (3.18). To find the closed-form expression for the integral
*x* and *n*.

The closed-form expressions for the Coates's integrals are given by the relations

- Received May 30, 2015.
- Accepted October 22, 2015.

- © 2015 The Author(s)