{'fraction_non_alphanumeric': 0.05099859728809031, 'fraction_numerical': 0.019003406586066394, 'mean_word_length': 4.113217213114754, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Antisymmetric Wilson loops in N = 4 SYM: from exact results to non-planar corrections

We consider the vacuum expectation values of 1/2-BPS circular Wilson loops in N=4 super Yang-Mills theory in the totally antisymmetric representation of the gauge group U(N) or SU(N). Localization and matrix model techniques provide exact, but rather formal, expressions for these expectation values. In this paper we show how to extract the leading and sub-leading behavior in a 1/N expansion with fixed 't Hooft coupling starting from these exact results. This is done by exploiting the relation between the generating function of antisymmetric Wilson loops and a finite-dimensional quantum system known as the truncated harmonic oscillator. Sum and integral representations for the 1/N terms are provided.


Introduction
Wilson loop operators have played a central role in the development of gauge/gravity dualities [1,2].
A key ingredient in the description of Wilson loop operators is the representation of the gauge group, typical ones for U (N ) and SU (N ) being the fundamental, totally symmetric and totally antisymmetric representations. Depending on the rank of the representation, the holographic dual corresponds to probe strings or D-branes [27,28] propagating on AdS 5 × S 5 , or to fully backreacted bubbling geometries [29][30][31][32]. In all cases, the on-shell action of the gravitational object agrees perfectly with the matrix model calculation [33] at leading order in 1/N and 1/ √ λ. We expect, however, that a thorough examination should yield a match at the next-to-leading order as well. It is then crucial to systematically extract these corrections on both sides of the duality.
From the gravitational perspective the calculation of sub-leading corrections amounts to analyzing semi-classical fluctuations of the background configurations and computing the corresponding one-loop partition functions. There has been considerable effort in this direction [34][35][36][37][38][39][40][41][42][43][44], although the calculations are plagued with ambiguities inherent to string theory in curved spaces.
On the gauge theory side, solving the matrix model is a non-trivial task, and computing subleading corrections is a conceptually clear, albeit technically difficult procedure. Recently, subleading corrections in 1/N to the Wilson loop in the totally antisymmetric representation have been found using loop equation techniques [45] and topological recursion in the Gausssian matrix model [46] (see also [47] for the case of symmetric representations). A particularly interesting development was reported in [48], where the authors managed to compute the exact vacuum expectation value of the circular Wilson loop in arbitrary irreducible representations using the method of orthogonal polynomials. In other words, they solved the matrix model exactly. While useful for some purposes, their results are rather formal, and it is not at all obvious how to extract the large-N (and large-λ) limit from them, much less any sub-leading corrections, except for simple cases. So far, they have only been used for numerical comparison with other approaches.
In this paper, we address and solve the problem of extracting the leading and first sub-leading terms in the 1/N expansion of the generating function of totally antisymmetric Wilson loops starting from the exact results in [48]. Our motivation is two-fold. First, it is obviously interesting to see how it can be done. Second, the techniques used in [45,46] solve the matrix model in the continuum limit, in which the discrete poles of the matrix model resolvent give rise to a cut singularity and a continuous eigenvalue density. In other words, the analyticity properties of the resolvent change in this limit, and it would be interesting to see whether and where this has any implications. We anticipate that we do not find any at order 1/N .
The paper is organized as follows. In section 2, we review the Gaussian matrix model and reproduce the formal result of [48] for the generating function of antisymmetric Wilson loops. At the same time, we establish a connection with the algebra of the truncated harmonic oscillator, which will be crucial for our developments. In section 3, we extract the leading large-N behaviour (fixed λ) of the generating function from the formal solution of the matrix model. The most substantial part of the paper is section 4, which contains the calculation of the 1/N terms of the results are shown to agree with [45,46], but turn out to be somewhat simpler in their final form.
Section 5 contains the conclusions. Some technical details are deferred to the appendices.

Gaussian Matrix Model
Localization techniques [5] map the expectation value of the circular Wilson loop in N = 4 SYM with gauge group U (N ) to an expectation value in a Gaussian matrix model. We begin this section by reviewing this model and developing some results that will be relevant for what follows.

Partition Function and Expectation Values
The Gaussian matrix model is defined by the partition function where X is a N × N hermitian matrix. If the gauge group is SU (N ) the matrices are also traceless, condition that can be implemented with a Lagrange multiplier. The expectation value of any quantity F (X) in the Hermitian ensemble is then given by When F is invariant under similarity transformations, i.e., F (V XV −1 ) = F (X), one can diagonalize the matrix X in terms of its eigenvalues x n , n = 1, . . . , N , and integrate out the remaining "angular" variables. This yields where the transformation Jacobian is known as the Vandermonde determinant.
The Vandermonde determinant enjoys several properties that we can exploit to our advantage.
For example, one can show that where P n (x) = h n x n + · · · is any family of polynomials that can be chosen to our convenience.
Using these facts, the partition function may be written as Here, y and z are arbitrary constants, and the coupling g is defined by . (2.8) We now introduce a quantity that will play a central role in our analysis of the circular Wilson loop. Let us define the N × N matrix Although this definition might seem spurious at this point, its relevance will become clear momentarily. Writing the determinant (2.6) explicitly as we find that The same procedure can be applied to the matrix model expectation values (2.2), which become (2.12) We emphazise that the integrals are independent of y and z by virtue of (2.5).
It should not come as a surprise that a very convenient choice of polynomials is where He n (x) are Hermite polynomials. 1 The family P n (x) is then orthonormal with respect to the Gaussian weight, namely, (2.14) 1 Even though we are physicists, we work with the probabilists' version of the Hermite polynomials in order to avoid some awkward factors of √ 2. Some basic properties of the Hermite polynomials are listed in appendix A.
Furthermore, the matrix elements (2.9) can be computed explicitly [49], yielding n (x) are Laguerre polynomials. Owing to the properties of these polynomials, the distinction between the cases m ≥ n and m < n is not necessary.

Properties of the matrix I(y, z)
We now proceed to discuss some noteworthy attributes of the matrix (2.15). First, one easily verifies that as can be shown using (2.18), (2.19) and (2.23). This implies that the eigenvalues of I(y, z) come in reciprocal pairs and that its characteristic polynomial is either palindromic or anti-palindromic.
The properties (2.21) are reminiscent of an exponential behavior. This is no coincidence. Indeed, one can verify the remarkable relation where A is the matrix given by A n,n+1 = √ n , (2.26) and all other entries vanishing. More explicitly, Note that any power of A is traceless and that A N = 0. The same is true for A T , of course.

Truncated harmonic oscillator
We promptly notice that A and A T are nothing more than the matrix representation of the ladder operators of the harmonic oscillator truncated to the first N energy eigenstates. Surely, the number is diagonal and satisfies However, A and A T themselves do not fulfill the Heisenberg algebra. Instead, their commutator which is traceless, as must be for finite-dimensional operators.
The N -dimensional quantum mechanical system formed by the operators A and A T , chiefly called the truncated harmonic oscillator, is well known in quantum optics. We recall here some of its properties, following the recent account [50] and refer the interested reader to more references in that paper. In the next sections we will exploit this connection to the truncated harmonic oscillator in order to extract the leading and sub-leading large N behavior of the circular Wilson loops.
In a fairly obvious notation, we denote by |n , n = 0, 1, . . . , N − 1 the column vector with a 1 in the n + 1-th position and zeros elsewhere. These states are eigenstates of the number operator N and whence are called the number basis. According to (2.27), the ladder operators act upon them as Of course, the matrix elements (2.15) are I mn (y, z) = m| e yA T e zA |n . (2.33) Actually, this same expression can be obtained by computing the matrix elements of the displacement operator D(y, z) = e ya † e za = e yz 2 e ya † +za in the full, infinite-dimensional, quantum harmonic oscillator [51], and then truncating to the first N number states.
We will now construct a basis for the N -dimensional Hilbert space which is more suited to our purposes. Consider the states where ζ is an arbitrary real parameter. By virtue of the recursion relations for the Hermite polynomials, the actions of A and A T on |ζ can be formally represented as We see that |ζ is an approximate eigenstate of A + A T . It becomes an exact eigenstate if ζ is a root of He N , This construction makes explicit the fact that A + A T is actually the companion matrix of the Hermite polynomial of degree N , namely, Now, the Hermite polynomial He N (ζ) has precisely N distinct roots, which we denote in increasing order by ζ i , i = 1, 2, . . . , N . One can show that the corresponding states |ζ i are linearly independent by computing their inner product with the aid of the Christoffel-Darboux formula, He n (η) He n (ζ) Then, 40) and the vectors |ζ i form an orthogonal basis of the Hilbert space called the position basis.
The matrix elements of I(y, z) in the position basis are We will use this expression extensively in the next sections. Before we move on, we would like to point out that it is possible to write an exact expression for I −1 (y, z) that takes a particularly simple form. Indeed, since A T acts on |ζ as a derivative, we have Even though we will not use this formula in what follows, we believe that it could lead to some simplifications in the analysis of sub-leading corrections to the expectation value of the circular Wilson loop.

