# Convert a two-ququart (16 x 16) density matrix into normal form--so that the components of the Bloch vectors of the two reduced systems are all zero

The two-ququart ($$16 \times 16$$) "Hiesmayr-Loffler" density matrix https://iopscience.iop.org/article/10.1088/1367-2630/15/8/083036/meta, (https://arxiv.org/abs/2004.06745 eq. (13)), What are the ranges of the four $q$ parameters in the magic simplex of Bell states formula?) $$\begin{equation} \rho_{HL}^{2qq}= \end{equation}$$ $$\begin{equation} \left( \begin{array}{cccccccccccccccc} \kappa _1 & 0 & 0 & 0 & 0 & \kappa _2 & 0 & 0 & 0 & 0 & \kappa _2 & 0 & 0 & 0 & 0 & \kappa _2 \\ 0 & Q_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & Q_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \kappa _3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \kappa _3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \kappa _2 & 0 & 0 & 0 & 0 & \kappa _1 & 0 & 0 & 0 & 0 & \kappa _2 & 0 & 0 & 0 & 0 & \kappa _2 \\ 0 & 0 & 0 & 0 & 0 & 0 & Q_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & Q_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & Q_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \kappa _3 & 0 & 0 & 0 & 0 & 0 & 0 \\ \kappa _2 & 0 & 0 & 0 & 0 & \kappa _2 & 0 & 0 & 0 & 0 & \kappa _1 & 0 & 0 & 0 & 0 & \kappa _2 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & Q_2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & Q_2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & Q_3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \kappa _3 & 0 \\ \kappa _2 & 0 & 0 & 0 & 0 & \kappa _2 & 0 & 0 & 0 & 0 & \kappa _2 & 0 & 0 & 0 & 0 & \kappa _1 \\ \end{array} \right), \end{equation}$$ where, $$\kappa_1=\frac{1}{4} \left(Q_1+3 Q_4\right),\kappa_2=\frac{1}{4} \left(Q_1-Q_4\right)$$ and $$\kappa_3=\frac{1}{4} \left(-Q_1-4 Q_2-4 Q_3-3 Q_4+1\right)$$ is not in normal form, as the Bloch vectors of the two reduced $$4 \times 4$$ subsystems both have a component $$\frac{1}{16} \sqrt{\frac{3}{2}} (Q_1+3 Q_4)$$ associated with the fifteen generator of $$SU(4)$$, $$\begin{equation} \lambda_{15}=\left( \begin{array}{cccc} \frac{1}{\sqrt{6}} & 0 & 0 & 0 \\ 0 & \frac{1}{\sqrt{6}} & 0 & 0 \\ 0 & 0 & \frac{1}{\sqrt{6}} & 0 \\ 0 & 0 & 0 & -\sqrt{\frac{3}{2}} \\ \end{array} \right). \end{equation}$$

The components associated with the fourteen other generators are all zero in both cases.

In other words, $$a_{\mu}= \mbox{Tr}[\rho_{HL}^{2qq} (\lambda_{\mu} \otimes \mathbb{1})]=0, b_{\mu}= \mbox{Tr}[\rho_{HL}^{2qq} (\mathbb{1} \otimes \lambda_{\mu})]=0,\mu=1,\ldots,14$$ and $$a_{15}= \mbox{Tr}[\rho_{HL}^{2qq} (\lambda_{15} \otimes \mathbb{1})]=b_{15}= \mbox{Tr}[\rho_{HL}^{2qq} ( \mathbb{1} \otimes \lambda_{15})]\frac{1}{16} \sqrt{\frac{3}{2}} (Q_1+3 Q_4)$$.

