Quantum Entanglement&Purity Testing: A Graph Zeta Function Perspective

We assign an arbitrary density matrix to a weighted graph and associate to it a graph zeta function that is both a generalization of the Ihara zeta function and a special case of the edge zeta function. We show that a recently developed bipartite pure state separability algorithm based on the symmetric group is equivalent to the condition that the coefficients in the exponential expansion of this zeta function are unity. Moreover, there is a one-to-one correspondence between the nonzero eigenvalues of a density matrix and the singularities of its zeta function. Several examples are given to illustrate these findings.


Introduction
One of the most interesting and often discussed properties arising in quantum information theory is that of quantum entanglement [12], wherein a bipartite quantum system is described by a joint state in such a way that the state of one subsystem cannot be described independently of the other, no matter the physical distance between them. In the restricted case of pure states that we consider here, a state |ψ⟩ is called entangled if it cannot be written as a product state |ϕ⟩ ⊗ |χ⟩. States which do not possess this property are called separable. In recent decades, a number of criteria for separability have been developed, including the positive partial transpose (PPT) criterion [18,24] and k-extendibility [7,26], and the problem of determining whether a state is separable or entangled has been shown to be NP-hard in many cases [13]. Additionally, quantum algorithms which test for separability, such as the SWAP test [1] and its generalizations [3], are under continuous development.
Meanwhile, another perspective has arisen in which quantum properties are framed in a graphtheorectic setting. This study began with the work of Braunstein, Ghosh, and Severini [5], who defined the density matrix of a graph as the normalized Laplacian associated to it. In their work, they give a graph-theorectic criterion for the entanglement of the associated density matrix; however, not all den-sity matrices can be encoded into a graph in this way, thus limiting the field of applicability of this criterion. The work of Hassan and Joag [16] addresses this problem by associating an arbitrary density matrix to a weighted graph. They then define a modified tensor product of graphs in such a way that an arbitrary quantum state is a product state if and only if it is the density matrix of a modified tensor product of weighted graphs.
The perspective of Braunstein et al. produces a graph-theorectic separability criterion known as the degree criterion [4], which was shown to be equivalent to the PPT criterion in [17]. In the limited case of density matrices which can be associated to a graph in this sense, the PPT criterion can be replaced by a simple graph-theorectical criterion. Here we will take another step in this direction by establishing a connection between a family of bipartite pure state separability tests, the cycle index polynomials of the symmetric group, and a new graph zeta function associated to a density matrix. The crux of the argument is that this zeta function is the generating function for the cycle index polynomial of the symmetric group evaluated at the moments of the density matrix, and this is exactly the acceptance probability of these tests. Moreover, we show that there is a one-to-one correspondence between the singularities of this zeta function and the nonzero eigenvalues of the associated density matrix, thereby establishing a variant of the Hilbert-Pólya conjecture associated to this zeta function.
In Section 2, we review the separability tests defined in [3,22] related to previous work in [8,9,26]. We then review previous constructions of the density matrix of a graph but ultimately take a different approach by naturally assigning an arbitrary density matrix to a weighted graph in Section 3. The relevant graph zeta function is then introduced in Section 4, and in Section 5, we prove that our separability tests are equivalent to the condition that the expansion of the corresponding zeta function has unit coefficients. This establishes a simple graphtheorectic test for the entanglement of pure bipartite states. We also derive a determinant representation for this function, and use it to prove the correspondence between its singularities and the nonzero eigenvalues of the associated density matrix in Section 6. Finally, we give concluding remarks in Section 7.