Generating Function
The precise relation between the circular Wilson loop and the Gaussian matrix model is Tr R e X , (2.43) where R denotes the representation of the U (N ) gauge group and Tr R the corresponding trace.
The vacuum expectation value on the left is defined on the N = 4 SYM theory; the right hand side corresponds to an insertion in the matrix model.
In this paper we will be concerned with the totally antisymmetric representation of rank k defined by the Young diagram The structure of the trace Tr A k can get increasingly complicated as the rank k grows. Indeed, for any N × N matrix X we have and so forth. Happily, there is a natural way to encode this structure into the generating function A straightforward calculation using (2.9), (2.10) and (2.12) reveals that the matrix model expectation value of this generating function is where the matrix I(g, g) was introduced in (2.15). The dependence on y and z disappears as a consequence of (2.21) and (2.20). This is essentially the answer reported by Fiol and Torrents in [48], except that they use a slightly different matrix, which we callĨ(g, g), obtained from I(g, g) by multiplying the m-th row by (m − 1)! g −m and dividing the m-th column by the same factor.
The determinant in (2.50) does not change under these operations, so which is the expression given in [48]. Note also that our generating functional differs by a factor of t N from theirs.
In the answer (2.50) we recognize the form of the generating function of antisymmetric traces (2.49) for the matrix I(g, g), which allows us, using (2.49) and (2.43), to obtain the formal yet remarkably simple result (2.53) For the gauge group SU (N ), the matrix model must be restricted to traceless matrices. This can be achieved either by using a Lagrange multiplier or by explicitly isolating the trace component in the matrix model integral, as was done in [45]. The outcome is (2.54) Because the representations A k and A N −k are conjugate to each other for the gauge group SU (N ), . This can be demonstrated by noting that, for an invertible matrix X, (2.49) implies Now, since I(g, g) is similar to its inverse, and det I(g, g) = 1, this gives which proves the assertion. Equation (2.56) also shows that the generating function F A (t; I(g, g)) is a palindromic polynomial.
The main object of study in the remainder of this paper is the function from which the traces Tr A k [I(g, g)] can be calculated by where we have introduced the ratio We are interested in the large N regime, keeping the 't Hooft coupling λ and the ratio κ fixed. In this regime, the integral in (2.58) is dominated by the saddle point value, 2 where t * satisfies the saddle point equation Moreover, F(t) admits an asymptotic expansion in 1/N , Our aim is to calculate the terms F 0 (t) and F 1 (t) from the exact, but formal, expression (2.57).

