Computing coherence vectors and correlation matrices, with application to quantum discord quantification

Coherence vectors and correlation matrices are important functions frequently used in physics. The numerical calculation of these functions directly from their definitions, which involves Kronecker products and matrix multiplications, may seem to be a reasonable option. Notwithstanding, as we demonstrate in this article, some algebraic manipulations before programming can reduce considerably their computational complexity. Besides, we provide Fortran code to generate generalized Gell Mann matrices and to compute the optimized and unoptimized versions of the associated Bloch's vectors and correlation matrix, in the case of bipartite quantum systems. As a code test and application example, we consider the calculation of Hilbert-Schmidt quantum discords.


I. INTRODUCTION
Correlation functions are fundamental objects for statistical analysis, and are thus ubiquitous in most kinds of scientific inquires and their applications [1,2]. In physics, correlation functions have an important role for research in areas such as quantum optics and open systems [3,4], phase transitions and condensed matter physics [5,6], and quantum field theory and nuclear and particle physics [7]. Another area in which correlation functions are omnipresent is quantum information science (QIS), an interdisciplinary field that extends the applicabilities of the classical theories of information, computation, and computational complexity [8][9][10].
In order to define coherence vectors and correlation matrices, let us consider a composite bipartite system with Hilbert space H ab = H a ⊗ H b . Hereafter the corresponding dimensions are denoted by d s = dim H s for s = ab, a, b. In addition, let Γ s j , with Tr(Γ s j ) = 0 and Tr(Γ s j Γ s k ) = 2δ jk , (1) * Electronic address: jonas.maziero@ufsm.br be a basis for the special unitary group SU (d s ). Any density operator describing the state of the system H ab can be written in the local basis Γ a j ⊗ Γ b k as follows: where j = 1, · · · , d 2 a − 1 and k = 1, · · · , d 2 b − 1, and I s is the identity operator in H s . One can readily verify that the components of the coherence (or Bloch's) vectors a = (a 1 , · · · , a d 2 a −1 ) and b = (b 1 , · · · , b d 2 b −1 ) and of the correlation matrix C = (c j,k ) are given by: It is worthwhile mentioning that the mean value of any observable in H s , for s = a, b, ab, can be obtained using these quantities. In https://github.com/jonasmaziero/LibForQ.git, we provide Fortran code to compute the coherence vectors, correlation matrices, and quantum discord quantifiers we deal with here. Besides these functions, there are other tools therein that may be of interest to the reader. The instructions on how to use the software are provided in the readme file. Related to the content of this section, the subroutine bloch_vector_gellmann_unopt(d s , ρ s , s) returns the coherence vectors a or b and the subroutine corrmat_gellmann_unopt(d a , d b , ρ, C) computes the correlation matrix C. Now, let us notice that if calculated directly from the equations above, for d a , d b ≫ 1, the computational complexity (CC) to obtain the coherence vectors a and b or the correlation matrix C is: The remainder of this article is structured as follows. In Sec. II, we obtain formulas for a, b, and C that are amenable for more efficient numerical computations. In Sec. III we test these formulas by applying them in the calculation of Hilbert-Schmidt quantum discords. In Sec. IV we make some final remarks about the usefulness and possible applications of the results reported here.

II. COMPUTING COHERENCE VECTORS AND CORRELATION MATRICES
The partial trace function [42] can be used in order to obtain the reduced states ρ a = Tr b (ρ) and ρ b = Tr a (ρ) and to write the components of the Bloch vectors in the form: Thus, when computing the coherence vectors of the parties a and b, we shall have to solve a similar problem; so let's consider it separately. That is to say, we shall regard a generic density operator written as where s j = 2 −1 d s Tr(ρ s Γ s j ). Now, and for the remainder of this article, we assume that the matrix elements of regarded density operator ρ in the standard computational basis are given. We want to compute the Bloch's vector [43]: s = (s 1 , · · · , s d 2 s −1 ). For the sake of doing that, a particular basis Γ s j must be chosen. Here we pick the generalized Gell Mann matrices, which are presented below in three groups [44]: (13) which are named the diagonal, symmetric, and antisymmetric group, respectively. The last two groups possess d s (d s − 1)/2 generators each. Any one of these matrices can be obtained by calling the subroutine gellmann(d s , g, k, l, Γ s(g) (k,l) ). For the first group, g = 1, we make j = k and, in this case, one can set l to any integer.
It is straightforward seeing that, for the generators above, the corresponding components of the Bloch's vector can expressed directly in terms of the matrix elements of the density operator ρ s as follows: These expressions were implemented in the Fortran subroutine bloch_vector_gellmann(d s , ρ s , s). With this subroutine, and the partial trace function [42], we can compute the coherence vectors a and b.
We observe that after these simple algebraic manipulations the computational complexity of the Bloch's vector turns out to be basically the CC for the partial trace function. Hence, from Ref. [42] we have that for d a , d b ≫ 1, One detail we should keep in mind, when making use of the codes linked to this article, is the convention we apply for the indexes of the components of s. For the first group of generators, Γ s(1) j , naturally, j = 1, · · · , d s − 1.
We continue with the second group of generators, Γ The same convention is used for the third group of generators, Γ Next we address the computation of the correlation matrix C = (c j,k ), which is a (d 2 a − 1)x(d 2 b − 1) matrix that we write in the form: with the sub-matrices given as shown below. For convenience, we define the auxiliary variables: ι := 2 j(j + 1) , κ := 2 k(k + 1) , and ς := The matrix elements of C (1,1) , whose dimension is (d a − 1)x(d b − 1), correspond to the diagonal generators for a and diagonal generators for b: The matrix elements of C (1,2) , whose dimension is , correspond to the diagonal generators for a and symmetric generators for b: The matrix elements of C (1,3) , whose dimension is , correspond to the diagonal generators for a and antisymmetric generators for b: The matrix elements of C (2,1) , whose dimension is − 1), correspond to the symmetric gen-erators for a and diagonal generators for b: The matrix elements of C (2,2) , whose dimension is , correspond to the symmetric generators for a and symmetric generators for b: The matrix elements of C (2,3) , whose dimension is , correspond to the symmetric generators for a and antisymmetric generators for b: The matrix elements of C (3,1) , whose dimension is 2 −1 d a (d a − 1)x(d b − 1), correspond to the antisymmetric generators for a and diagonal generators for b: The matrix elements of C (3,2) , whose dimension is − 1), correspond to the antisymmetric generators for a and symmetric generators for b: The matrix elements of C (3,3) , whose dimension is −1), correspond to antisymmetric generators for a and antisymmetric generators for b: We remark that when implementing these expressions numerically, for the sake of mapping the local to the global computational basis, we utilize, e.g., The subroutine corrmat_gellmann(d a , d b , ρ, C) returns the correlation matrix C = (c j,k ), as written in Eq. (18), associated with the bipartite density operator ρ and computed using the Gell Mann basis, as described in this section. The convention for the indexes of the matrix elements c j,k is defined in the same way as for the coherence vectors. The computational complexity for C, computed via the optimized expressions obtained in this section, is, for d a , d b ≫ 1, By generating some random density matrices [45], we checked that the expressions and the corresponding code for the unoptimized and optimized versions of a, b, and C agree. Additional tests shall be presented in the next section, where we calculate some quantum discord quantifiers.