Review of Separability Tests
The separability tests outlined in [3,22] are examples of G-Bose symmetry tests, where G is some finite group. Let U : G → U(H) be a unitary representation of G on the Hilbert space H. A state ρ is called G-Bose symmetric if where Π G := 1 |G| g∈G U (g) is the projection onto the G-symmetric subspace, or the space of states |ψ⟩ such that U (g) |ψ⟩ = |ψ⟩ for all g ∈ G. To be clear, in general both mixed and pure states may demonstrate this property. In both [14,Chapter 8] and [22], it was shown how to test a state for G-Bose symmetry using a quantum computer, and we review a special case of this procedure here.  The G-Bose symmetry property is equivalent to the condition Π G S |ψ S ⟩ 2 = 1. To test for this, we generate a superposition control state wherein we take a superposition over some set of computational basis elements labelled with a corresponding group element G, and this can be done in a general sense with a quantum Fourier transform, as discussed in [22]. This control state is used to implement a corresponding unitary U (g) if the control qubit is in the state |g⟩. Afterwards, applying an inverse quantum Fourier transform and measuring all control qubits, accept if the outcome |0⟩⟨0| C occurs and reject otherwise. The acceptance probability, given a pure state, is then given by This result is easily generalized to mixed states by convexity, and we find that the acceptance probability is equal to Tr[Π G S ρ S ]. Consider the circuit in Fig 1. To test for the separability of a bipartite pure state ψ AB , we prepare k copies of the state ψ AB and label the k copies of the A system by A 1 · · · A k and the k copies of the B system by B 1 · · · B k . We then perform an S k -Bose symmetry test on the state ψ ⊗k AB , wherein we identify S with A 1 B 1 · · · A k B k and U S (π) with is the projection onto the symmetric subspace. The acceptance probability p (k) = 1 for all k if and only if ψ AB is separable. Thus, the S k -Bose symmetry test is indeed a separability test for pure bipartite states.
In [3], this test is generalized to any finite group G. Moreover, the acceptance probability of this generalized test is shown to be given by the cycle index polynomial of G, which is defined for permutation groups, where c j (g) is the number of j-cycles in the cycle decomposition of g. This definition is easily extended to any finite group by Cayley's theorem [10]. When we specialize to the case G = S k , this acceptance probability takes the form where the sum is taken over the partitions of k.

Density Matrices and Weighted Graphs
The study of the density matrix of a graph was initiated by Braunstein et al. in [5]. A graph X is a set V (X) = {v 1 , . . . , v n } of labeled vertices along with a set E(X) of pairs {v i , v j } of vertices called edges.
The (vertex) adjacency matrix A X of the graph is given by and it encodes the information about the edges of the graph. The degree d(v i ) of a vertex v i is the number of edges which include the vertex v i . We define the degree matrix ∆ X of the graph to be the diagonal matrix consisting of the degrees of each vertex; that is, The combinatorial Laplacian of the graph is defined to be the difference between the degree matrix and the adjacency matrix: which is both positive semi-definite and symmetric; however, it does not have unit trace. For this reason, we define the density matrix of a graph to be We may prescribe an arbitrary orientation to the edges so that we may label them by e 1 , . . . , e |E(X)| and their respective inverses by e |E(X)|+1 , . . . , e 2|E(X)| . A path in the graph X is a sequence of edges a 1 · · · a s such that the origin vertex of a j+1 is the terminal vertex of a j . The primes of a graph are defined to be the equivalence classes of closed, backtrackless, tailless, primitive paths. By this we mean equivalence classes of paths such that the origin vertex of a 1 is the terminal vertex of a s , a s ̸ = a −1 1 , a j+1 ̸ = a −1 j , and such that the path is not the power of another path. The equivalence classes are given by the cyclic shifts of a 1 · · · a s . That is, we define [a 1 · · · a s ] = [a s a 1 · · · a s−1 ] = · · · = [a 2 · · · a s a 1 ].
Clearly, there are density matrices which do not fit the graph-theoretic description outlined above. This point has been addressed by Hassan and Joag [16] by generalizing the association of a density matrix to a graph to include a subset of weighted graphs. A weighted graph is a graph equipped with a weight function ω : E(X) → C. In what follows, we use the notation ω ij := ω({v i , v j }) for the weight of the edge connecting v i to v j . Define the adjacency matrix of a weighted graph (X, ω) by The degree matrix is defined the same way as before.
An approach similar to that in [16] is to define the density matrix of a weighted graph by where L X,ω := ∆ X,ω − A X,ω is the weighted Laplacian of the graph. For this to make sense, we need the weighted Laplacian to be nonzero, and graphically this is equivalent to the condition that the graph does not consist only of loops. In order to insure that the density matrix is positive semi-definite, we could restrict our attention to weight functions with magnitudes unchanged by a swap of the indices. That is, weight functions satisfying |ω ij | = |ω ji |. Note that this condition is satisfied by both symmetric and conjugate symmetric weight functions.
Proof. Let u ∈ C n be an arbitrary vector upon which L X,ω acts and denote its components by u i . Then Now the first term in (2) can be written as Similarly, the second term in (2) is given by Thus, we have from which it follows that the Laplacian is positive semi-definite.
Since (1) obviously has unit trace, it does indeed define a density matrix. Notice, however, that the entries of this density matrix are positive real numbers. By following the proof of Proposition 1, the reader will find that if |ω ij | is replaced by ω ij in the weighted adjacency matrix, then no conclusion can be made about the positive semi-definiteness of ρ X,ω . It is therefore not obvious how to extend this construction to take into account density matrices with complex coefficients. For this reason, we will take a different approach by naturally assigning an arbitrary density matrix to a weighted graph. This shortcoming is noted by Hassan and Joag in their version of the construction; although, they show that the property of having a density matrix is invariant under isomorphism of graphs. Moreover, several conditions for which a weighted graph does or does not have a density matrix are given in their work [16].
For our purposes, a simpler approach can be taken which includes all density matrices at the expense of having to exclude further graphs. Indeed, it is not hard to construct a weighted graph (X, ω) from a density matrix ρ. We construct it as follows: 1. The number of vertices is the dimension n of the Hilbert space upon which ρ acts. Label them 1, . . . , n. 2. If ρ ij is nonzero, include an edge from vertex i to vertex j with weight ω ij := ρ ij .
Notice that in the second point, an edge is included for both ρ ij and ρ ji , so that primes involving an edge from i to j followed by an edge from j to i can be considered (see Fig 2). It is this graph which we will use to construct our zeta function in the next section. In fact, the above procedure works for any matrix, not just density matrices. However, we will restrict our attention to those matrices which are positive semidefinite with unit trace so that relationships with quantum information can be developed. From this perspective, the only weighted graphs that we associate to a density matrix are those for which the matrix with entries ω ij given by the weight function is a density matrix. While it is immediately apparent that the weights of the loops satisfy ω ii ≥ 0 with i ω ii = 1, and the remaining weights satisfy ω ij = ω * ji , a complete classification of this subset of weighted graphs remains an open question.

