The Trace Method for Cotangent Sums

This paper presents a combinatorial study of sums of integer powers of the cotangent which is a popular theme in classical calculus. Our main tool the realization of cotangent values as eigenvalues of a simple self-adjoint matrix with integer matrix. We use the trace method to draw conclusions about integer values of the sums and an explicit evaluation. It is remarkable that throughout the calculations the combinatorics are governed by the higher tangent and arctangent numbers exclusively. Finally we indicate a new approximation of the values of the Riemann zeta function at even integer arguments.


Introduction
In the present paper is concerned with cotangent sums of the form (1.1) Spm, n, αq " n´1 ÿ k"0 cot m α`kπ n for α ‰ kπ, n, m P N and n ě 2. Such sums arose in connection with certain central limit theorems arising in free probability see our paper [23], where the matrices considered below arise in a natural way. We apply the trace method to evaluate such expressions into closed form. For some specific values of α P R such formulas can be found in the literature, see, e.g., [20,15,47,8], in particular the question for which values of the parameters the sum (2.1) yields integer values is intriguing. The most general formula so far was given in [20], where the cotangent values cot α`kπ n were realized as roots of a polynomial and expressed via Cramer's rule applied to the Newton relations between elementary and power sum symmetric functions. In the present paper we go one step further and show that the polynomial found in [20] is in fact the characteristic polynomial of a simple matrix. We then exploit the trace method to draw certain conclusions about the sum (2.1) and compute closed formulas. It is a well known fact that the trace of a matrix equals the sum of its eigenvalues counting algebraic multiplicities. This relation is respected by functional calculus and the identity Tr f pAq " ÿ f pλ i q holds for arbitrary holomorphic (and other functions in the case of self-adjoint matrices), in particular, powers and polynomials. The trace method consists in the evaluation of this identity for particular matrices in order to obtain nontrivial combinatorial relations. It is perhaps interesting to note that the papers [12,13] evaluate certain trigonometric sums using matrices with trigonometric entries and integer eigenvalues, while in the present paper we exploit integer matrices with trigonometric eigenvalues.