Wilson loops at leading order
In this section, we shall calculate the leading order term F 0 (t) of the generating function. Of course, F 0 (t) is known both from a matrix model calculation and from the holographic dual [28,33] (the latter implying large λ in addition to large N ). Here, we reproduce it starting from the exact expression (2.57) by means of two different calculations, which exploit the relation of the matrix I to the truncated harmonic oscillator.
Let us start in the number basis. This calculation will lead to the expression for F 0 (t) that Okuyama [46] attributes to an unpublished note by Beccaria. Starting from (2.57) and Taylorexpanding the logarithm we may write In order to compute the powers of the matrix I(g, g) we resort to the representation (2.25). To leading order in 1/N , commutators of the matrices A and A T can be neglected, so that we can essentially normal order the product in the trace. This leads to We can now expand the exponentials and calculate the traces. We defer this little calculation to appendix B. With the results (B.3) and (B.5), this yields The leading term in 1/N is found to be where I 1 (x) denotes a modified Bessel function. This is precisely J 0 in (2.21) of [46].
A slightly simpler way of arriving at (3.4) is to recognize (3.2) as The trace follows easily from (2.15) and the identity [52] N m=1 Hence, (3.5) becomes For large order, the Laguerre polynomials satisfy [52] lim with I α (x) being a modified Bessel function. Recalling that g = λ 4N , applying the limit (3.8) to (3.7) gives (3.4).
For a comparison with the matrix model calculation in the saddle point approximation [33], one can use the following integral representation of the modified Bessel function [52], After substituting (3.9) into (3.4), the summation can be carried out to give (3.10) A simple change of variables transforms this integral into the one in (2.26) of [33].
The computation of F 0 (t) turns out to be considerably simpler in the position basis. The reason for this is that, to leading order in 1/N , we can also write where all commutators in the Baker-Campbell-Hausdorff formula contribute to the O(1/N ) terms.
According to (2.37), the operator in the exponential is diagonal in the position basis, so that (2.57) simply becomes The sum is taken over the roots of the Hermite polynomial He N (ζ). For large N , the distribution of zeros becomes dense and is described by the Wigner semi-circle law. Using formulas (A.11) and (A.12) to determine the integral measure and dropping the 1/N contributions, we get which agrees with (3.10).