I would like to convert $$\rho_{HL}^{2qq}$$ to normal form (https://journals.aps.org/pra/abstract/10.1103/PhysRevA.68.012103), that is, where all fifteen of the components of the Bloch vectors of the two reduced ququart ($$4 \times 4$$) density matrices become zero.

Given such a normal form, I can then pursue the application to it of the necessary and sufficient conditions for separability of Li and Qiao

So doing, would require the construction of the $$15 \times 15$$ correlation matrix $$T_{\mu \nu}= \mbox{Tr}[\rho_{HL}^{2qq} (\lambda_{\mu} \otimes \lambda_{\nu})]$$, and, then, the determination of its singular-value decomposition.

Here is some Matlab code--which I will try to implement--written by Nathaniel Johnston http://www.qetlab.com/FilterNormalForm for this explicit purpose (however, this may work only with numerical--not symbolic--input):

%%  FILTERNORMALFORM    Computes the filter normal form of an operator
%   This function has one required argument:
%     RHO: a density matrix
%
%   XI = FilterNormalForm(RHO) is a vector of the coefficients in RHO's
%   filter normal form (see Section IV.D of ), which are useful for
%   showing that RHO is entangled.
%
%   The filter normal form is not guaranteed to exist if RHO is not full
%   rank. If a filter normal form can not be found, an error is returned.
%
%   This function has two optional input arguments:
%     DIM (default has both subsystems of equal dimension)
%     TOL (default sqrt(eps))
%
%   This function has four optional output arguments:
%     GA,GB: cells of mutually orthonormal matrices
%     FA,FB: invertible matrices
%
%   [XI,GA,GB,FA,FB] = FilterNormalForm(RHO,DIM,TOL) returns XI, GA, GB,
%   FA, FB such that (eye(length(RHO)) + TensorSum(XI,GA,GB))/length(RHO)
%   equals kron(FA,FB)*RHO*kron(FA,FB)'. In other words, FA and FB are
%   matrices implementing the local filter, XI is a vector of coefficients
%   in the filter normal form, and GA and GB are cells of matrices in the
%   tensor-sum decomposition of the filter normal form.
%
%   DIM is a 1-by-2 vector containing the dimensions of the subsystems on
%   which RHO acts. TOL is the numerical tolerance used when constructing
%   the filter normal form.
%
%   URL: http://www.qetlab.com/FilterNormalForm
%
%   References:
%    O. Gittsovich, O. Guehne, P. Hyllus, and J. Eisert. Unifying several
%   separability conditions using the covariance matrix criterion. Phys.
%   Rev. A, 78:052319, 2008. E-print: arXiv:0803.0757 [quant-ph]

%   requires: OperatorSchmidtDecomposition.m, OperatorSinkhorn.m,
%             opt_args.m, PartialTrace.m, PermuteSystems.m,
%             SchmidtDecomposition.m, Swap.m
%
%   author: Nathaniel Johnston (nathaniel@njohnston.ca)
%   package: QETLAB
%   last updated: October 3, 2014

function [xi,GA,GB,FA,FB] = FilterNormalForm(rho,varargin)

lrho = length(rho);

% set optional argument defaults: dim=sqrt(length(rho)), tol=sqrt(eps)
[dim,tol] = opt_args({ round(sqrt(lrho)), sqrt(eps) },varargin{:});

% allow the user to enter a single number for dim
if(length(dim) == 1)
dim = [dim,lrho/dim];
if abs(dim(2) - round(dim(2))) >= 2*lrho*eps
error('FilterNormalForm:InvalidDim','If DIM is a scalar, it must evenly divide length(RHO); please provide the DIM array containing the dimensions of the subsystems.');
end
dim(2) = round(dim(2));
end

try
[sigma,F] = OperatorSinkhorn(rho,dim);
catch err
% Operator Sinkhorn didn't converge.
if(strcmpi(err.identifier,'OperatorSinkhorn:LowRank'))
error('FilterNormalForm:NoFNF','The state RHO can not be transformed into a filter normal form. This is often the case if RHO is not of full rank.');
else
rethrow(err);
end
end

% Do some post-processing to make the output more useful and consistent
% with the literature.
pd = prod(dim);
[xi,GA,GB] = OperatorSchmidtDecomposition(sigma - trace(sigma)*eye(pd)/pd,dim);
xi = pd*xi;
FA = F{1};
FB = F{2};