III. COMPUTING HILBERT-SCHMIDT QUANTUM DISCORDS
The calculation of quantum discord functions (QD) usually involves hard optimization problems [46,47]. In the last few years, a large amount of effort have been dedicated towards computing QD analytically, with success being obtained mostly for low dimensional quantum systems [48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65]. Although not meeting all the required properties for a bona fide QD quantifier [66], the Hilbert-Schmidt discord (HSD) [67], is drawing much attention due to its amenability for analytical computations, when compared with most other QD measures. In the last equation, the minimization is performed over the classical-quantum states with p j being a probability distribution, |a j an orthonormal basis for H a , ρ b j generic density operators defined in In this article, as a basic test for the Fortran code provided to obtain coherence vectors and correlation matrices, we shall compute the following lower bound for the HSD [68]: where λ a j are the eigenvalues, sorted in non-increasing order, of the (d 2 a − 1)x(d 2 a − 1) matrix: In the equation above t stands for the transpose of a vector or matrix. We observe that the other version of the HSD, D b hs , can be obtained from the equations above simply by exchanging a and b and using C t C instead of CC t .
It is interesting regarding that, as was proved in Ref. [69], a bipartite state ρ, with polarization vectors a and b and correlation matrix C, is classical-quantum if and only if there exists a (d a − 1)-dimensional projector Π a in the space R d 2 a −1 such that: Π a a = a and Π a C = C, Based on this fact, an ameliorated version for the Hilbert-Schmidt quantum discord (AHSD) was proposed [69]: with the matrix Υ a defined as where f and g are arbitrary functions of b ≡ ||b|| 2 . Then, by setting f (b) = g(b) = P (ρ b ) and using the purity, to address the problem of non-contractivity of the Hilbert-Schmidt distance, the following analytical formula was presented [69]: Thus both discord quantifiers D a hs and D a hsa are, in the end of the day, obtained from the eigenvalues λ a j . And the computation of these eigenvalues requires the knowledge of the coherence vector a (or b) and of the correlation matrix C. These QD measures were implemented in the Fortran functions discord_hs(ssys, d a , d b , ρ) and discord_hsa (ssys, d a , d b , ρ), where ssys = 's', with s = a, b, specifies which version of the quantum discord is to be computed. As an example, let us use the formulas provided in this article and the associated code to compute the HSD and AHSD of Werner states in H a ⊗ H b (with d a = d b ): where w ∈ [−1, 1] and F = da j,k=1 |jk kj|. The reduced states of ρ w are I s /d s , whose purity is P (ρ s ) = 1/d s . The results for the HSD and AHSD of ρ w are presented in Fig. 1. D a hsa (ρ w ) = daD a hs (ρ w ) = (daw − 1) 2 /((da − 1)(da + 1) 2 ). Due to the symmetry of ρ w , here D a hsa (ρ w ) = D b hsa (ρ w ). In the inset is shown the difference between the times taken by the two methods to compute the AHSD for a fixed value of da. We see clearly here that our optimized algorithm gives a exponential speedup against the brute force calculation of the Bloch's vectors and correlation matrix.

IV. CONCLUDING REMARKS
In this article, we addressed the problem of computing coherence vectors and correlations matrices. We obtained formulas for these functions that make possible a considerably more efficient numerical implementation when compared with the direct use of their definitions. We provided Fortran code to calculate all the quantities regarded in this paper. As a test for our formulas and code, we computed Hilbert-Schmidt quantum discords of Werner states. It is important observing that, although our focus here was in quantum information science, the tools provided can find application in other areas, such as e.g. in the calculation of order parameters and correlations functions for the study of phase transitions in discrete quantum or classical systems.

Conflict of Interests
The author declares that the funding mentioned in the Acknowledgments section do not lead to any conflict of interest. Additionally, the author declares that there is no conflict of interest regarding the publication of this manuscript.