Wilson loops at next-to-leading order
In this section, we extend the calculation of the generating function F(t) to order 1/N , i.e., we obtain the subleading term F 1 (t) in the expansion (2.62). Our results agree with [46], although we also managed to make a slight simplification. In the holographic regime of large λ, we will find agreement with [45].
Exploiting the relation between the matrix I(g, g) in the position basis and the zeros of Hermite polynomials appears to be the most promising path to compute the subleading corrections. Starting again from (3.1), we will proceed as follows. In subsection 4.1, we compute the matrix elements of  The leading term I (0) is known from the calculation at the end of the previous section and is diagonal. The main effort in subsection 4.1 consists in finding I (1) . We then need to compute powers of I(g, g). This is not as straightforward as one might think because I (1) is a dense matrix.
Consider, for example, The last term, which naively seems to be of order 1/N 2 , is, in fact, of order 1/N , because the matrix multiplication contains a sum over N terms of order 1. The same happens for higher powers. The  In the sum that must be subtracted (the light shaded area), we introduce the variable s = Notice that, until now, it was not necessary to require that ζ be a root of He N . From now on, we will assume that it is. This implies that the s = 0 term in the sum in  We shall not investigate this recursion further, because it will turn out that only a s,0 is relevant for our purposes. After substituting (4.8) into (4.7), we can reorder the summations over s and l by introducing r = s − 2l and obtain With these results, let us now consider the matrix element of e gA † e gA between two (nonnormalized) position eigenstates. Remember that the sum over k and l in (4.4) is given by (4.6) minus (4.11). We find ζ i | e gA † e gA |ζ j = e gζ j − 1 2 g 2 ζ i |ζ j (4.12) where we have used the definition of the states |ζ (2.34), extended the summation over k (this can be done because of the binomial) and let s = N − n. Furthermore, rewriting the sums over l and k as a sum over n = l + k and k gives ζ i | e gA † e gA |ζ j = e gζ j − 1 2 g 2 ζ i |ζ j (4.13) He N −s (ζ i ) He N −r (ζ j ) ∞ n=0 g r+s+2n (r + s + 2n)! n k=0 r + s + 2n k a r+2(n−k),n−k (N − s − k)! .
In the first term on the right hand side, one can recognize the leading order result, corrected by the exponential e − 1 2 g 2 . Therefore, the second term is of order 1/N , so that leading order relations can be used to manipulate it. A short inspection of the sum over k in (4.13) for some low values of n shows that the leading order contribution in 1/N comes from the term with k = n. Therefore, using (4.10), we get where subleading terms in 1/N have been omitted. Then, after normalizing the matrix elements (4.13) by means of (2.40) and substituting (4.14), we obtain I ij (g, g) = e gζ j − 1 2 g 2 δ ij (4.15) It is possible to obtain simple expressions for the ratios He N −s (ζ i )/ He N +1 (ζ i ). We include the calculation in Appendix D, where we show that, to leading order in 1/N , Hence, we obtain our final result for the matrix I(g, g) in the position basis,