Preliminaries on Linear Algebra and the Tangent function
The main role in this paper is played by a certain matrix and its intricate relations to the tangent and cotangent functions.
2.1. A matrix. For scalars a, b, c P C we denote by r c a b c s n P M n pCq the matrix whose diagonal elements are equal to c, whose upper-triangular entries are equal to a and whose lower-triangular elements are equal to b, respectively. For simplicity of notation, we use the same letter J n and B n for the following matrices (2.1) χ n pα; λq " pcot α`iqpλ´iq n´p cot α´iqpλ`iq n 2i " Impcot α`iqpλ´iq n (assuming λ real) and the eigenvalues are given by Proof. The spectrum of the matrix C n can be computed from its characteristic polynomial χ n pα; λq " detpλI´C n q using the following recurrence relation. Let w " a`i , then we have 0´s w´s w´s w . . . λ´a " p2λ´2a`w`s wqχ n´1 pα; λq´pλ´a`wqpλ´a`s wqχ n´2 pα; λq and the solution of this recurrence equation (with initial values χ 0 pα; λq " 1 and χ 1 pα; λq " λ) is wpλ´a`s wq n´s wpλ´a`wq n w´s w " pa`iqpλ´iq n´p a´iqpλ`iq n 2i " Impa`iqpλ´iq n .
Thus we have to solve the equation To compute the zeros, write a`i " r 0 e iα , i.e., a " cot α and assume λ´i " re´i θ . Then equation (2.2) becomes Im r 0 e iα r n e´i nθ " 0 and is equivalent to the equation sinpα´nθq " 0, that is, α´nθ "´kπ for some k P Z. Thus the solutions of (2.2) can be written as λ k " i`r k e´i θ k with θ k " α`kπ n . Now our matrix is selfadjoint, all roots of the characteristic polynomial (2.2) are real and hence´1 " Impλ k´i q "´r k sin θ k ; we conclude that r k " 1 sin θ k and λ k " Repλ k´i q " r k cos θ k " cot θ k . Consequently 2. An alternative formula for this polynomial can be found in [20,Formula (4)]. Indeed the coefficients of this polynomial are as follows (2.4) χ n pα; λq " Impa`iqpλ´iq n " Im cf. [20,Formula (4b)].
In fact the discussion of [20] starts by showing that the characteristic polynomial χ n pα; xq is related to the expression sin arccot x. Indeed evaluation of the polynomial (2.1) at λ " cot θ and few elementary manipulations yield the identity χ n pα; cot θq " sinpnθ´αq sin α sin n θ .

2.2.
Formulas for tanpnxq. A simple manipulation of the addition formulae for sine and cosine show that the tangent function obeys the addition rule This rule is not practical for iteration and the following equivalent elegant formula proposed by Szmulowicz [45] is a convenient alternative It follows immediately from the identity and in particular, tanpn arctan zq is a rational function. Indeed " i p1´i tan xq n´p 1`i tan xq n p1´i tan xq n`p 1`i tan xq n (2.10) and cotpnxq " i pcot x`iq n`p cot x´iq n pcot x`iq n´p cot x´iq n . (2.11) Thus we obtain the well known formula [5, item 16] (2.12) tanpn arctan zq " i p1´izq n´p 1`izq n p1´izq n`p 1`izq n ; comparing with the reciprocal polynomial of (2.1) at a " cot α " 0 which is p n pzq " z n χ n p0; 1{zq " p1´izq n`p 1`izq n 2 we see that (2.13) tanpn arctan zq "´1 n`1p The reciprocal provides the following crucial identity for cot (2.17) cotpnx´αq "´i pcot α`iqp1´izq n`p cot α´iqp1`izq n pcot α`iqp1´izq n´p cot α´iqp1`izq n which after comparison with the reciprocal polynomial χ n pα; zq " z n χ n pα; 1{zq " pcot α`iqp1´izq n´p cot α´iqp1`izq n 2i identifies to (2.18) cotpn arctan z´αq " 1 n`1χ 1 n`1 pα; zq χ n pα; zq 2.4. Derivatives of tan and cot. The higher derivatives of tan z and cot z are closely related, since cot z " tan`π 2´z˘. It is easy to see that there exist polynomials P n pzq such that d n dz n tan z " P n ptan zq; indeed these derivative polynomials satisfy the recursion P n`1 pxq " p1`x 2 qP 1 n pxq and can be used to efficiently compute tangent and Bernoulli numbers [33]. Explicitly, these polynomials can be expressed via the geometric polynomials [10, (2.1)] as follows, see [10, (3.10-11)]: (2.20) On the other hand (see [ The higher order tangent numbers [14] are defined as coefficients of the series Their bivariate generating function On the other hand, from the addition formula (2.5) we infer the exponential generating function of the derivative polynomials to be x`tan z 1´x tan z .
Comparing the two generating functions we find the relation On the other hand let us denote by A pkq n the arctangent numbers (see [16, p. 260] or [19]) defined by their exponential generating function up to sign these are the same as the coefficients of the hyperbolic arctangent function pkq n n! z n the latter are nonnegative and (2.28) A pkq n " p´iq k i nÃpkq n . 2.6. Derivatives of arctan. The derivatives of arctan z are rational functions and it is easy to verify by induction that they are given by the following formulas and thus In this section we briefly recall the combinatorics behind the composition of exponential generating functions. Let pa n q ně1 and pb n q ně1 be sequences and define a new sequence by their combinatorial convolution Then Fa di Bruno's formula [42, Theorem 5.1.4] asserts that their exponential generating functions F a pzq " Equivalently, given smooth functions f and g, the m-th derivative of the composed function is

Trace formula
Finite trigonometric sums of the sort (2.1) with offset α " 0, i.e., n´1 ÿ k"1 cot m kπ n are a recurrent theme in the mathematical literature. Sums of this type arise in number theory in connection with Dedekind sums and topology [49,29], and more recently were used to evaluate the Riemann zeta function, see [48,Problem 141ff] for the apparently first occurrence of this connection and later rediscoveries [32,46,37,6,3,21]; Berndt and Yeap [8] attribute the first occurence of our power sums to [44, p. 155]. In more recent times the sum (2.1) was evaluated for the offsets α " 0 [8,21,2,22,26,47,28] and α " π{4 ( [11], see the discussion in Remark 3.3 below). Alternatively, for quadratic forms finite Fourier analysis is applicable [4]. The first to investigate the sum (2.1) for general offsets α was [20] and this is also the subject of the present paper. To this end it turns out to be convenient to express (2.1) as a trace of a certain matrix, which results in a very compact formula and allows to draw a number of conclusions about (2.1). For example, if cot α is an integer, e.g., α " π 4 , it follows trivially that (2.1) evaluates to an integer, as was observed by different means in [11]. For an evaluation of cosecant sums via the trace method see [43].
Proof. It is clear that the trace (3.1) is a polynomial of degree at most n in cot α. Moreover since the entries of the matrices J n and B n are integers, the coefficients p m,m´2k pnq are integers as well. For positivity, we show that the mixed moments of J n and B n are positive. To see this, note that P n " 1 n J n is a self-adjoint projection of rank 1. It follows that for any matrix C n the compression P n C n P n lies in the 1-dimensional algebra generated by P n , more precisely, P n C n P n " ξ T C n ξP n where ξ " 1 ? n p1, 1, . . . , 1q T spans the image of P n . For our matrix B n clearly ξ T B n ξ " ř b ij " 0 and by antisymmetry, also for odd powers ξ T B k n ξ " p´1q k ξ T B k n ξ " 0. It follows that any mixed moment n¨¨¨J kr n B lr n q " n k 1`¨¨¨`kr TrpP k 1 n B l 1 n P k 2 n B l 2 n¨¨¨P kr n B lr n q " n k 1`¨¨¨`kr TrpP n B l 1 n P n B l 2 n P n¨¨¨Pn B lr n P n q " (1) In particular, Spm, n, αq evaluates to an integer (natural number) whenever cot α is an integer (natural number). It was shown in [11] to the surprise of the authors that ř n k"1 p´1q k´1 cot 2m´1 p2k´1qπ 4n and ř n k"1 cot 2m p2k´1qπ 4n can be represented as integer-valued polynomials in n of degrees 2m´1 and 2m respectively. We will show that these facts are elementary consequences of the trace formula. Applying Lemma 2.1 to the matrix " 1 1`i 1´i 1 ‰ n , we obtain its eigenvalues as λ k " cotˆπ 4n`k n π˙, for k P t1, . . . , nu, because α " arccotp1q " π 4 . The identities of Byrne and Smith [11] involve the sum ř n k"1 cot r`π 4n`k n π˘; indeed cot rˆπ 4n`k n πu sing cotp π 4n`k n πq "´cotp´π 4n`n´k n πq for the second sum, we get if r is even.
We will see later that even for noninteger values of cot α the sum may evaluate to an integer, e.g., for n " 2 and cot α " 1 2 , Lucas numbers appear, see (5.1) below. (2) The second part of Theorem 3.1 could be seen as a very special case the BMV conjecture [40]: if A and B are positive semi-definite matrices, then for all positive integers m, the polynomial in t, TrpA`tBq m , has only non-negative coefficients. The proof above shows that the assertion is also true whenever A is an orthogonal projection of rank one and B is a positive or antisymmetric self-adjoint matrix. (3) From the Newton identities between power sum and elementary symmetric polynomials we conclude in situation when |B| " 1, then we reduce to Theorem 6.1 with m " 1 or when For |B| " n this is confirmed by the well known identities n´1 ź k"0 sinˆk π n`z˙" 2 1´n sinpnzq and n´1 ź k"0 cosˆk π n`z˙" 2 1´n sin´nz`π 2 n¯.
For other literature about trigonometric multiple cotangent sum similar to those in (3.3), we refer the reader to [8, Section 6] and [49].

Generating functions
In the present section we compute the generating function of the cotangent sums (2.1), for fixed n, i.e., which is the moment generating function of the matrix cot αJ n`Bn with respect to the nonnormalized trace. Moreover we will compute the moment generating function of the matrix B n with respect to the nonnormalized trace and with respect to the state ω with density matrix P n " 1 n J n , that is, ωpCq " TrpP n Cq " 1 n where as above by ξ we denote the unit vector ξ " 1 ? n p1, 1, . . . , 1q T and C " rc i,j s n i,j"1 P M n pCq. The moment generating functions M xJn`Bn pzq " TrppI´zpxJ n`Bn qq´1q, M Bn pzq " TrppI´zB n q´1q, with respect to the trace are easy to compute directly through the characteristic polynomials. On the other hand, direct computation of M Bn pzq " ωppI´zB n q´1q " TrpP n pI´zB n q´1q requires information about the eigenvectors which we could not obtain. It will therefore be computed indirectly. The tangent function and its inverse will play a major role in these computations and we collect some facts about these functions first.
where θ k " α`kπ n . More generally, the moment generating function of the matrix pencil xJ n`Bn is x`tanpn arctan zq 1´x tanpn arctan zq˙.
Proof. Once we have realized cot θ k as roots of a polynomial, it is easy to write down the generating function of the sequence (2.1) as a logarithmic derivative. Indeed, let g n pzq " n´1 ÿ k"0 1 z´cot θ k " χ 1 n pα; zq χ n pα; zq " n pcot α`iqpz´iq n´1´p cot α´iqpz`iq n´1 pcot α`iqpz´iq n´p cot α´iqpz`iq n then the ordinary generating function is pcot α`iqp1´izq n´p cot α´iqp1`izq n " n 1`z 2ˆ1`i z pcot α`iqp1´izq n`p cot α´iqp1`izq n pcot α`iqp1´izq n´p cot α´iqp1`izq n" n 1`z 2 p1´z cotpn arctan z´αqq where in the last step we used identity (2.17). The general formula (4.4) follows by substituting α " arccot x and the addition formula for tangent (2.5).

4.2.
A functional relation. In this section we indicate an algorithm to calculate the coefficients p m,m´2k pnq, which is the main contribution of this paper. The following lemma is a special case of cyclic Boolean convolution [35]; we reproduce the calculation here for the reader's convenience.  TrppxJ n`Bn q m q " Tr´pxJ n q m`Bm ǹ ÿ kě1 p 0 ě0 p 1 ,p 2 ,...,p k ě1 q 1 ,q 2 ,...,q k ě1 p 0`q1`p1`¨¨¨`qk`pk "m B p 0 n pxJ n q q 1 B p 1 n pxJ n q q 2 B p 2 n¨¨¨p xJ n q q k B p k ǹ ÿ kě1 q 0 ě0 p 1 ,p 2 ,...,p k ě1 q 1 ,q 2 ,...,q k ě1 q 0`p1`q1`¨¨¨`pk`qk "m pxJ n q q 0 B p 1 n pxJ n q q 1 B p 2 n pxJq q 2¨¨¨B p k n pxJ n q q k" TrpB m n q`TrppxJ n q m q ÿ kě1 p 0 ě0 p 1 ,p 2 ,...,p k ě1 q 1 ,q 2 ,...,q k ě1 p 0`q1`p1`¨¨¨`qk`pk "m pxnq q 1`q2`¨¨¨`qk TrpP B p 1 n q TrpP B p 2 n q¨¨¨TrpP B p k´1 n q TrpP B p k`p0 n q ÿ kě1 q 0 ě0 p 1 ,p 2 ,...,p k ě1 q 1 ,q 2 ,...,q k ě1 q 0`p1`q1`¨¨¨`pk`qk "m pxnq q 0`q1`¨¨¨`qk TrpP B p 1 n q TrpP B p 2 n q¨¨¨TrpP B p k n q   with initial condition gp0q " 1 has unique solution gpzq " zM Bn pzq " 1 n tanpn arctan zq. Proof. Observe that the considered expression can be rewritten to the first order linear equations on the standard form g 1 pzq`qpzqgpzq " ppzq. The method involves construction of an explicit solution to show that the associated integral equation has unique solution. So if we substitute gpzq " tanpn arctan zq n we see that this function satisfy equation from Lemma 4.4.

Combinatorial interpretation
In this section we indicate explicit combinatorial interpretations of the coefficients of polynomials (3.2) which express the value of trace of matrices in terms of Dyck paths and rooted binary trees. We emphasize that these coefficients p m,k pnq are nonzero, whenever m and k have the same parity.

Dimension 2.
First let us record that for n " 2 at offset cot α 0 " 1 2 we recover the well known sequence Lucas numbers (A000032 in the encyclopedia of integer sequences [39]). Indeed, the characteristic polynomial (2.1) is and the roots are the golden ratios φ˘" 1˘?5 2 with moments (5.1) Spm, 2, α 0 q " L m " φ m`φḿ satisfying the recurrence relation

5.2.
Dyck paths -interpretation of p 2m`1,1 pnq. For the general case we establish some recurrence relations. An explicit formula will be established in Corollary 6.6 below.
Proposition 5.1. The moments which is reminiscent of the recurrence relations for the Motzkin numbers.
Proof. The function Q n pzq " tanpn arctan zq z is rational by . Indeed Q 1 pzq " 1 and for higher order the addition theorem for the tangent function (2.5) yields the recurrence Q n pzq " 1 z tanparctan z`pn´1q arctan zq " z`tanppn´1q arctan zq z´z 2 tanppn´1q arctan zq " 1`Q n´1 pzq 1´z 2 Q n´1 pzq or eqivalently From Lemma 4.4 we infer that ř 8 m"0 d n,m z 2m " Q n pzq and we can readily calculate the required recurrence for the moments d n,m .
The continued fraction of the rational function Q n pzq is finite and was computed in [36]: We can thus infer from Flajolet's theory of continued fractions [24] the following formula for the moments d n,m .  where the sum runs over Dyck paths of length at most 2m with weights a k´1 " n´k 2k´1 , b k " n`k 2k`1 , k " 1, 2, . . . , n.
Example 5.3. For n " 3 the generating function is which is reminiscent of the Catalan recurrence relations.
Definition 5.4. A rooted binary tree is a rooted tree in which each node has at most two children, one of which we distinguish as firstborn. We use the convention that the root is not a child and therefore does not count as a firstborn; our trees are unordered but we take the convention that firstborns are always drawn on the right. We denote by T n,m the set of rooted binary trees with m leaves, such that each leaf has a brother and every path emanating from the root contains at most n´1 firstborns. We note that the set T n,m is empty unless m ď n. For a rooted binary tree τ P T n,m we denote by Pathspτ q the set of maximal rooted paths. For such a path p P Pathspτ q we denote by rppq be the number of firstborn nodes occuring in p and its weight ωppq " n´rppq which is a number between 1 and n. Proof. Let us denote the right-hand side of (5.5) by c n,m . If m " 1 the root is the only node and does not count as a firstborn, therefore e n,1 " c n,1 " n. Moreover T 2,2 only contains one tree of weight e 2,2 " c 2,2 " 2. More generally T n,2 contains pn´1q trees and e n,2 " c n,2 " npn´1q`pn´1qpn´2q`¨¨¨`2ˆ1. So, it is sufficient to verify that e n,m " c n,m for n, m ě 3. Notice that any rooted binary tree can be viewed as one or two (non-empty because n, m ě 3) rooted binary trees grafted onto a common root; see T n´1,m T n,m´k T n´1,k Case 1. Assume that the root has only one child v 0 . Then every path from the root to a leaf with at most n´1 firstborns consists of the first step and a path from v 0 with at most n´2 firstborns. So the weight remains the same and the number of leaves remains m. Case 2a). Let τ be such a tree and p a path passing through the firstborn child v 0 . Then we can consider the latter as root vertex of new binary tree in T n´1,k with k leaves for k P t1, . . . , m´1u. Denote by p 1 the restriction of the path p to this subtree. Observe that p 1 contains at most n´2 firstborns because v 0 already counts as a firstborn and the weights of p and p 1 coincide. Indeed rppq " rpp 1 q`1 and so ωppq " n´rppq " n´1´rpp 1 q " ωpp 1 q. Case 2b). Let now p be a path passing through the other child, that is, p 1 is a path in a tree from T n,m´k and again the weight does not change.

3ˆ2 2ˆ1
Finally we have and now we can write Thus we see that c n,m " e n,m .

5.4.
Interpretation of p m,k pnq for k ě 2. The combinatorial objects that we consider now are called circular binary forests.
Definition 5.6. Assume that m, k P N have the same parity. For k P t2, . . . , mu a circular binary forest T E n,m,k of degree k is a set of k binary trees as above arranged on a circle with a total number of m leaves, see Figure 5.3 for an example. The weight of a forest F " pτ 1 , τ 2 , . . . , τ k q is the product ωpF q " ź τ PF ωpτ q.
Proposition 5.7. Assuming that m and k pk ‰ 0q have the same parity, then Proof. From the proof of Theorem (3.1), we see that p m,k pnq " ÿ l 0 ,l 1 ,...,l k ě0 l 1 ,...,l k´1 ,l 0`lk even ř l i "m´k TrpB l 0 n J n B l 1 n J n B l 2 n¨¨¨J n B l k n q ÿ l 0 ,l 1 ,...,l k ě0 l 1 ,...,l k´1 ,l 0`lk even ř l i "m´k This can be visualized in terms of forests, see Figure 5.3.

5.5.
Interpretation of TrpB 2m n q -constant term. The moment generating function of B n is M Bn pzq " np1`z tanpn arctan zqq 1`z 2 " n`nz 2 Q n pzq 1`z 2 .
If we expand the generating function in powers of z, then we obtain TrpB 2m n q " nd n,m´1´T rpB 2m´2 n q " nd n,m´1´n d n,m´2`n d n,m´3`¨¨¨`n d n,0 p´1q m´1`n p´1q m Proof. We start by expressing the generating function (4.3) in terms of the functions f pzq " lnp|sinpz´αq|q and gpzq " n arctan z. Indeed observe that n 1`z 2 p1´z cotpn arctan z´αqq " and moreover the Leibniz rule of order m implies moreover if π is odd then |π| " m mod 2 and we obtain finally the Taylor coefficients of n 1`z 2 contribute ni m for even m and the claim follows.
Corollary 6.2. The cotangent sums (2.1) evaluate to where A pkq m are the arctangent numbers (2.26); note that these are alternating (2.28). Proof. We extract the essential part of the formula (6.1) and arrive at the expression ÿ πPP odd pmq P |π|´1 pcot αqpniq |π| µp0 m , πq " hence by (2.31) F g pF f pzqq " e t atanh z´1 and the coefficient of t k yields the desired coefficient c m,k "Ã pkq m and from (2.28) we gather the correct sign. Remark 6.3. Comtet [16, p. 260] asserts that the arctangent numbers are inverse to the derivative polynomials. This means that the standard monomials can be expanded as a linear combination of tangent polynomials as follows [19, Formula (2.14)]: Let us explain now that the similarity of this formula with (6.3) is not a coincidence. Indeed using the property that the derivative polynomials linearize the cotangent power and the simple formula Sp1, n, αq " Tr C n " n cot α allow for the following alternative straightforward proof: Finally let us evaluate formula (6.3) at offsets α " π{4 and α " π{2.
Proof. The evaluation of the generating function (2.24) yields P p0, zq " tan z P p1, zq " 1`tan z 1´tan z " tanp2zq`secp2zq and we conclude that P n p0q " T n and P n p1q " 2 n E n " E B n which are also known as Euler numbers of type B [34]. Corollary 6.6. Extracting the linear coefficient of (6.3) we can obtain an explicit expression for the moments (5.2) (6.7) TrpJ n B 2m n q " 1 p2mq!˜A

7.2.
Asymptotic analysis and derivative. In order to investigate asymptotic properties formula from Theorem 6.1 it is sufficient to consider the singleton partition and we obtain lim nÑ8 1 n m n´1 ÿ k"0 cot m α`kπ n " " 1 pm´1q! P m´1 pcot αq if m ą 1 cot α if m " 1.
Approximation of the Riemann zeta function for even values by powers of cotangent is well studied, see [46,3,21].