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_4right),kappa_2=frac{1}{4} left(Q_1-Q_4right)$ and $kappa_3=frac{1}{4} left(-Q_1-4 Q_2-4 Q_3-3 Q_4+1right)$
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

https://www.nature.com/articles/s41598-018-19709-z https://idp.springer.com/authorize/casa? https://idp.springer.com/authorize/casa?

https://idp.springer.com/authorize/casa?redirect_uri=https://link.springer.com/article/10.1007/s11128-018-1862-5&casa_token=9VD5vRx0Ns4AAAAA:m0_N47thT2-Kg0sFZaSDDJWwojwCy5Wx9e-UBFPV_6N5yBqAz9WB9KG0bFzKi4tkP7UFpxzs49xRg0De

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 [1]), 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:
%   [1] 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 ([email protected])
%   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};