Powers of I(g, g) and their traces
In this subsection, we return to calculate the quantities needed in the generating function (3.1), namely the powers of the matrix I(g, g) and their traces. As mentioned at the beginning of this section, calculating the powers of I(g, g) to order 1/N is not straightforward. Fortunately, after calculating a few powers, e.g., I 2 (g, g) and I 3 (g, g), one recognizes a pattern that can then be proven by induction for all powers. The calculation is somewhat tedious and involves the facts that the modified Bessel functions have the generating function and that they satisfy the addition theorem [52] I n (x + y) =  3 In the second term, which is already of order 1/N , this care was not needed.
Here, M (m) rs denote some (infinite) matrices, which satisfy the recursion relations and M (1) rs = 0. Let us evaluate the trace of (4.21). This calculation is nearly identical to the one resulting in (C.6), in particular with respect to the cancellation of various 1/N contributions, and yields The first term on the right hand side reproduces the leading order expression (3.4), with an exponential correction that can be traced back to the exponential in (2.52), which is missing in (2.57).
The second term, which is the substantial part of the 1/N contributions, will be evaluated in the remainder of this subsection.
Henceforth, let z = √ λ in order to simplify the notation. For m = 2, we simply have from where we have re-ordered the summations over r and t. For m > 2, using the same reordering, where S a (v; z) stands for and This yields where we have recognized the summation formula (4.20).
Unfortunately, the same trick does not suffice to easily obtain the nested sums for a > 2.
However, using a generating function for S a (v; z), we prove in appendix E that This is a remarkable result, which we have not found in the literature. With (4.28), we can return to (4.24), which simplifies to

Generating function
Putting together (3.1), (4.23) and (4.29), we obtain the generating function F(t) to order 1/N , which appears in the integrand in J 1 . The step from the second to the third line consists in symmetrizing the summands with respect to a → m − a. We note that our result for the 1/N term is slightly simpler than Okuyama's, because it does not involve an integral.
In the remainder of this section, we will find an integral representation of F 1 (t). 4 Consider the first term in brackets in (4.30). Because it contains F 0 (t), but also corrections in 1/N , we will denote it by F 0 (t). After using the integral representation of the modified Bessel function (3.9), the sum over m can be performed, which yields Now, consider the second term in brackets in (4.30), which we shall denote by F 1 (t)/N . Using the standard integral representation for I 0 and I 1 [52], the sum over a can be done, which gives The last term on the right hand side, which does not depend on φ, integrates to zero in the φ-integral. For the remaining term, performing the sum over m in (4.33) leads to Thus, combining (4.32) with (4.34), we obtain the integral representation for F 1 (t) as

Holographic regime
Our aim in this subsection is to evaluate the saddle point value of (4.35) in the regime of large λ.
We remind the reader that the saddle point, in this regime, is given by where the angle θ * satisfies while the fraction in front of the logarithm is approximately −1. Therefore, the contribution of the region A is In the region A (φ < θ * < θ), we have . (4.43) The term in parentheses precisely cancels the denominator of the term in front of the logarithm, and what remains of the integrand is again exponentially suppressed for large λ. Thus, In the remaining region C, where θ < θ * and φ < θ * , we have Therefore, we can write the contribution from the region C to the integral as Because the integration domain in (4.46) is symmetric with respect to φ and θ, we can symmetrize the integrand and rewrite (4.46) as The first term is readily integrated and yields In the second term in (4.47), we write cos 2 θ − cos 2 φ = sin 2 φ − sin 2 θ and realize that, using the symmetry of the integrand and the integration domain, we can write Here, the φ-integral must be interpreted as the principle value because of the pole for φ = θ, whereas there was no pole with the symmetric integrand. However, after rewriting (4.49) as the φ-integral is easily done using integration by parts, which also takes care of the principal value.

Conclusions
In this paper we tackled the problem of non-planar corrections to the expectation value of We have provided both sum and integral representations of the 1/N terms and have evaluated them explicitly in the holographic large-λ regime, which allows for easier comparison with the holographic dual picture. This term should match with the gravitational backreaction of the D-brane on the gravity side. A particularly interesting observation is the connection between the Wilson loop generating function and the finite-dimensional quantum system known as the truncated harmonic oscillator. This system, which is familiar to the Quantum Optics community, provides a description of the problem that seems to be more amenable to an asymptotic 1/N expansion. En route, we obtained interesting mathematical relations and sum rules involving the Hermite polynomials. It would be interesting to see these formulas, which we proved in appendices D and E, in different applications.
One can envisage two main lines of generalization of the present work. First, it would be interesting to extend the methods developed here to other representations of the gauge group.
A particularly interesting case is the totally symmetric representation S k (see also [47]), whose generating function is slightly more complicated than the antisymmetric one but still quite simple, so the problem seems tractable. Second, one may investigate how the present approach extends to higher orders in 1/N . At first sight, there are a number of technical obstacles that must be overcome, because order-1/N 2 terms have been neglected at many points of the calculation. So, the question whether our approach lends itself to a systematic 1/N expansion is highly non-trivial.
This problem is closely related to the fact that the large, but finite-N matrix model solution differs in its analyticity properties from the continuum limit. Moreover, corrections to the saddle point calculation of the Wilson loop expectation values become relevant at order 1/N 2 . We leave these interesting questions for the future.