Graph Zeta Functions
The most famous zeta function is that of Bernhard Riemann [11,23], which is defined as the function for Re(s) > 1 and its analytic continuation elsewhere. Associated to this function is the Riemann hypothesis, which states that all nontrivial zeros of the zeta function are contained on the vertical line with Re(s) = 1 2 . Proving this hypothesis is one of the most important problems in mathematics as many results rely on the assumption of its truth. One potential approach to solving this problem is the Hilbert-Pólya conjecture, which states that the imaginary parts of the nontrivial zeros of ζ are the eigenvalues of a self-adjoint operator. This function was connected to the theory of prime numbers by Euler upon proving the identity Subsequently, several analogous functions have been defined, including the Ihara zeta function [6,21,25] associated to a graph X, which is defined by where the product is over the primes in X, and ν(p) denotes the length (number of edges) of P . There is a determinant formula for this function which can be traced back to Bass [2] and Hashimoto [20] given by where r = |E(X)| − |V (X)| + 1 is the rank of the fundamental group of the graph.
A generalization of the Ihara zeta function called the edge zeta function can be defined as follows: For a graph X, there is an associated edge matrix W with entries given by the variable ω ab if the terminal vertex of edge a is the origin vertex of edge b and b ̸ = a −1 , and zero otherwise. The edge zeta function is defined by whereÑ E (C) = ω a 1 a 2 ω a 2 a 3 · · · ω asa 1 is the edge norm and C = a 1 · · · a s is a closed path in X.
One can formulate analogs to the prime number theorem and the Riemann hypothesis for these functions, making them interesting in their own right. Here we will define a different but related zeta function ζ ρ (u) associated to a density matrix. Indeed, we define where the product is over equivalence classes of primes in the weighted graph associated to ρ, ν(P ) again denotes the length of the prime representative P , and N E (P ) is the product of the weights of the edges which make up P . If the weights of the graph were all unity, then the Ihara zeta function would be recovered. Of course, this cannot be true for a density matrix with rank greater than one since its trace is unity. However, this point is worth noting when (3) is extended to all matrices, in which case the condition can certainly hold. Note that while (3) is a generalization of the Ihara zeta function, it differs from the natural generalization seen in the book by Terras [25]. However, by setting ω ab = ρ ij u where i denotes the origin vertex of edge a and j denotes the origin vertex of edge b (the terminal vertex of a), the density matrix zeta function is recovered. Indeed, the edge norm becomes the product of the weights of the closed path in the weighted graph associated to the density matrix multiplied by an extra factor of u ν(P ) . Therefore, the density matrix zeta function defined by (3) is a special case of the edge zeta function.
In the next section, we will show that the separability tests in [3] and [22] corresponding to the symmetric group are equivalent to the conditions for the zeta function of the graph associated to the reduced density matrix of a pure state.

Equivalence of Test with Zeta Function Criterion
Our main result relies on the fact that the zeta function associated to a density matrix is the generating function for the cycle index polynomial of the symmetric group evaluated at x j = Tr[ρ j ] for j = 1, . . . , n. To see this, we will derive an exponential expression for this function from (3). We will also prove the following theorem, which gives a determinate formula for this zeta function. Throughout, we will assume that |u| is small enough to force the convergence of the series involved. The proof is similar to that in [19] for the edge zeta function. Theorem 1. Let ρ be a density matrix and let ζ ρ denote the associated zeta function. Then Proof. Starting from the definition (3), we have ζ ρ (u) = [P ] (1 − N E (P )u ν(P ) ) −1 .

Now taking the logarithm, this becomes
where in the second equality we have expanded the logarithm, and in the third, the inner sum is now over all prime paths of length m (the equivalence class of such a prime contains m elements, explaining the factor of m in the denominator). Let us define an operator by L := k,l ρ kl ∂ ∂ρ kl and observe that N E (P ) = ρ i 1 j 1 · · · ρ i ν(P ) j ν(P ) for i 1 , . . . , i ν(P ) , j 1 , . . . , j ν(P ) the vertices in P . Then L(N E (P )) j = jν(P )(N E (P )) j , so that we have where we have used the fact that (N E (P )) j = N E (P j ) since P j is just j copies of each edge in P . Now, we are doing nothing more than summing over all the primes of any given length and of any given power. Thus, (5) is equivalent to a sum over all closed, backtrackless, tailless paths. That is, we can drop the primitive assumption and write (5) as where the sum is over the paths mentioned above. The ij-th component of the m-th power of uρ is given by and now we recognize the summand as the value of N E (C) for some path C = e ii 2 e i 2 i 3 · · · e i m−1 im e imj .
Letting i = j, this becomes a closed, backtrackless, tailless path; therefore, we have shown where o(C) denotes the origin vertex of C. Now summing over i, we have so that summing over m produces On the other hand, To finish off the proof, we note that with ρ ij = 0 for all i, j, we have ζ ρ (u) = 1 = det(I − uρ) −1 , so that applying the method of characteristics yields log ζ ρ (u) det(I − uρ) −1 = 0, which implies ζ ρ (u) = det(I − uρ) −1 .
Theorem 1 tells us that the reciprocals of the nonzero eigenvalues of ρ coincide with the zeros of 1/ζ ρ and therefore the singularities of ζ ρ . It thereby establishes a variant of the Hilbert-Pólya conjecture wherein, instead of the imaginary parts of the nontrivial zeros of a zeta function, we consider the singularities of the zeta function. Then the corresponding statement is that the singularities of this zeta function are given by the reciprocal eigenvalues of a self-adjoint operator, namely the matrix ρ. In fact, it is clear from (4) that ζ ρ has no zeros; although, it does vanish asymptotically.
Note that the matrix representation of a quantum state is basis dependent, and under a change of basis, the weighted graph associated to this matrix will change too. However, the zeta function is unchanged, and this fact follows from Theorem 1. Indeed, if P is a change of basis matrix and ρ ′ = P ρP −1 , then ζ ρ ′ (u) = det(I − uρ ′ ) −1 = det(P P −1 − uP ρP −1 ) −1 = det(P ) −1 det(I − uρ) det(P ) = det(I − uρ) = ζ ρ (u). We can therefore rest assured that ζ ρ (u) is well-defined. The next theorem shows that ζ ρ (u) is the generating function for the cycle index polynomial Z(S n )(1, Tr[ρ 2 ], . . . , Tr[ρ n ]). This is the key to the relationship between this zeta function and the separability tests discussed in Section 2.
Proof. By Corollary 1, we have the identity Let us split up this exponential into a product and then expand. This gives us Then the n-th coefficient is given by the sum over all terms where the exponent ∞ k=1 kj k of u is equal to n, each of which is given by a partition of n. That is, we have that the n-th coefficient is du n ζ ρ (u) u=0 = 1 for all n ∈ Z ≥0 . Proof. In Section 2 (and in [3,22]) it was shown that a bipartite pure state ψ AB is separable if and only if the acceptance probability of the symmetric group separability algorithm reviewed there is 1 for all n, and that this is equivalent to the statement that the state is separable if and only if Z(S n )(1, Tr[ρ 2 ], . . . , Tr[ρ n ]) = 1 for all n. By Theorem 2, the corollary then follows.
By considering the density matrix zeta function in the original form (3), Corollary 2 exchanges the computation and evaluation of the cycle index polynomial of S n for a graph-theorectic calculation. It follows that the tests in [3,22] can be viewed from a graph-theorectic perspective just as the PPT criterion was shown to be equivalent to the graphtheoretic degree criterion [4]. To illustrate this point, let us compute a few of the coefficients. For n = 0, the condition for separability holds trivially since (1 − N E (P )) −1 , but there are no primes with zero edges, so that the product evaluates to unity. For n = 1, we use the fact that d du ζ ρ (u) = ζ ρ (u) d du log(ζ ρ (u)) and note that d du log(ζ ρ (u)) = − d du Now letting u = 0, the only terms that survive are those with ν(P ) = 1, and we have The only primes with ν(P ) = 1 are the loops, which correspond to the diagonal entries of the density matrix. Therefore, [ d du ζ ρ (u)] u=0 is the sum of the diagonal entries of ρ; that is, [ d du ζ ρ (u)] u=0 = Tr[ρ] = 1, which is consistent with the cycle index polynomial calculation.
For the n = 2 case, the coefficient is given by evaluating 1 2 d 2 du 2 ζ ρ (u) at u = 0. Observe that 1 2 d 2 du 2 ζ ρ (u) is given by so that evaluating at u = 0 yields where we have used the previous results for n = 0, 1 to simplify. We obtain the graph-theorectic version of the second term from (3) in a similar way to the n = 1 case. Indeed, using the standard quotient rule and evaluating at u = 0 yields (N E (P )) 2 , and the n = 2 coefficient is therefore given by Now, the corresponding cycle index polynomial calculation produces the value 1 2 (1+Tr[ρ 2 ]). Therefore, it must be the case that (N E (P )) 2 , and this is easy to see when ρ is written in the diagonal basis (which exists by the spectral theorem). Indeed, in this case the graph consists only of loops so that the ν(P ) = 2 term vanishes and the ν(P ) = 1 term is the sum of the squares of the eigenvalues of ρ.
Consider the density matrix which is the pure state |+⟩⟨+|, so that all coefficients should be equal to unity. The weighted graph associated to ρ is given in Fig 2. As expected, the n = 1 coefficient is and the n = 2 coefficient is Next consider the maximally mixed state ρ = 1 2 I on one qubit. Since this state is mixed, every purification of it is entangled. The Bell state 1 √ 2 (|00⟩+|11⟩) is such a purification since tracing out either subsystem gives the state 1 2 (|0⟩⟨0| + |1⟩⟨1|) = 1 2 I. Therefore, by computing the n = 2 coefficient of ζ ρ (u), we get a graph-theorectic proof that the Bell state is entangled. Indeed, the weighted graph associated to ρ is given in Fig 3 and consists only of loops since the density matrix is diagonal. There are therefore no primes with ν(P ) = 2 and the only two primes with ν(P ) = 1 are the loops themselves. Thus, the n = 2 coefficient is is entangled. Note that the n = 2 case is the graph-theorectic equivalent to the well-known SWAP test [1,15]. It is perhaps worth noting that these tests for each n double as tests for the purity of a density matrix. The n = 2 case measures the purity exactly. In fact, this is really the mechanism behind the separability test since a pure bipartite state is separable if and only if its reduced state is pure.