A Some properties of the Hermite polynomials
In this appendix, we list a number of formulae regarding the (probabilists') Hermite polynomials, which are useful for the analysis in the main text. These relations can be found in standard references [49,52]. Sometimes a translation from the physicists' version of the polynomials is necessary. They are related by The Hermite polynomials He n (x) satisfy the differential equation and the asymptotic expansion of the parabolic cylinder function and u s (t) are polynomials The function g(µ) is irrelevant for our purposes. Setting µ = √ 2N + 1 and t = cos θ, these relations imply that the N roots of He N are given approximately by

B Trace of normal ordered products
In this appendix, we calculate the trace of the normal ordered product A T m A n . First, consider A T n A n . Using the definition of the number operator (2.28) and the commutators (2.29), we get where we have iterated the first two steps to arive at the final expression. where θ i are defined in (A.11) in terms of the zeros of the Hermite polynomial He N (ζ). In the large N limit, it is justified to convert such a sum into an integral. Here we describe here how to do it correctly to order 1/N .
We start by setting x = i − 1 2 and use the Euler-Maclaurin formula in mid-point form [53]. The mid-point form has the advantages that the integration domain lies manifestly symmetric within the inveral (0, N ) and that the boundary terms in the Euler-Maclaurin formula contain only derivatives of the integrand. The latter turn out to contribute at least of order 1/N 2 and are, therefore, irrelevant for our purposes. Hence, we have In the second equality, we have changed the integration variable using (A.11). The angle θ 1/2 , which marks the tiny edges missing from the interval (0, π), also follows from (A.11), To obtain the second term on the right hand side we have assumed that the function f is regular at 0 and π.
Let us check the formula (C.4) by calculating the trace of the matrix I(g, g) in the position basis In the position basis, we have from ( We especially point out the presence of the last term on the second line, which comes from the edge terms of (C.4) and is crucial for cancelling other 1/N contributions. The bracket on the penultimate line vanishes by means of a well-known summation formula of the modified Bessel functions [52].
In these expressions, we have dropped all contributions of order 1/N 2 . Obviously (C.6) agrees with (C.5) to order 1/N . D Proof of (4. 16) In this appendix, we provide a proof of the asymptotic formula where we have dropped the term s of order 1/N in front of U s−3 (cos θ) and used the recursion relation for the Chebychev polynomials [49] in the last step.

(E.5)
We rearrange the sums over v and r by introducing s = v − r and t a = r, such that sin[(s + t a )θ] I ta−t a−1 (z) · · · I t 2 −t 1 (z) I s+t 1 (z) . (E.6) Now, we can extend the summation over t a to all integers without changing the value of the sum.
The reason for this is that if M rs is symmetric, by the antisymmetry of the sine. It is easy to see that this is the case for the matrix involving the rest of the t-summations and the modified Bessel functions in (E.6). Therefore, sin[(s + t a )θ] I ta−t a−1 (z) · · · I t 2 −t 1 (z) I s+t 1 (z) {sin[(s + t a−1 )θ] cos(t a θ) + cos[(s + t a−1 )θ] sin(t a θ)} × I ta (z) I t a−1 −t a−2 (z) · · · I t 2 −t 1 (z) I s+t 1 (z) = e z cos θ ∞ v=1 S a−1 (v; z) sin(vθ) . (E. 8) In the step to the last line, we have carried out the sum over t a using (E.2) and used the fact that the sum of sin(t a θ) I ta (z) vanishes. Continuing recursively, we obtain ∞ v=1 S a (v; z) sin(vθ) = e (a−2)z cos θ ∞ v=1 S 2 (v; z) sin(vθ) = z 2 sin θ e az cos θ , (E.9) knowing that (E.4) is true for a = 2. Thus, we have proven (E.4) for all a ≥ 2, which implies (E.1).