Singularity-Eigenvalue Correspondence
Observe that det(I − uρ) is a polynomial in u of degree at most equal to the dimension of the Hilbert space upon which ρ acts. Then the number of zeros of the reciprocal zeta function 1/ζ ρ is given by the fundamental theorem of algebra, and this coincides with the number of singularities of ζ ρ . Suppose that ψ AB is a pure bipartite state and let ρ B denote its reduced density matrix given by tracing out the A-subsystem. If ρ B is pure, it is a rank one operator such that the only nonzero eigenvalue is λ = 1. Since the reciprocals of the nonzero eigenvalues of ρ B coincide with the zeros of 1/ζ ρ by Theorem 1, we have another criterion for the separability of pure bipartite states.
Corollary 3. Let ψ AB be a pure bipartite state and let ρ B denote its reduced density matrix by tracing out one of the systems. Then ψ AB is separable if and only if the only singularity of ζ ρ (u) is at u = 1.
Proof. By Theorem 1, we have Notice that if u = 0, then ζ ρ ̸ = 0. Then the condition that 1/ζ ρ (u) = 0 is equivalent to the eigenvalue problem Now, ψ AB is separable if and only if ρ B is a pure state, which is true if and only if the only nonzero eigenvalue is 1 u = 1, so that we have u = 1. This is equivalent to the only zero of the reciprocal zeta function being at u = 1, which is equivalent to the only singularity of ζ ρ B being at u = 1. To illustrate this point, we now examine the plots of the zeta function associated to different choices of density matrices. In Fig 4, we have an archetypal example of a pure state, |+⟩ = 1 √ 2 (|0⟩ + |1⟩), and the singularity at u = 1 is marked by a vertical line. Let us consider something a little more interesting. Since the W -state given by |W ⟩ = 1 √ 3 (|001⟩+|010⟩+|100⟩) is a pure state, its zeta function should look the same as that in Fig 4. However, this state is entangled, so that we expect to pick up at least one singularity which is not at u = 1 in the zeta function after the first system is traced out. This is exactly what we see in Fig 5. If we instead take ψ AB to be the GHZ-state given by |GHZ⟩ = 1 √ 2 (|000⟩ + |111⟩), we find a similar result to |W ⟩, except that the zeta function looks somewhat different. Here it turns out that both zeros of the reciprocal zeta function of the reduced state ρ B are located at u = 2, so that this is the only singularity that appears in  Another interesting choice for ρ B is given by the isotropic state that continuously transforms between the maximally entangled and maximally mixed states. Indeed, we define from which we recover the maximally entangled state on two qubits at p = 0 and the maximally mixed state at p = 1. A three dimensional plot of ζ ρ (u, p) is given in Fig 7. The maximally entangled state is a pure state, so that the cross section at p = 0 looks   We see from Fig 7 that as p varies from zero to one, a second pole comes from infinity and recombines with the original u = 1 pole when p = 1, where the state is maximally mixed; however, the merging of the poles does not happen at u = 1, as the first pole shifts to the right as p increases. Instead, the merging happens at the dimension of the Hilbert space (in this case 4), and this can be seen from (6) with p = 1, where the only non-zero eigenvalue is 1 4 . It should be noted that neither Corollary 2 nor Corollary 3 apply to mixed states in general. To see   By tracing out the two extra copies of σ, we recover the original state σ, which is mixed. The corresponding plots for ζ σ and its coefficients are shown in Fig 11 and Fig 12, respectively. This shows that the tests we have developed for pure states do not hold for mixed states. Therefore, another method will have to be used to check for entanglement in this case.

Conclusion
A connection between the zeta function ζ ρ defined in (3) and the separability of bipartite pure states was derived. However, this result does not carry over to the more general case where the joint state is mixed. Establishing a similar result in the mixed state setting is worth an investigation of its own. On paper, to test whether a pure bipartite state is entangled, it suffices to compute the n = 2 case since the reduced state has unit purity if and only if the joint state is separable. From the graphtheoretic perspective, this means that the only relevant primes when testing for entanglement are those with length ν(P ) = 1 or 2. The higher order tests in [3,22] are introduced because the acceptance probability decays as n → ∞, so when noise is introduced, it may be beneficial to perform a higher order test.
Through the determinant formula (4), it was shown that the nonzero eigenvalues of a density matrix ρ are in correspondence with the singularities of ζ ρ . Such a connection has the potential to transform the study of the spectra of random density matrices to a graph-theorectic setting, and we leave this to future work. The reader is directed to [25] for more about the connections between random matrix theory and the various graph zeta functions.

Data Availability Statement
The Mathematica code used in this work is available in the following GitHub repository: https://github.com/mlabo15/ZetaFunctions.