# Prediction of the liquid-crystal phase behavior of hard right triangles from
fourth-virial DFT theories

Enrique Velasco<EMAIL_ADDRESS>Departamento de Física Teórica de la
Materia Condensada, Instituto de Física de la Materia Condensada (IFIMAC) and
Instituto de Ciencia de Materiales Nicolás Cabrera, Universidad Autónoma de
Madrid, E-28049, Madrid, Spain  Yuri Martínez-Ratón<EMAIL_ADDRESS>Grupo
Interdisciplinar de Sistemas Complejos (GISC), Departamento de Matemáticas,
Escuela Politécnica Superior, Universidad Carlos III de Madrid, Avenida de la
Universidad 30, E-28911, Leganés, Madrid, Spain

###### Abstract

We have used an extended Scaled-Particle Theory that incorporates four-body
correlations through the fourth-order virial coefficient to analyse the
orientational properties of a fluid of hard right-angle triangles. This fluid
has been analysed by computer simulation studies, with clear indications of
strong octatic correlations present in the liquid-crystal phase, although the
more symmetric order tetratic phase would seem to be the most plausible
candidate. Standard theories based on the second virial coefficient are unable
to reproduce this behaviour. Our extended theory predicts that octatic
correlations, associated to a symmetry under global rotations of the oriented
fluid by $45^{\circ}$, are highly enhanced, but not enough to give rise to a
thermodynamically stable phase with strict octatic symmetry. We discuss
different scenarios to improve the theoretical understanding of the elusive
octatic phase in this intriguing fluid.

Liquid crystals, hard right triangles, virial coefficients, density functional
theory.

##  I Introduction

Fluids of elongated particles in two dimensions (2D) deGennes  continue to
unveil surprising behaviours Musevic . Since the discovery of crystallisation
in the hard-disc systems, fluids made of hard elongated particles, such as
ellipsoids, have been seen to stabilise nematic phases with uniaxial symmetry
Frenkel . As common in 2D, these phases possess quasi-long-range orientational
order Frenkel1 ; Cuesta ; Odriozola , but mean-field theories can still
describe large regions where the tensor order parameter exhibits persistent
values Kleman ; Wensink . Since particles interact through purely overlap
interactions, order in 2D phases made of hard elongated particles is solely
governed by entropy, which shows here its most subtle nature.

More complicated 2D particle shapes have been explored more recently from a
theoretical perspective Schlacken ; MAR0 ; MAR ; Schmidt ; Triplett ; Wittmann
; Escobedo0 ; Sabi ; Quintana ; Lowen ; cinacchi , motivated by the
possibility to fabricate colloidal particles of virtually any shapeexp1 ; exp2
; exp3 . Particles with regular polygonal shapes have been demostrated to
exhibit mesophases Schlacken ; Anderson , with $4$-atic (or tetratic, with two
directors and a global symmetry under rotation by $90^{\circ}$) phases for
squares, $6$-atic (with three directors and a symmetry under $60^{\circ}$) for
equilateral triangles and hexagons. No further mesophases appear to get
stabilised for polygonal shapes with more edges: they crystallise directly
from the isotropic fluid through a KTHN-type transition, as in the case of
discs Anderson ; Mak . Simulation of mixtures of particles have also found
interesting behaviours Escobedo .

Nonregular polygons open up new possibilities MAR ; dani . Here entropy plays
an even more subtle role: particles tend to form local clusters of oriented
particles that can be viewed as ‘superparticles’, with symmetries sometimes
very different from that of the ‘monomers’ and therefore from the symmetry of
the bulk liquid-crystal phase that would trivially follow from the monomers.
Such is the case in fluids made of low-aspect-ratio rectangles, which tend to
form highly stable square clusters that stabilise a global $4$-atic phase
Granada . The same behaviour is found in vibrated monolayers of granular
grains Muller ; granos_nosotros  and in experiments on colloidal particles
exp3 . Hard-kite shaped particles have been studied by simulation sim  and
theory kites , and $4$-atic phases were also found. The basic understanding of
this phase lies in the excluded area between particles (second-order virial
coefficient), an essential ingredient of the Scaled-Particle Theory (SPT)
extension of classical Onsager theory. Three-body correlations can be
incorporated into the theory through the third-order virial coefficients
Mederos , and the ensuing corrections are important: the stability region of
the $4$-atic phase is extended to larger aspect ratios and lower densities.

Recently we have studied a fluid made of hard right-angle triangles (HRT) MAR1
. Motivated by Monte Carlo (MC) simulations by Gantapara et al. Gan , we
analysed the fluid using the standard SPT theory, based on the second virial
coefficient (which is analytic), and an extension that includes the third-
virial coefficient, calculated using MC integration. It turns out that none of
these theories can reproduce the behaviour predicted from the simulations: as
the isotropic fluid is compressed, clustering of particles in clusters of
various shapes give rise to strong $8$-atic correlations, and an orientational
distribution function with $4$-atic symmetry but high secondary peaks at
$45^{\circ}$ with respect to the main peaks results. The equilibrium
orientational function from the theories, by contrast, shows no hint of the
high-order $8$-atic symmetry.

The HRT fluid seems to be the first case where an Onsager-type theory fails to
give even a qualitative picture of a mesophase. In previous work we have
discussed this problem. On the one hand, the theory is strictly valid for
infinitely long rods, while here we are considering particle of low aspect
ratio. On the other, the sequence of virial coefficients in 2D fluids is known
to have peculiar behaviour, and the condition that scaled higher-order
coefficients are small is not fulfilled. Finally, a crucial property of the
HRT fluid and other 2D fluids made of hard polygonal particles is the strong
clustering tendency of the fluid MAR ; MAR1 , with at least five types of
clusters that form in the fluid at intermediate phases (before
crystallisation). Four of these clusters involve two particles (dimers), of
square, triangular and rhomboidal shapes (with two enantiomers in the latter
case), and one involve four particles arranged in a square. The fluid can
somehow be viewed as a multicomponent mixture of dynamic ‘superparticles’ and
might be more quantitatively described by association theories than by
particle (monomer)-based theories. Before undertaking such a programme, we
speculated MAR ; MAR1  that a theory based on four-body correlations (i.e. on
the fourth virial coefficient) might give some indication as to whether high-
order particle correlations, involved in clustering tendencies of the
particles, might be important to understand the equilibrium structure of the
fluid.

In the present work we show the predictions of such a theory. A resummed SPT
theory is developed using the standard technique Padilla ; Mederos ; MAR1 ,
which allows to systematically incorporate an arbitrary number of virial
coefficients. These objects are generalised virial coefficients in the sense
that they are functionals of the orientational distribution function. The
third and fourth virial coefficients are computed numerically, and the
instability of the isotropic (I) phase against orientational orders of
different symmetries is investigated. This process allows to analyse the
effect of increasing low-order, from two- to four-, particle correlations on
the onset of bulk orientational order. Focusing on the $8$-atic (or octatic)
orientational symmetry, we explore the tendency of the fluid to stabilise
orientational order through a bifurcation analysis. Our conclusion is that
four-particle correlations do enhance octatic symmetry. More definite
conclusions as to the orientational distribution of particles would require a
full minimization of the theories. However, a more quantitative theory should
incorporate particle clustering MAR1 , which is not possible with the present
SPT scheme, and further studies will have to await until a proper treatment of
clusters can be formulated.

The article is organised as follows. In Section II we present the theory and a
method to extend SPT to include the fourth virial coefficient. Also, we
provide some details on the numerical calculation of this coefficient. In
Section III we present the results, along with the results obtained from
extrapolations of virial coefficients and resummations of the virial theory in
the isotropic phase, which may help to explain the role of many-particle
correlations. Some conclusions are drawn in Section IV. The Appendix collects
some further numerical details and additional results.

##  II Theory

Our theory uses the same strategy followed in our previous works on the third
virial coefficient Mederos ; MAR1 , but now the theory is extended to an
arbitrary number of virial coefficients. We start from the following analytic
expression for the equation of state, valid for a fluid of oriented particles:

$\displaystyle\beta
pa=\frac{\eta}{1-\eta}+\frac{\displaystyle\sum_{k=2}^{n}c_{k}[h]\eta^{k}}{\displaystyle(1-\eta)^{q}}.$
(1)

In this expression $\eta=\rho a$ is the packing fraction, with $\rho$ and $a$
the number density and particle area, respectively. $h(\phi)$ is the
orientational distribution function, which quantifies the average orientation
of the particles, where $\phi$ is the angle between the averaged orientation
of the main symmetry axis of a typical particle and the main director. This
function is normalised,

$\displaystyle\int_{0}^{2\pi}d\phi h(\phi)=1.$ (2)

The coefficients $c_{n}[h]$ are functionals of the orientational distribution
function and can be related to the standard virial coefficients by equating
the low-density expansion of (1) with the exact virial expansion up to $n$th
order (note that the expansion in Eqn. (1) is truncated at this order):

$\displaystyle
c_{n}[h]=b_{n}[h]+\sum_{k=1}^{n-2}(-1)^{k}\binom{q}{k}b_{n-k}[h],$ (3)

where we have defined

$\displaystyle b_{n}[h]\equiv\frac{B_{n}[h]}{a^{n-1}}-1,$ (4)

with $B_{n}[h]$ the standard $n$th virial coefficient. Note that these
coefficients are also functionals of $h(\phi)$. The parameter $q$ in Eqn. (1)
controls the divergence of the second term at close-packing, $\eta=\eta_{\rm
cp}\equiv 1$, and will be important to compare the resulting equation of state
with MC simulations. Fixing $q=0$ the exact virial expansion up to the $n$th
order is recovered. Also, selecting $q=2$ and $n=2$ we obtain the SPT
approximation, while for the same $q$ and $n=3$ the equation of state turns
out to be the same as the one used in MAR1  to study fluid orientational
ordering close to the I-$p$–atic bifurcation point, henceforth called
$B_{3}$-SPT. This shows that Eqn. (1) is a versatile starting point to explore
nonstandard approximations to the equation of state of oriented fluids in two
dimensions.

In this article we start by exploring the case $q=2$ and $n=4$ (henceforth
called $B_{4}$-SPT), and then use extrapolated values for the first three
virial coefficients, $B_{2}$, $B_{3}$ and $B_{4}$ of the I phase to predict
the value of $B_{5}$, giving rise to the $B_{5}^{*}$-SPT theory (the asterisk
denotes an extrapolated value for $B_{5}$). Not only the values of the virial
coefficients for the completely disordered fluid are calculated, but also
relevant Fourier components with respect to weak orientational order of
particular $p$-atic symmetries. This allows for a bifurcation analysis of the
I phase with respect to these symmetries. Finally, the effect of $q$ has also
been assessed for $n=3$, $4$ and $5$, which can be used to explore the effect
of the divergence of the equation of state on the relative locations of the
I-$2$–atic and I-$8$-atic bifurcations, and in general on the performance of
the theory when compared with MC simulations.

In order to proceed we first need to derive the free energy. As usual it is
more convenient to obtain the excess free energy by integrating the excess
pressure, which can be written as $\beta p_{\rm exc}a=\beta pa-\eta$ ($\beta
p_{\rm id}=\rho$ is the ideal pressure). Using the thermodynamic relation
$\displaystyle{\beta p_{\rm exc}a=\eta^{2}\frac{\partial\varphi_{\rm
exc}}{\partial\eta}}$, with $\varphi_{\rm exc}$ the excess part of the
Helmholtz free-energy per particle, we can integrate the expression above with
respect to $\eta$ to obtain

$\displaystyle\varphi_{\rm
exc}[h]=-\log(1-\eta)+\sum_{k=2}^{n}c_{k}[h]\operatorname{{\vphantom{\mathnormal{F}}}_{2}\kern-0.5pt\mathnormal{F}_{1}}(k-1,q;k;\eta)\frac{\eta^{k-1}}{k-1}.$
(5)

The integration constant was set to zero to ensure that the excess free energy
be zero at the low-density limit.
$\operatorname{{\vphantom{\mathnormal{F}}}_{2}\kern-0.5pt\mathnormal{F}_{1}}$
is the hypergeometric function,

$\displaystyle\operatorname{{\vphantom{\mathnormal{F}}}_{2}\kern-0.5pt\mathnormal{F}_{1}}\left(n,q;n+1;\eta\right)\frac{\eta^{n}}{n}=\int_{0}^{\eta}\frac{u^{n-1}}{(1-u)^{q}}du.$
(6)

The ideal part of the free-energy per particle has the exact form

$\displaystyle\varphi_{\rm id}[h]=\log\eta-1+\int_{0}^{2\pi}d\phi
h(\phi)\log\left(2\pi h(\phi)\right),$ (7)

giving the total free-energy functional as $\varphi[h]=\varphi_{\rm
id}[h]+\varphi_{\rm exc}[h]$.

Since we will only be interested in assessing the contribution of higher-order
correlations to the stability of the I phase against $p$-atic symmetries, we
restrict here to a stability analysis of the functional against orientational
fluctuations of a given symmetry. In practice this means that the full Fourier
expansion of $h(\phi)$ can be truncated, keeping only the term with the
required symmetry. More specifically, close to the I-$2m$-atic ($m=1,\dots,4$)
bifurcation point, we approximate

$\displaystyle h(\phi)\simeq\frac{1}{2\pi}\left[1+h_{m}\cos(2m\phi)\right],$
(8)

with $h_{m}$ the first Fourier amplitude with $m$-th symmetry. Explicitely,
$m=1$ (2$-$atic), $2$ (4$-$atic), $3$ (6$-$atic) and $4$ (8$-$atic). Inserting
this expression in the total free energy functional $\varphi[h]$ we obtain, to
lowest order in $h_{m}$,

$\displaystyle\varphi[h]=\varphi_{\rm I}+\Delta\varphi[h],$ (9)

where

$\displaystyle\varphi_{\rm
I}=\log\left(\frac{\eta}{1-\eta}\right)-1+\sum_{k=2}^{n}{\cal
C}_{k}^{(0)}[h]\frac{\eta^{k-1}}{k-1}\operatorname{{\vphantom{\mathnormal{F}}}_{2}\kern-0.5pt\mathnormal{F}_{1}}(k-1,q;k;\eta),$
(10)

is the free energy of the I phase, and

$\displaystyle\Delta\varphi[h]=\frac{h_{m}^{2}}{2}\chi_{m}(\eta),\ m\geq 1$
(11)

with

$\displaystyle\chi_{m}(\eta)=1+2\sum_{k=2}^{n}\frac{{\cal
C}_{k}^{(m)}\eta^{k-1}}{(k-1)!}\operatorname{{\vphantom{\mathnormal{F}}}_{2}\kern-0.5pt\mathnormal{F}_{1}}(k-1,q;k;\eta),\quad
m=1,\dots,4,$ (12)

is the extra contribution associated to an orientational fluctuation with
symmetry of order $m$ and amplitude $h_{m}$. In the expressions above we have
defined

$\displaystyle{\cal C}_{n}^{(m)}={\cal
B}_{n}^{(m)}+\sum_{k=1}^{n-2}(-1)^{k}\binom{q}{k}{\cal B}_{n-k}^{(m)},$ (13)

while the coefficients ${\cal B}_{n}^{(m)}$ denote the Fourier components of
the scaled virial coefficients. To quadratic order

$\displaystyle b_{n}[h]={\cal B}_{n}^{(0)}+\frac{h_{m}^{2}}{2}{\cal
B}_{n}^{(m)}.$ (14)

The procedures used to calculate the ${\cal B}_{n}^{(m)}$ coefficients are
described in Appendix A, while the numerical values for the bifurcations are
discussed in Sec. III. We only note here that, in order to find the value of
packing fraction $\eta_{m}$ for the I-$2m$-atic bifurcation we need to solve
the equation $\chi_{m}(\eta)=0$.

In this work we limit the maximum order of the virial coefficient to $m=4$.
However, it will be worthwhile to explore the effect of an extrapolated fifth
and possibly higher-order virial coefficients. From Eqn. (1) with $n=4$, i.e.
the $B_{4}$-SPT approximation, we can obtain the relation between $b_{m}[h]$
for $m\geq 5$ and $b_{2}[h]$, $b_{3}[h]$ and $b_{4}[h]$:

$\displaystyle b_{m}[h]$ $\displaystyle=$
$\displaystyle\left\\{2(m-2)(m-3)b_{4}[h]-2(m-2)(m-4)(q-1)b_{3}[h]\right.$
(15) $\displaystyle+$
$\displaystyle\left.(m-3)(m-4)(q-1)(q-2)b_{2}[h]\right\\}\frac{q(q+1)\cdots(q+m-5)}{2(m-2)!}.$

##  III Results

We start by comparing the values of packing fraction at the bifurcations from
the I phase to the different $p$-atic phases, using extended SPT theories with
increasing number of virial coefficients. Table 1 presents the results. Except
$B_{5}^{*}$-SPT (which is based on an extrapolated fifth virial coefficient),
all of the theories predict a first bifurcation to the $2$-atic phase.
However, the next bifurcation is invariably the I$-8$-atic bifurcation, which
shows the tendency of the HRT fluid to develop octatic correlations. Also, the
difference $\Delta\eta$ in $\eta_{n}$ between the I$-8$-atic and I$-2$-atic
bifurcations tends to dramatically decrease when the theory is extended with
the fourth virial coefficient (cf. $\Delta\eta=0.1195$ for the standard SPT
theory, based on the second virial coefficient, with $\Delta\eta=0.0118$ for
$B_{4}$-SPT, but with the difference between the $B_{3}$-SPT and SPT theories,
$\Delta\eta=0.1028$, hardly changing). This clearly demonstrates the
importance of four-particle correlations in this system. These correlations
are involved in the formation of very stable tetramers of triangles with
square shape, which may dominate the properties of the system at high
densities.

The HRT fluid has been examined in detail by Gantapara et al. Gan  using MC
simulation. This work predicts a first-order phase transition from the I phase
to a 4-atic liquid-crystal phase. Our own simulations MAR1  indicate that the
nature of the liquid-crystal phase strongly depends on the protocol (either
compression or expansion) and starting configurations used in the simulations,
to the extent that the equilibrium configurations may exhibit strong octatic
or purely tetratic correlations. The value of packing fraction at which the I
phase changes to the liquid-crystal phase was obtained to be $\eta=0.733$.
From Table 1 we can see that the $B_{4}$-SPT approximation is very good at
predicting the correct density. Whether the symmetry of the liquid-crystal
phase is tetratic or octatic is a more delicate question that demands further
analysis. Note that the predicted values for $\eta_{2}\sim 1$ are large in all
cases, due to the small values of the ${\cal B}_{n}^{(2)}$ coefficients. This
may indicate that the I-T transition is of first order, as shown in Appendix B
in the framework of SPT.

Bifurcation | I-N | I-T | I-TR | I-O
---|---|---|---|---
$\eta_{n}$ from SPT | 0.8249 | 0.9928 | 0.9821 | 0.9444
$\eta_{n}$ from $B_{3}$-SPT | 0.7325 | 0.9794 | 0.9328 | 0.8353
$\eta_{n}$ from $B_{4}$-SPT | 0.7281 | 0.9681 | 0.8631 | 0.7399
$\eta_{n}$ from $B_{5}^{*}$-SPT | 0.7255 | 0.9590 | 0.8304 | 0.7091
$\eta_{n}$ from $B_{4}$ | 1.0038 | 3.3121 | 1.8485 | 1.2049
$\eta_{n}$ from $B_{5}^{*}$ | 0.8715 | 2.1817 | 1.3071 | 0.9326
Table 1: Values of the packing fractions $\eta_{n}$ at I-N ($n=1$), I-T
($n=2$), I-TR ($n=3$) and I-O ($n=4$) bifurcations from the SPT, $B_{3}$-SPT
and $B_{4}$-SPT theories, and from the truncated virial expansions $B_{4}$ and
$B_{5}^{*}$.

To complete the picture, we have also calculated the bifurcations using
truncated virial expansions instead of the resummed virial series based on the
SPT approximation. The results are shown in Table 1. In common with the
standard Onsager-like theories for low aspect-ratio particles in two
dimensions, the values of packing fractions are absurdly large and generally
unphysical for the theories based on $B_{2}$ and $B_{3}$ virial coefficients
(not shown). But this is because density correlations are too grossly
represented while angular correlations are expected to be more faithfully
captured by the virial coefficients. Therefore the trends in packing fraction
values as more virial coefficients are added may be relevant. In line with the
SPT-based theories, the bifurcation of the I$-$8-atic transition becomes
closer to the I$-$2-atic bifurcation as more virial coefficients are added,
again demonstrating the importance of the higher (especially the fourth)
virial coefficients to represent the structure of the HRT fluid.

To further explore the effect of higher-order virial coefficients, we have
obtained extrapolated values for the fifth virial coefficient. Note that, even
though the third- and fourth-order virial coefficients are obtained
numerically with a numerical effort which is acceptable, the fifth virial
coefficient is much harder to obtain (especially their necessary Fourier
projections are highly fluctuating with the relative angle between particles
and require very detailed and costly MC integrations). Therefore we simply
extrapolate the lower-order coefficients. This is shown in Fig. 1, where the
values of ${\cal B}_{n}^{(k)}$, for $n=2,3$ and $4$, and $k=0,\cdots,4$, are
shown. We can see that all the absolute values of the coefficients are
increasing functions of $n$. The figures include parabolic interpolations to
the data. We take the corresponding extrapolations up to $n=5$ and examine the
consequences. Coming back to the bifurcation densities, we see from Table 1
that the trend for $\Delta\eta$ as more virial coefficients are taken into
account continues to be decreasing, actually becoming negative. This indicates
that fifth-order correlations further promote octatic ordering. This is just a
trend, since the extrapolated values for ${\cal B}_{4}^{(k)}$, and the
resulting $B_{5}^{*}$-SPT theory, need not be accurate. The conclusion is the
same if one looks at the truncated virial series results, Table 1, with
$\eta_{2}<\eta_{4}$ but very close. As a side comment, it is interesting that
the first two bifurcations of the truncated virial expansion of the fifth
order are at physical packing fraction values (i.e. below unity), in contrast
to the lower-order expansions.

Figure 1:  Virial Fourier components ${\cal B}_{n}^{(k)}$ for $k=1,4$ (a) (the
inset shows the case $k=0$) and $k=2,3$ (b) together with parabolic
extrapolation to $n=5$.

We now turn to a discussion on the thermodynamics of the different
approximations by examining the equation of state in the I phase. In the
results presented above the exponent $q$ was set to a value of $2$. Now we
consider flexible choices for $q$ and examine the consequences for the
equation of state in the whole range of densities. In all cases we compare
with the MC results of Gantapara et al. Gan  from the compression runs (these
data have been digitized from the original article and presented as a smooth
curve to aid in visualizing the data in comparison with the different
theoretical equations of state)

Figure 2: (a): EOS from the 3th, 4th and 5th virial expansions. (b): EOS from
Eq. (1) selecting $q=2$ and fixing the correct virial coefficients up to
second, third, fourth and fifth orders (the later using its extrapolated
value). (c): EOS from Eq. (1) using the correct $B_{n}$ ($n=2,\dots,4$) and
using different values of $q$.

Fig. 2 presents the equations of state in the I, for the different theories
and also with different values of $q$. Note that, according to the
simulations, the I phase is stable up to $\eta=0.733$, but the compression
results from the simulations could be taken as reproducing a metastable I
phase beyond the liquid-crystal transition. As we have shown previously MAR1 ,
the structure of the fluid in this regime appears to be quite complex, with
different clusters of particles which give rise to sampling problems in the
simulations. Therefore, the comparison beyond this density should be taken
with caution, but we extend the density interval shown in the figures in order
to assess the impact of different diverging behaviours.

We start by looking at Fig. 2(a), which shows the equations of state from
truncated virial expansion. As expected the addition of more virial
coefficients improve the results (note that only the ${\cal B}_{n}^{(0)}$
coefficients are needed in the I phase). But it is remarkable that the
expansion based on $B_{5}^{*}$ is very accurate in the whole density interval
where the I phase is stable. Clearly the addition of the extrapolated fifth
virial coefficient corrects almost completely the equation of state, which
demonstrates that the extrapolation may be accurate.

Let us consider the equations of state as derived from the resummed expansions
based on SPT theory, i.e. selecting $q=2$. This is shown in Figs. 2(b).
Although the pressure is correctly reproduced at moderate densities by all
theories, the results are somewhat disappointing at high density when they are
compared with the simulations. This is otherwise expected, as resummed
theories are known to be relatively accurate in predicting densities of phase
transitions to liquid-crystal phases, but no so in reproducing the equations
of state quantitatively.

Figure 3: Packing fractions $\eta_{1}(q)$ and $\eta_{4}(q)$ resulting from (a)
$B_{3}$-SPT, (b) $B_{4}$-SPT and (c) $B_{5}^{*}$-SPT theories. In all cases
the bifurcation condition $\chi_{n}(\eta)=0$ has been applied.

Finally, we present an investigation which aims to explore the effect of
different diverging behaviours of the equation of state on the pressure and
the bifurcation points. This is controlled by the parameter $q$ introduced in
Eqn. (1). We present some results for the cases $q=0.5$, $0.75$ and $1$. Fig.
2(c) show the equations of state for these cases. Clearly, as $q$ is decreased
the results improve significantly. The best quantitative agreement is obtained
for $q=0.5$.

The changes in the bifurcation points are shown in Fig. 3, which are focused
on the I-T and I-O bifurcations, represented by $\eta_{1}$ and $\eta_{4}$,
respectively. In the different panels we represent the variations of both
parameters as $q$ is changed within a wide interval, for the three resummed
theories incorporating the third, fourth and (extrapolated) fifth virial
coefficients. The crossover between $\eta_{1}$ and $\eta_{4}$ is induced as
$q$ increases, but the crossover value decreases with the order of the theory.
Note that, for the $B_{4}$-SPT theory, the critical value of $q$ is close to
$2$. Note that an optimal equation of state requires values of $q$ which are
not optimal for the bifurcation point, but the inclusion of higher virial
coefficients may improve this situation, as shown by the resummed theory based
on an extrapolated fifth virial coefficient.

##  IV Conclusions

In this work we have examined the effect of the high-order virial coefficients
on the structure and thermodynamics of the HRT fluid. Virial coefficients and
relevant Fourier projections up to fourth order have been evaluated by MC
integration. In addition, the fifth-order virial coefficient has been
extrapolated. When compared to simulation, the fifth-order truncated virial
expansion, which requires the extrapolated value of the zeroth-order Fourier
component ${\cal B}_{5}^{(0)}$, gives a good approximation for the equation of
state of the I phase up to a packing fraction of $\eta\simeq 0.75$. Note that
this is slightly beyond the predicted range of stability for the I phase
before the liquid-crystal phase becomes stable. In addition, it is the lowest-
order truncated virial expansion for which the I-N and I-O bifurcation points
occur at packing fractions below unity. The resummed $B_{n}$-SPT approaches
overestimate the pressure, but packing fractions at bifurcation are correct
when $n>2$. To investigate the impact of different diverging behaviours of the
equation of state, we explored values of the exponent $q$ different from the
standard choice $q=2$. As $q$ is decreased from $1$ to $0.5$ the resulting
equation of state compares quite reasonably with simulation. However, the
values of the bifurcation packing fractions $\eta_{1}$ and $\eta_{4}$ deviate
quite significantly from simulation, and the difference
$\Delta\eta=\eta_{4}-\eta_{1}$ increases as $q$ decreases.

As far as the resummed theories are concerned, the $B_{4}$-SPT theory gives
$\eta_{1}\lesssim\eta_{4}$, while the $B_{5}^{*}$-SPT theory (with an
extrapolated fifth virial coefficient) already predicts $\eta_{1}>\eta_{4}$.
The cross-over point is located between 1 and 1.5, which in turn implies that
the equation of state will be better than in the standard case. An obvious
outcome of this investigation is that minor changes in the form of the
equation of state may have important quantitative consequences. Overall, the
results presented in this article point to the necessity to include virial
coefficients beyond the fourth if the equation of state and the I-O
bifurcation point are to be reproduced correctly. Unfortunately, such an
effort would incur a high computational cost, rendering the approach
impractical. Clearly, an alternative approach is needed to explain the
behaviour of the HRT fluid.

##  Appendix A Fourier components of virial coefficients

In this section we give details on the calculation of the Fourier components
of the virial coefficients $B_{n}[h]$, $n=2,\dots,4$. These components are
needed in the analysis of the bifurcation points from the I fluid to the
$p$–atic phase. The virial coefficients are in general functionals of the
orientational distribution function $h(\phi)$ and can be written in the usual
diagrammatic form as Hoover

(16) (17) (18)

Nodes and bonds label particles and Mayer functions, respectively. The
reference system at node 1 which means that its position and angle are
$({\bm{r}}_{1},\phi_{1})=({\bf 0},0)$. Therefore the nodes labelled as 2, 3
and 4 (in a clockwise direction) have relative coordinates
${\bm{r}}_{2}-{\bm{r}}_{1}\equiv{\bm{r}}$,
${\bm{r}}_{3}-{\bm{r}}_{1}\equiv{\bm{r}}^{\prime}$,
${\bm{r}}_{4}-{\bm{r}}_{1}\equiv{\bm{r}}^{\prime\prime}$,
$\phi_{2}-\phi_{1}\equiv\phi$, $\phi_{3}-\phi_{1}\equiv\phi^{\prime}$,
$\phi_{4}-\phi_{1}\equiv\phi^{\prime\prime}$. The spatial and angular
integrations with respect to ${\bm{r}}_{1}$ and $\phi_{1}$ in $B_{n}[h]$ can
then be performed trivially. In particular, a factor equal to the total area
$A$ cancels out, see Eqs. (16)-(18). Integration over $\phi_{1}$ allows to
define the following angular functions:

$\displaystyle\Psi_{2}(\phi)\equiv\int_{0}^{2\pi}d\phi_{1}h(\phi_{1})h(\phi_{1}+\phi),$
(19)
$\displaystyle\Psi_{3}(\phi,\phi^{\prime})\equiv\int_{0}^{2\pi}d\phi_{1}h(\phi_{1})h(\phi_{1}+\phi)h(\phi_{1}+\phi^{\prime}),$
(20)
$\displaystyle\Psi_{4}(\phi,\phi^{\prime},\phi^{\prime\prime})\equiv\int_{0}^{2\pi}d\phi_{1}h(\phi_{1})h(\phi_{1}+\phi)h(\phi_{1}+\phi^{\prime})h(\phi_{1}+\phi^{\prime\prime}),$
(21)

Close to the I-$2n$-atic bifurcation point we can approximate $h(\phi)$, up to
first order, as

$\displaystyle
h(\phi)\simeq\frac{1}{2\pi}\left[1+h_{n}\cos(2n\phi)\right],\quad
n=1,\dots,4.$ (22)

To lowest order of the angular functions become

$\displaystyle\Psi_{2}(\phi)=\frac{1}{2\pi}\left\\{1+\frac{h_{n}^{2}}{2}\cos(2n\phi)\right\\},$
(23)
$\displaystyle\Psi_{3}(\phi,\phi^{\prime})=\frac{1}{(2\pi)^{2}}\left\\{1+\frac{h_{n}^{2}}{2}\left[\cos(2n\phi)+\cos(2n\phi^{\prime})+\cos(2n\phi)\cos(2n\phi^{\prime})\right]\right\\},$
(24)
$\displaystyle\Psi_{4}(\phi,\phi^{\prime},\phi^{\prime\prime})=\frac{1}{(2\pi)^{3}}\left\\{1+\frac{h_{n}^{2}}{2}\left[\cos(2n\phi)+\cos(2n\phi^{\prime})+\cos(2n\phi^{\prime\prime})\right.\right.$
$\displaystyle\left.\left.+\cos(2n\phi)\cos(2n\phi^{\prime})+\cos(2n\phi)\cos(2n\phi^{\prime\prime})+\cos(2n\phi^{\prime})\cos(2n\phi^{\prime\prime})\right]\right\\}.$
(25)

The spatial integrals involved in $B_{n}[h]$ can be defined to be angular
kernels,

$\displaystyle{\cal K}_{2}(\phi)=-\frac{1}{2}\int d{\bm{r}}f({\bm{r}},\phi),$
(26) $\displaystyle{\cal K}_{3}(\phi,\phi^{\prime})=-\frac{1}{3}\int
d{\bm{r}}\int
d{\bm{r}}^{\prime}f({\bm{r}},\phi)f({\bm{r}}^{\prime},\phi^{\prime})f({\bm{r}}-{\bm{r}}^{\prime},\phi-\phi^{\prime}),$
(27) $\displaystyle{\cal
K}^{(1)}_{4}(\phi,\phi^{\prime},\phi^{\prime\prime})=-\frac{3}{8}\int
d{\bm{r}}\int d{\bm{r}}^{\prime}\int
d{\bm{r}}^{\prime\prime}f({\bm{r}},\phi)f({\bm{r}}^{\prime\prime},\phi^{\prime\prime})f({\bm{r}}-{\bm{r}}^{\prime},\phi-\phi^{\prime})$
$\displaystyle\hskip 199.16928pt\times
f({\bm{r}}^{\prime}-{\bm{r}}^{\prime\prime},\phi^{\prime}-\phi^{\prime\prime}).$
(28) $\displaystyle{\cal
K}^{(2)}_{4}(\phi,\phi^{\prime},\phi^{\prime\prime})=-\frac{3}{4}\int
d{\bm{r}}\int d{\bm{r}}^{\prime}\int
d{\bm{r}}^{\prime\prime}f({\bm{r}},\phi)f({\bm{r}}^{\prime},\phi^{\prime})f({\bm{r}}^{\prime\prime},\phi^{\prime\prime})f({\bm{r}}-{\bm{r}}^{\prime},\phi-\phi^{\prime})$
$\displaystyle\hskip 199.16928pt\times
f({\bm{r}}^{\prime}-{\bm{r}}^{\prime\prime},\phi^{\prime}-\phi^{\prime\prime}).$
(29) $\displaystyle{\cal
K}^{(3)}_{4}(\phi,\phi^{\prime},\phi^{\prime\prime})=-\frac{1}{8}\int
d{\bm{r}}\int d{\bm{r}}^{\prime}\int
d{\bm{r}}^{\prime\prime}f({\bm{r}},\phi)f({\bm{r}}^{\prime},\phi^{\prime})f({\bm{r}}^{\prime\prime},\phi^{\prime\prime})f({\bm{r}}-{\bm{r}}^{\prime},\phi-\phi^{\prime})$
$\displaystyle\hskip 199.16928pt\times
f({\bm{r}}^{\prime}-{\bm{r}}^{\prime\prime},\phi^{\prime}-\phi^{\prime\prime})f({\bm{r}}-{\bm{r}}^{\prime\prime},\phi-\phi^{\prime\prime}),$
(30)

Note that ${\cal K}_{2}(\phi)$ is just half the excluded area. The superindex
$m$ in the definition of ${\cal K}_{4}^{(m)}(\cdots)$ labels the empty square-
diagram ($m=1$), the square-diagram with one diagonal ($m=2$), and the square-
diagram with two diagonals ($m=3$), respectively. To implement the bifurcation
analysis we need the following Fourier components of these kernels:

$\displaystyle{\cal K}_{2,n}=\int_{0}^{2\pi}d\phi\cos(2n\phi){\cal
K}_{2}(\phi),$ (31) $\displaystyle{\cal
K}_{3,n,m}=\int_{0}^{2\pi}d\phi\cos(2n\phi)\int_{0}^{2\pi}d\phi^{\prime}\cos(2m\phi^{\prime}){\cal
K}_{3}(\phi,\phi^{\prime}),$ (32) $\displaystyle{\cal
K}^{(k)}_{4,n,m,l}=\int_{0}^{2\pi}d\phi\cos(2n\phi)\int_{0}^{2\pi}d\phi^{\prime}\cos(2m\phi^{\prime})\int_{0}^{2\pi}d\phi^{\prime\prime}\cos(2l\phi^{\prime\prime}){\cal
K}^{(k)}_{4}(\phi,\phi^{\prime},\phi^{\prime\prime}),$ (33)

Finally, the virial coefficients $B_{n}[h]$ can be approximated, close to the
bifurcation point, by using Eqns. (16)-(18), (23)-(25), and (26)-(30), as

$\displaystyle B_{2}[h]$ $\displaystyle=$
$\displaystyle\prod_{k=1}^{2}\left(\int_{0}^{2\pi}d\phi_{k}h(\phi_{k})\right){\cal
K}_{2}(\phi)=\int_{0}^{2\pi}d\phi\Psi_{2}(\phi){\cal
K}_{2}(\phi)=B_{2}^{(0)}+\frac{h_{n}^{2}}{2}B_{2}^{(n)},$ $\displaystyle
B_{3}[h]$ $\displaystyle=$
$\displaystyle\prod_{k=1}^{3}\left(\int_{0}^{2\pi}d\phi_{k}h(\phi_{k})\right){\cal
K}_{3}(\phi,\phi^{\prime})$ $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}d\phi\int_{0}^{2\pi}d\phi^{\prime}\Psi_{3}(\phi,\phi^{\prime}){\cal
K}_{3}(\phi,\phi^{\prime})=B_{3}^{(0)}+\frac{h_{n}^{2}}{2}B_{3}^{(n)}$
$\displaystyle B_{4}[h]$ $\displaystyle=$
$\displaystyle\prod_{k=1}^{4}\left(\int_{0}^{2\pi}d\phi_{k}h(\phi_{k})\right)\left(\sum_{m=1}^{3}{\cal
K}^{(m)}_{4}(\phi,\phi^{\prime},\phi^{\prime\prime})\right)$ (34)
$\displaystyle=$
$\displaystyle\int_{0}^{2\pi}d\phi\int_{0}^{2\pi}d\phi^{\prime}\int_{0}^{2\pi}d\phi^{\prime\prime}\Psi_{4}(\phi,\phi^{\prime},\phi^{\prime\prime})\left(\sum_{m=1}^{3}{\cal
K}^{(m)}_{4}(\phi,\phi^{\prime},\phi^{\prime\prime})\right)=B_{4}^{(0)}+\frac{h_{n}^{2}}{2}B_{4}^{(n)},$

with

$\displaystyle B_{2}^{(0)}=\frac{{\cal K}_{2,0}}{2\pi},\hskip
11.38092ptB_{2}^{(n)}=\frac{{\cal K}_{2,n}}{2\pi}$ $\displaystyle
B_{3}^{(0)}=\frac{{\cal K}_{3,0,0}}{(2\pi)^{2}},\hskip
11.38092ptB_{3}^{(n)}=\frac{1}{(2\pi)^{2}}\left(2{\cal K}_{3,n,0}+{\cal
K}_{3,n,n}\right)$ $\displaystyle
B_{4}^{(0)}=\frac{1}{(2\pi)^{3}}\sum_{m=1}^{3}{\cal K}^{(m)}_{4,0,0,0},$
$\displaystyle
B_{4}^{(n)}=\frac{1}{(2\pi)^{3}}\left\\{\sum_{m=1}^{2}\left[2{\cal
K}^{(m)}_{4,n,0,0}+{\cal K}^{(m)}_{4,0,n0}+2{\cal K}^{(m)}_{4,n,n,0}+{\cal
K}^{(m)}_{4,n,0,n}\right]+3\left[{\cal K}^{(3)}_{4,n,0,0}+{\cal
K}^{(3)}_{4,n,n,0}\right]\right\\}.$

From these approximations we define the scaled Fourier components of the
virial coefficients:

$\displaystyle{\cal B}_{k}^{(0)}\equiv\frac{B_{k}^{(0)}}{a^{k-1}}-1,\quad{\cal
B}_{k}^{(n)}\equiv\frac{B_{k}^{(n)}}{a^{k-1}},$ (36)

The second-order coefficients ${\cal B}_{2}^{(n)}$ can be calculated
analytically:

$\displaystyle{\cal
B}_{2}^{(n)}=-\frac{8}{\pi(4n^{2}-1)}\cos^{2}\left[\frac{(2n-1)\pi}{8}\right]\cos^{2}\left[\frac{(2n+1)\pi}{8}\right],$
(37)

The remaining coefficients have to be computed numerically. We have used MC
integration for the spatial integrals and Gaussian quadratures for the angular
integrals, using special tricks to deal with the rapidly varying trigonometric
functions of high index. $10^{8}$ configurations were used to evaluate the
spatial integrals. The results are collected in Table 2.

$k$ | ${\cal B}_{k}^{(0)}$ | ${\cal B}_{k}^{(1)}$ | ${\cal B}_{k}^{(2)}$ | ${\cal B}_{k}^{(3)}$ | ${\cal B}_{k}^{(4)}$
---|---|---|---|---|---
$2$ | $1.8552$ | $-0.1061$ | $-0.0036$ | $-0.0091$ | $-0.0294$
$3$ | $4.9158$ | $-0.3597$ | $-0.0148$ | $-0.0516$ | $-0.1662$
$4$ | $8.6307$ | $-0.6296$ | $-0.0336$ | $-0.1876$ | $-0.5897$
Table 2: Values of the coefficients ${\cal B}_{k}^{(n)}$, obtained
analytically for $k=2$ from Enq. (37), and numerically from MC integration and
Gaussian quadrature for $k>2$.

##  Appendix B Relative stability of the 4–atic phase

For those $n$-th virial theories which predict an I-2–atic bifurcation bellow
the I-8–atic one one may wonder which of the following scenarios takes place
at densities above the I-8–atic bifurcation: (i) The 8–atic and 2–atic free-
energy branches cross each other at some density, or (ii) the 2–atic branch
continues to be the lowest one. In order to investigate this point, we have
minimised the B2-SPT functional, considering a subset of Fourier coefficients
$\\{h_{k}\\}$ with $k=4j$:

$\displaystyle h(\phi)=\frac{1}{2\pi}\left[1+\sum_{k=0}^{n_{\rm
max}}h_{k}\cos(2k\phi)\right].$ (38)

This choice gives a distribution $h(\phi)$ with perfect 8–atic symmetry. A
free-energy branch was generated for a density interval starting at the
I-8–atic bifurcation point and up to densities such that the 8–atic order
parameter is $Q_{8}\simeq 0.97$ (we checked that with this condition the
truncated Fourier series still gives correct results). The 2–atic branch was
also calculated up to densities such that $Q_{2}\simeq 0.97$. Fig. 4(a) shows
the free-energy differences
$\Delta\varphi\equiv\varphi_{\rm\alpha}-\varphi_{\rm I}$ ($\alpha=$2,6,8–atic)
between $\alpha$ and I phases calculated from their respective bifurcation
points. As expected, the 8–atic branch is always metastable. However, the
results indicate that the first scenario above can be discarded, as the
difference between the 8–atic and 2–atic free-energy branches is huge (note
that the latter bifurcates from the isotropic at a much lower density). The
situation is even worse in the case of the metastable 6–atic phase, as it
bifurcates at even higher packing fractions (see Table 1). We note that the
4–atic phase also bifurcates from the I phase but, unlike the 2–atic, 8–atic
and 6–atic phases, it does so via a first-order transition. This can be
demonstrated, not via Fourier-amplitude minimization (we were unable to obtain
a metastable 4–atic solution with the proper restrictions over $\\{h_{k}\\}$),
but from a simple one-parameter minimization of the orientation distribution,

$h(\phi)=\dfrac{e^{\lambda\cos(2n\phi)}}{2\pi I_{0}(\lambda)},\hskip
17.07182pt\phi\in[0,2\pi].$ (39)

$I_{0}(x)$ the zeroth-order modified Bessel function of the first kind. The
results are shown in Fig. 4(b): from the bifurcation point (open circle in the
figure), located at a packing fraction $\eta=0.9928$ (see Table 1), an
unstable 4–atic branch departs towards lower densities, with a free energy
higher than that of the I phase. This branch terminates at $\eta\approx 0.955$
(open square in the figure), where the T phase becomes metastable for the
first time, with a high order parameter $Q_{4}$. From this point a second T
branch develops towards higher densities, $\Delta\varphi$ eventually becoming
negative at $\eta\simeq 0.968$, indicating that the 4–atic phase is more
stable than the I phase. This is the usual scenario for a first-order phase
transition. When $\eta$ further increases from this value, the 4–atic free-
energy branch also crosses the 8–atic branch at $\eta\simeq 0.988$. In any
case, the free energy of the 2–atic phase has by far the lowest value, as can
be seen in Fig. 4. The crossing of the 4–atic and 8–atic branches at high
packing fractions is interesting. It can be understood by invoking the limit
value of the scaled second-virial coefficient as $\lambda\to\infty$:

$\displaystyle\tilde{b}_{2}^{(n)}\equiv\lim_{\lambda\to\infty}b_{2}[h]$
$\displaystyle=\frac{1}{2an}\left[\frac{A_{\rm excl}(0)+A_{\rm
excl}(\pi)}{2}+\sum_{k=1}^{n-1}A_{\rm
excl}\left(\frac{k\pi}{n}\right)\right]-1.$

Inserting the known analytic expression for the excluded area $A_{\rm
excl}(\phi)$, we obtain
$\tilde{b}_{2}^{(1)}<\tilde{b}_{2}^{(2)}\lesssim\tilde{b}_{2}^{(4)}<\tilde{b}_{2}^{(3)}$
for the 2–atic ($n=1$), 4–atic ($n=2$), 6–atic ($n=3$) and 8–atic ($n=4$)
symmetries. This explains the reason for the crossing behavior: the double-
averaged excluded area with respect to $h(\phi)$ for the 4–atic symmetry,
although similar in magnitude, is lower than that obtained for the 8–atic
symmetry. Obviously this occurs only at very high densities, when the
orientational order is almost perfect and the above asymptotic expression can
be justified. The fact that the free energy of the 8–atic phase free-energy is
lower than that of the 4–atic phase at lower densities implies that the
opposite behavior is true when the orientational distribution function is less
sharply peaked. $n$th-order virial theories with an I-8–atic bifurcation below
the I-2–atic one, with the free-energy of the former below that of the latter,
are expected to support the second scenario: After the I-8–atic bifurcation at
high packing fractions, the 4–atic energy branch crosses the 8–atic one, and a
$8$-atic–$4$-atic first-order phase transition takes place. However, if the
crossing point is relatively close to the I-8–atic bifurcation, the two-phase
coexistence might involve the I and the $4$-atic phases, as shown by the MC
simulations.

Figure 4: Free-energy differences
$\Delta\varphi\equiv\varphi_{\alpha}-\varphi_{I}$, with $\alpha$= 2–atic (dot-
dashed), 8–atic (solid), 4–atic (dashed) and 6–atic (dotted), as a function of
packing fraction $\eta$, obtained from (a) the Fourier-coefficient method, and
(b) the one-parameter minimizations. In (a) the dot-dashed line with a shorter
step is a simple extrapolation of the $2$–atic branch to higher densities. The
inset in (b) shows the complete $2$–atic branch. Filled circles in (a) and (b)
show the bifurcation points at second order transitions, while the open circle
in (b) indicates the first-order counterpart in the 4–atic branch. The open
square is the location of the first metastable solution with 4–atic symmetry.

###### Acknowledgements.

Financial support from Grant No. PID2021-126307NB-C21 (MCIU/AEI/FEDER,UE) is
acknowledged.

## References

  * (1) P. G. de Gennes, Symp. Faraday Soc. 5, 165 (1971).
  * (2) U. Tkalec and I. Muševič, Soft Matter 9, 8140 (2013).
  * (3) D. Frenkel and B. M. Mulder, Mol. Phys. 55, 1171-1192 (1985).
  * (4) D. Frenkel and E. Eppenga, Phys. Rev. A 31, 1776 (1985).
  * (5) J. A. Cuesta abd D. Frenkel, Phys. Rev. A 42, 2126 (1990).
  * (6) G. Bautista-Carbajal and G. Odriozola, J. Chem. Phys. 140, 204502 (2014).
  * (7) M. Kleman and G. Ryschenkow, J. Chem. Phys. 64, 4136 (1976).
  * (8) H. H. Wensink and R. L. C. Vink, J. Phys.: Condens. Matter 19, 466109 (2007).
  * (9) H. Schlacken, H.-J. Mogel and P. Schiller, Mol. Phys. 93, 777 (1998).
  * (10) Y. Martínez-Ratón, E. Velasco and L. Mederos, J. Chem. Phys. 122, 064903 (2005).
  * (11) R. Wittmann, C. E. Sitta, F. Smallenburg, and H. Lowën, J. Chem. Phys. 147, 134908 (2017).
  * (12) Y. Martínez-Ratón and E. Velasco, Phys. Fluids 34, 037110 (2022).
  * (13) Y. Geigenfeind, S. Rosenzweig, M. Schmidt, and D. de las Heras, J. Chem. Phys. 142, 174701 (2015).
  * (14) D. A. Triplett and K. Fitchhorn, Phys. Rev. E 77, 011707 (2008).
  * (15) C. Avendaño, F. A. Escobedo, Soft Matter 8, 4675-4681 (2012).
  * (16) S. Mizani, P. Gurin, R. Aliabadi, H. Salehi, and S. Varga, J. Chem. Phys. 153, 034501 (2020).
  * (17) J. A. Martínez-González, J. C. Armas-Pérez, J. Quintana-H, J. Stat. Phys. 150, 559 (2013).
  * (18) P. A. Monderkamp, R. S. Windisch, R. Wittmann, and H. Löwen, arXiv:2301.09541 (2023).
  * (19) J. P. Ramírez González and G. Cinacchi, Phys. Rev. E 102, 042903 (2020).
  * (20) K. Zhao, C. Harrison, D. Huse, W. B. Russel, and P. M. Chaikin, Phys. Rev. E 76, 040401(R) (2007).
  * (21) K. Zhao, R. Bruinsma, and T. G. Mason, Proc. Natl. Acad. Sci. USA 108, 2684 (2011).
  * (22) K. Zhao, R. Bruisma, and T. G. Mason, Nat. Communn. 3, 801 (2012).
  * (23) J. A. Anderson, J. Antonaglia, J. A. Millan, M. Engel, and S. C. Glotzer, Phys. Rev. X 7, 021001 (2017).
  * (24) C. H. Mak, Phys. Rev. E 73, 065104R (2006).
  * (25) B. P. Prajwal y F. A. Escobedo, Phys. Rev. Materials 5, 024003 (2021).
  * (26) T. Geigenfeind y D. de las Heras, J. Chem. Phys. 150, 184906 (2019).
  * (27) Y. Martínez-Ratón and E. Velasco, Phys. Rev. E 79, 011711 (2009).
  * (28) , T. Müller, D. de las Heras, I. Rehberg, and K. Huang, Phys. Rev. E 91, 062207 (2015).
  * (29) M. González-Pinto, F. Borondo, Y. Martínez-Ratón, and E. Velasco, Soft Matter 13, 2571 (2017).
  * (30) Z. Hou, Y. Zong, Z. Sun, F. Ye., T. G. Mason, and K. Zhao, Nat. Commun. 11, 2064 (2020).
  * (31) Y. Martínez-Ratón and E. Velasco, Phys. Rev. E 102, 052128 (2020).
  * (32) Y. Martínez-Ratón, E. Velasco and L. Mederos, J. Chem. Phys. 125, 014501 (2006).
  * (33) Y. Martínez-Ratón and E. Velasco, Phys. Rev. E 104, 054132 (2021).
  * (34) A. P. Gantapara, W. Qi and M. Dijkstra, Soft Matter 11, 8684 (2015).
  * (35) P. Padilla and E. Velasco, J. Chem. Phys. 106, 10299 (1997).
  * (36) F. H. Ree and W. G. Hoover, J. Chem. Phys. 40, 939 (1964).


# Efficient nonlinear witnessing of non-absolutely separable states with lossy
detectors

Ayan Patra1, Shiladitya Mal1,2, Aditi Sen(De)1 1 Harish-Chandra Research
Institute and HBNI, Chhatnag Road, Jhunsi, Allahabad - 211019, India
2 Department of Physics and Center for Quantum Frontiers of Research and
Technology (QFort), National Cheng Kung University, Tainan 701, Taiwan

###### Abstract

Entangled states are undoubtedly an integral part of various quantum
information processing tasks. On the other hand, absolutely separable states
which cannot be made entangled under any global unitary operations are useless
from the resource theoretic perspective, and hence identifying non-absolutely
separable states can be an important issue for designing quantum technologies.
Here we report that nonlinear witness operators provide significant
improvements in detecting non-absolutely separable states over their linear
analogs, by invoking examples of states in various dimensions. We also address
the problem of closing detection loophole and find critical efficiency of
detectors above which no fake detection of non-absolutely separable (non-
absolutely positive partial transposed) states is possible.

##  I Introduction

With the advent of quantum information science, shared entanglement
horodecki'09  turns out to be one of the main resource for quantum
technologies nielsen  which include quantum key distribution ekert'91 , dense
coding dc , teleportation tele , clock synchronisation clock , remote state
preparation rsp , measurement-based computation computation  to name a few. In
the theory of entanglement, developing efficient methods for the generation,
detection and preservation of entangled states is an important enterprise
horodecki'09 ; guhne'09 ; bookchap  while finding unprofitable resources like
unentangled or separable states which can be decomposed in terms of a convex
combination of pure product states werner'89  also plays a crucial role. Over
the years, several criteria have been proposed for identifying entangled
states although according to the computational complexity class, the so-called
entanglement-separability problem is NP-hard gurvits'03 ; ghariban'10 . A
prominent mathematical detection method is the partial transposition map, a
necessary and sufficient criterion for a bipartite system having dimensions up
to six, which is through the existence of all positive but not completely
positive maps peres'96 ; horo'96 .

Apart from these theoretical endeavors, entanglement witness (EW) provides a
powerful tool in the domain of identification of entangled states. Since the
set of the separable state is convex and compact, the Hahn Banach theorem
ensures that there exists a witness operator for each entangled states horo'96
; terhal'20 . EWs are Hermitian operators, whose expectation value with
respect to all separable states is non-negative while it gives negative value
for at least one entangled state. Importantly, they provide an efficient
method of detecting entanglement in laboratories via local measurements,
thereby establishing themselves as useful guhne'02 .

As far as production of entanglement is concerned, separable states can be
made entangled by suitable joint unitary operations horodecki'09 . On the
contrary, there exists a class of separable states which cannot be made
entangled by the application of any joint unitary gate, known as absolutely
separable (AS) states or separable states from spectrum knill'03 ; kusz'01 ;
vers'01 ; hilde'07 ; slater'09 ; johnston'13 ; nirman'14 ; aruna'15 ; jivu'15
. In a similar fashion, absolutely positive partial transposed (PPT) states
are introduced. Although beyond qubit-qutrit states, PPT and separability are
not equivalent due to the existence of PPT bound entangled states horostate'97
; horo'99 , interestingly, it was found that absolute-PPT and AS are same in
these dimensions hilde'07 ; johnston'13 . From the perspective of resource
theory, AS states or absolute-PPT state are a kind of free or useless states
and hence detecting them are significant to identify resources. In the recent
past, witnesses for non-AS states have been proposed in a similar spirit of
linear EWs nirman'14 .

In the domain of detecting entanglement, another interesting twist comes from
the possibility of improving linear witness operators by finding the way of
obtaining nonlinear EWs guhne'06 ; guhne'07 ; zhang'07 ; moroder'08 ; koto'10
; nl0 ; nl1 ; nl2 ; ks'19 ; nl3 . It was shown that every EW can be upgraded
by adding nonlinear term(s) and the novel technique enables us to show that
the set of separable states cannot have facets guhne'07 , thereby shed light
in the geometry of quantum states beng'06  (see recent results on the boundary
of the set of AS states jivu'15a ; hmd ).

In the present work, we explore nonlinear improvement of linear witnesses for
determining non-AS states. To this end, some expectation values in quadratic
form have to be subtracted from the linear witness operators by maintaining
their essential properties. We consider two ways of modifying linear witnesses
guhne'06  – one is by subtracting a single while in the other case, terms
corresponding to a full basis set are subtracted. The method is demonstrated
by considering a class of two-qubit, qubit-qudit, and two-qutrit states. In
all these cases, we show that nonlinear witnesses perform better for detecting
non-AS or non-absolutely PPT states in comparison with the linear ones.
Specifically, in higher dimensions, we propose a class of absolutely separable
as well as absolutely PPT states by mixing bound entangled states with white
noise and construct nonlinear entanglement witnesses explicitly to detect non-
absolutely PPT states. Towards obtaining the results, we also provide classes
of global unitary operators which can transform PPT states to the one having
non-positive partial transposition (NPPT).

Finally, we discuss the change in behavior of these nonlinear witnesses for
non-AS states in a more realistic condition i.e., under inefficient detectors.
Specifically, we consider the scenario when the detectors may not click
indicating lost events. Addressing detection loophole in Bell test is an old
problem btdl . In the context of EW, critical detector efficiency above which
no fake detection of entanglement is possible was derived in Ref. dl'07 .
Recently similar conditions for closing detection loophole in the context of
nonlinear and measurement device-independent EWs are reported ks'19 ; mdl'21 .
In the context of identification of non-AS states, we also show that nonlinear
witnesses provide substantial improvements in critical detector efficiencies
than that obtained via linear witness operators.

We organize the paper in the following way. In Sec. II, the mathematical
condition for absolutely separable as well as absolutely PPT states and
nonlinear entanglement witness operators are presented. To demonstrate the
power of nonlinear witnesses, a class of two-qubit non-AS states is considered
in Sec. III while the classes of non-absolutely PPT states are constructed and
their identification methods via nonlinear witnesses are exhibited in Sec. IV.
In the next section (Sec. V), we provide a method to overcome inefficiencies
in detection procedure and finally we conclude in Sec. VI.

##  II Preliminary concepts, definitions and notations

Before going into the main results, we first present the known criteria for
detecting absolutely separable and PPT states. Let us also briefly introduce
the nonlinear witness operators for the set of non-absolutely separable as
well as non-absolutely PPT states.

###  II.1  Criteria for absolutely separable and positive partial transpose
states and their witness operators

A subset of separable states, which can not be made entangled using any global
unitary operation are called absolutely separable states and similarly
absolutely PPT state are those which remains PPT even after applying any
global unitary operation.

If a bipartite state $\rho_{AB}$ in $2\otimes n$ dimension dimexplain  has
eigenvalues $\lambda_{1}$,$\lambda_{2},\cdots,\lambda_{2n}$ in descending
order, the condition for the absolute separability reads as hilde'07

$\lambda_{1}-\lambda_{2n-1}-2\sqrt{\lambda_{2n-2}\lambda_{2n}}\leq 0.$ (1)

Interestingly, note that in $2\otimes n$, the set containing absolutely PPT
states coincides with the set of absolutely separable states johnston'13
which is not true for the partial transposition criteria in the entanglement-
separability paradigm for $n>3$.

On the other hand, in higher dimensions, i.e., in $3\otimes n$, let the
eigenvalues of $\rho_{AB}$ be $\lambda_{1}$,$\lambda_{2},\cdots,\lambda_{3n}$
organized in the descending order. The states are absolutely PPT hilde'07 ,
when they satisfy the following conditions:

$\begin{split}\begin{vmatrix}2\lambda_{3n}&\lambda_{3n-1}-\lambda_{1}&\lambda_{3n-3}-\lambda_{2}\\\
\lambda_{3n-1}-\lambda_{1}&2\lambda_{3n-2}&\lambda_{3n-4}-\lambda_{3}\\\
\lambda_{3n-3}-\lambda_{2}&\lambda_{3n-4}-\lambda_{3}&2\lambda_{3n-5}\end{vmatrix}\geq
0,\end{split}$ (2)

and

$\begin{split}\begin{vmatrix}2\lambda_{3n}&\lambda_{3n-1}-\lambda_{1}&\lambda_{3n-2}-\lambda_{2}\\\
\lambda_{3n-1}-\lambda_{1}&2\lambda_{3n-3}&\lambda_{3n-4}-\lambda_{3}\\\
\lambda_{3n-2}-\lambda_{2}&\lambda_{3n-4}-\lambda_{3}&2\lambda_{3n-5}\end{vmatrix}\geq
0.\end{split}$ (3)

Notice that there exists no simple criteria for absolutely separable states in
higher dimensions. Moreover, with the increase of dimensions, finding all the
eigenvalues for checking the above criteria requires full tomography tomo
which becomes cumbersome. Hence the witness operators can play a crucial role
to detect non-absolutely PPT states, both theoretically and experimentally.

###  II.2 Nonlinear witnesses for non-absolutely PPT states

In functional analysis, a celebrated theorem, known as Hahn Banach separation
theorem, states that if $S_{1}$ and $S_{2}$ be two nonempty, convex disjoint
subsets of a normed linear space $V$, there is a hyperplane that separates
$S_{1}$ and $S_{2}$ func . Moreover, if $A$ and $B$ are two nonempty disjoint
subsets of a normed linear vector space $V$, where one of them, say, $A$ is
convex, there exists a hyperplane, serving as an witness operator which can
separate the entire subset $A$ from any point of $B$. Since separable states
as well as absolutely separable states both form compact and convex sets, a
Hermitian operator, $W$ for which $\mbox{Tr}(\sigma W)\geq 0$ for all
separable states (absolutely separable states), $\sigma$, and $\mbox{Tr}(\rho
W)<0$ for at least one entangled state (non-AS state), $\rho$, is called the
linear witness operator. They are linear witnesses since linear expression of
mean values of $W$ is involved in the definition. Note also that the set of
absolutely PPT states is also compact, and hence one can construct witness
operators which can detect non-absolutely PPT states. Given two EWs, $W_{1}$
and $W_{2}$, $W_{2}$ is said to be finer than $W_{1}$ if it can detect all
states which are also identified by $W_{1}$ while an EW is called optimal if
there is no other witness finer than it oew .

On the other hand, it was shown that the linear witness operators can always
be upgraded according to its capability of identifying non-separable states by
introducing nonlinearity guhne'06 ; guhne'07 . For example, let us consider
states whose entanglement can be detected by NPPT. In this case, non-
separability can be witnessed by $|\phi\rangle\langle\phi|^{T_{B}}$ where
$|\phi\rangle$ is the eigenvector corresponding to the negative eigenvalue of
$\rho^{T_{B}}$. As prescribed in guhne'06 , one can introduce nonlinearity in
the following ways.

$F^{1}(\rho)=\langle|\phi\rangle\langle\phi|^{T_{B}}\rangle-\dfrac{1}{S(\psi)}\langle{X}^{T_{B}}\rangle\langle({X}^{T_{B}})^{\dagger}\rangle$
(4)

and

$F^{2}(\rho)=\langle|\phi\rangle\langle\phi|^{T_{B}}\rangle-\sum_{i=1}^{k}\langle{X}_{i}^{T_{B}}\rangle\langle({X}_{i}^{T_{B}})^{\dagger}\rangle$
(5)

where all the expectation values are taken with respect to the given state
$\rho$. Here in Eq. (4), $X$ is given by $|\phi\rangle\langle\psi|$, where
$|\psi\rangle$ is any arbitrary state and $S(\psi)$ is the square of the
largest Schmidt coefficient of the state $|\psi\rangle$ while in Eq. (5),
$X_{i}$ is defined by $|\phi\rangle\langle\psi_{i}|$ ($i=1,2,\cdots,d$) with
an orthonormal basis being $\\{|\psi_{i}\rangle\\}$.

Let us now adopt the similar procedure to detect non-absolutely separable
(non-absolutely PPT) states, say, $\rho^{\prime}$ which can be made entangled
from a separable (PPT) one by global unitary operator $U$ nirman'14 . In this
situation, to detect $\rho^{\prime}$, the witness operator,
$|\phi\rangle\langle\phi|^{T_{B}}$ should be modified as
$U^{\dagger}|\phi\rangle\langle\phi|^{T_{B}}U$, where $|\phi\rangle$ is the
eigenvector corresponding to the negative eigenvalue of
$({U}\rho^{\prime}{U}^{\dagger})^{T_{B}}$.

Let us now incorporate the nonlinear terms into witness operators for
determining non-AS (non-PPT) states. Like Eqs. (4) and (5), the detection
method for non-AS (non-absolutely PPT) states take the form as

$\displaystyle F^{1}(\rho^{\prime})$ $\displaystyle=$
$\displaystyle\langle{U}^{\dagger}|\phi\rangle\langle\phi|^{T_{B}}U\rangle-\dfrac{1}{S(\psi)}\langle{U}^{\dagger}{X}^{T_{B}}U\rangle\langle(U^{\dagger}{X}^{T_{B}}U)^{\dagger}\rangle,$
$\displaystyle F^{2}(\rho^{\prime})$ $\displaystyle=$
$\displaystyle\langle{U}^{\dagger}|\phi\rangle\langle\phi|^{T_{B}}U\rangle$
(6) $\displaystyle-$
$\displaystyle\sum_{i=1}^{k}\langle{U}^{\dagger}{X_{i}}^{T_{B}}U\rangle\langle(U^{\dagger}{X_{i}}^{T_{B}}U)^{\dagger}\rangle,$

where $X^{T_{B}}$ is replaced by $U^{\dagger}{X}^{T_{B}}U$ and all the
expectation values have to be taken with respect to the given state
$\rho^{\prime}$. As we will show that in the detection process of non-
absolutely separable (PPT) states, finding a nontrivial global $U$ is, in
general, difficult.

Figure 1: (Color online.) Witness operators, $\langle W\rangle$, vs. $p$ of
the generalized Werner state, $\rho_{gW}$ in Eq. (7). Solid (blue), dashed
(red) and dotted (green) lines correspond to nonlinear witness operators,
$F^{2}$, $F^{1}$ and linear witnesses respectively. Here $\alpha=\pi/12$ in
$|\xi\rangle$ and $\theta=0.352$ involved in the witness operator,
$|\phi\rangle$, in radians. Both the axes are dimensionless.

##  III Nonlinear improvement in witnessing two-qubit non-absolutely separable
states

We will now explicitly show advantages of nonlinear witness over their linear
counterparts by constructing them and by invoking a class of two-qubit states.

We now use the criteria, given in (1) to find the range of the state parameter
for which the state is absolutely separable. To illustrate the method, let us
consider a generalized Werner state werner'89 ,

$\rho_{gW}=p|\xi\rangle\langle\xi|+\dfrac{1-p}{4}I_{4\times 4}.$ (7)

where $|\xi\rangle=\cos\alpha|00\rangle+e^{i\phi}\sin\alpha|11\rangle$, $p$ is
the mixing parameter and $I_{4\times 4}$ is the identity operator with the
subscript representing the size of the matrix, thereby indicating the
dimension of the system. Following the condition written in (1) with $2n=4$,
we find that $\rho_{gW}$ is absolutely separable in the range $0\leq
p\leq\frac{1}{3}$ which is independent of $\alpha$ and $\phi$, while the state
is entangled when $\dfrac{1}{1+2\sin 2\alpha}<p\leq 1$. Except $\alpha=\pi/4$,
there exists a range of $p$ in which the state is separable but not absolutely
separable. For example, for $\alpha=\frac{\pi}{12}$, i.e., $\sin
2\alpha=\frac{1}{2}$, the state $\rho_{gW}$ is entangled with
$\dfrac{1}{2}<p\leq 1$ and separable but not absolutely separable in
$1/3<p\leq 1/2$. It implies that there exists a global unitary operator which
can convert the state $\rho_{gW}$ to an entangled state (or, NPPT in $2\otimes
2$ ), say $\rho_{e}$ in the range $\frac{1}{3}<p\leq\frac{1}{1+2\sin
2\alpha}$.

Let us consider a global unitary operator, given by

$\displaystyle U=\dfrac{1}{\sqrt{2}}\begin{pmatrix}1&0&0&1\\\
0&\sqrt{2}&0&0\\\ 0&0&\sqrt{2}&0\\\ -1&0&0&1\end{pmatrix}.$ (8)

After applying this unitary operator, $\rho_{gW}$ becomes $\rho_{e}$ for
$p>0.366$, i.e. when $0.366<p<0.5$, the state is entangled although it was
initially separable, i.e., PPT.

By choosing $|\phi\rangle=\frac{1}{\sqrt{2}}(|01\rangle+|10\rangle)$, we can
show that $U^{\dagger}|\phi\rangle\langle\phi|^{T_{B}}U$ can detect the
$\rho_{gW}$ as non-absolutely separable in the range $0.366<p<0.5$. Let us
take $|\phi\rangle=\cos\theta|01\rangle+\sin\theta|10\rangle$,
$(0\leq\theta\leq\pi)$ and in that case,
$U^{\dagger}|\phi\rangle\langle\phi|^{T_{B}}U$ can detect the state
$\rho_{gW}$ as non-absolutely separable in the range $\dfrac{1}{1+\sqrt{3}\sin
2\theta}<p<0.5$, with the condition being

$\langle W(\rho_{1})\rangle=\frac{-\sqrt{3}p\sin 2\theta-p+1}{4}<0.$ (9)

To make a relatively weaker witness than optimal with $\theta=\pi/4$, consider
any value of $\theta$ between $\pi/10.21$ and $\pi/2.4872$, to detect the
range for which the state $\rho_{gW}$ is non-absolutely separable. Let us now
introduce nonlinearity to improve the range of detection. By considering
$|\psi\rangle=|01\rangle$, and by using Eq. (II.2), we can reach to the
condition for the state $\rho_{gW}$ to be non-absolutely separable if

$\displaystyle F^{1}_{|01\rangle}(\rho_{1})=\frac{-\sqrt{3}p\sin
2\theta-p+1}{4}$ (10)
$\displaystyle-\frac{(\cos\theta(1-p)-\sqrt{3}p\sin\theta)^{2}}{16}<0,$

while using $|\psi\rangle$=$|10\rangle$, the condition gets modified as

$\displaystyle F^{1}_{|10\rangle}(\rho_{1})=\frac{-\sqrt{3}p\sin
2\theta-p+1}{4}$ (11)
$\displaystyle-\frac{(\sin\theta(1-p)-\sqrt{3}p\cos\theta)^{2}}{16}<0.$

Instead of product states, if we use entangled state as
$|\psi\rangle=\frac{1}{\sqrt{2}}(|01\rangle\pm|10\rangle)$, we cannot provide
any advantage except for some range of $\theta$. On the other hand, for
obtaining $F^{2}$, we choose an orthonormal basis,
$\\{|\psi\rangle_{i}=|00\rangle,|01\rangle,|10\rangle,|11\rangle$}, and obtain
the criteria for witnessing non-absolutely separable state as

$\displaystyle F^{2}(\rho_{1})=\dfrac{-\sqrt{3}p\sin 2\theta-p+1}{4}$ (12)
$\displaystyle-\dfrac{(\cos\theta(1-p)-\sqrt{3}p\sin\theta)^{2}}{16}$
$\displaystyle-\dfrac{(\sin\theta(1-p)-\sqrt{3}p\cos\theta)^{2}}{16}<0.$

Notice that instead of the computational basis, $\\{|\psi_{i}\rangle\\}$, if
we take the Bell basis (i.e., entangled states as basis elements), we reach to
the same condition as above. Fig. 1 depicts the comparison between linear and
two kinds of nonlinear witness operators in case of detecting non-absolutely
separable state, $\rho_{gW}$ for fixed values of $\alpha$ and $\theta$.

##  IV Construction of Non-absolutely PPT states with the conjunction of PPT
bound entangled states and their witnesses

In this section, we show the usefulness of nonlinearity in witness operators
for recognizing non-absolutely separable as well as PPT states in higher
dimensions. To illustrate this, we construct absolutely separable (PPT) states
by mixing white noise with PPT bound entangled states in $2\otimes 4$ and
$3\otimes 3$. Such examples also shed light on the boundaries of the set of
PPT bound entangled states and absolutely PPT states.

Figure 2: (Color online.) Mixing parameter, $p$ against $b$ in
$\rho_{2}=p\rho_{b}+\dfrac{1-p}{8}I_{8\times 8}$ in Eq. (13). Circles (red),
squares (blue) and triangles (green) represent the boundaries of absolutely
separable states and non-AS states, PPT states and NPPT states after applying
unitary operators $U(\pi/3,\pi)$ and $U_{1}$ in Eqs. (15) and (23
respectively. Both the axes are dimensionless.

####  IV.0.1 Witnessing non-absolutely separable states via PPT bound
entangled states

In $2\otimes 4$ system, let us concentrate on the mixture of the PPT bound
entangled state with white noise, given by

$\rho_{2}=p\rho_{b}+\dfrac{1-p}{8}I_{8\times 8},$ (13)

where $\rho_{b}$ is the bound entangled state horostate'97  represented as

$\rho_{b}=\dfrac{1}{7b+1}\begin{pmatrix}b&0&0&0&0&b&0&0\\\ 0&b&0&0&0&0&b&0\\\
0&0&b&0&0&0&0&b\\\ 0&0&0&b&0&0&0&0\\\
0&0&0&0&\dfrac{1+b}{2}&0&0&\dfrac{\sqrt{1-b^{2}}}{2}\\\ b&0&0&0&0&b&0&0\\\
0&b&0&0&0&0&b&0\\\
0&0&b&0&\dfrac{\sqrt{1-b^{2}}}{2}&0&0&\dfrac{1+b}{2}\end{pmatrix},$ (14)

with $b\in[0,1]$. Except for $b=0,1$ where $\rho_{b}$ is separable, it is PPT
bound entangled. The condition for this state to be absolutely separable,
thereby absolutely PPT in this case is given by (following inequality (1) with
$2n=8$), $\lambda_{1}-\lambda_{7}-2\sqrt{\lambda_{6}\lambda_{8}}\leq 0$, where
$\lambda_{i}$s are eigenvalues of $\rho_{2}$ in descending order. In Fig. 2,
circles indicates the boundary of absolutely separable and non-AS states in
the $(b,p)$-plane and hence the region below the boundary represents the
absolute separability or absolute PPT of $\rho_{2}$. Moreover, we find that
$\rho_{2}$ is always PPT for any values of b and p and therefore, it is
important to identify the states lying above the boundary are non-AS states
which can be converted to NPPT states with the help of global unitary
operators. For example, we consider a global unitary operator, $U$, given by

$\displaystyle U(\phi_{1},\phi_{2})$ $\displaystyle=$ $\displaystyle
a_{1}\;[\sigma_{x}\otimes\sigma_{y}\otimes\sigma_{z}]+a_{2}\;[\sigma_{y}\otimes\sigma_{z}\otimes\sigma_{x}]$
(15) $\displaystyle+$ $\displaystyle
a_{3}\;[\sigma_{z}\otimes\sigma_{x}\otimes\sigma_{y}]$

where $\sigma$’s are Pauli spin matrices with some parameters $a_{1}$, $a_{2}$
and $a_{3}$, such that $a_{1}^{2}+a_{2}^{2}+a_{3}^{2}=1$. We parametrize them
as $a_{1}=\cos\phi_{1},\,a_{2}=\sin\phi_{1}\sin\phi_{2}\text{ and
}a_{3}=\sin\phi_{1}\cos\phi_{2}$, $0\leq\phi_{1}\leq\pi,\,0\leq\phi_{2}\leq
2\pi$. Taking $\phi_{1}$ and $\phi_{2}$ as $\pi/3$ and $\pi$ respectively, we
observe that some non-AS states become NPPT which are marked in Fig. 2 by
squares. We observe that there still exists some non-absolute PPT states
(lying between the envelopes of circles and squares) which we cannot make NPPT
by this unitary operator.

Figure 3: Witness operators with $p$ for $\rho_{2}$ in Eq. (13). Choices of
$\\{\theta_{i}\\}$s in $|\phi^{\prime}\rangle$ are mentioned in the text. Here
$b=0.7$. All other specifications are same as in Fig. 1. Both the axes are
dimensionless.

At this point, it should be noted that different values of $b$ require
different witnesses to detect the state properly. Let us illustrate the entire
process for $b=0.7$. We find that
$W=U(\pi/3,\pi)^{\dagger}|\phi\rangle\langle\phi|^{T_{B}}U(\pi/3,\pi)$ can
detect $\rho_{2}$ as non-AS state, when $0.62<p\leq 1$ as shown in Fig. 3.
Here $|\phi\rangle$ is the column vector
$(-0.13348,\;0.67743,\;-0.09738,\;0.02271,\;0.00333,\;0.04054,\\\
-0.71427,\;0.03788)^{T}$ . It is the best linear witness operator for $b=0.7$
and for a given $U(\pi/3,\pi)$, since $W$ constructed by using $|\phi\rangle$
can detect all the non-AS states when $p\in(0.62,1]$.

Let us now take a general $|\phi^{\prime}\rangle$ =
$(\sin\theta_{1}\sin\theta_{2}\sin\theta_{3}\sin\theta_{4}\sin\theta_{5}\sin\theta_{6}\sin\theta_{7},\\\
\sin\theta_{1}\sin\theta_{2}\sin\theta_{3}\sin\theta_{4}\sin\theta_{5}\sin\theta_{6}\cos\theta_{7},\;\sin\theta_{1}\sin\theta_{2}\\\
\sin\theta_{3}\sin\theta_{4}\sin\theta_{5}\cos\theta_{6},\;\sin\theta_{1}\sin\theta_{2}\sin\theta_{3}\sin\theta_{4}\cos\theta_{5},\\\
\sin\theta_{1}\sin\theta_{2}\sin\theta_{3}\cos\theta_{4},\;\sin\theta_{1}\sin\theta_{2}\cos\theta_{3},\;\sin\theta_{1}\cos\theta_{2},\\\
\cos\theta_{1})^{T}$ with $0\leq\theta_{i}\leq\pi,\,\,i=1,\ldots 6$ and
$0\leq\theta_{7}\leq 2\pi$. In case of linear witness operator,
$U(\pi/3,\pi)^{\dagger}|\phi^{\prime}\rangle\langle\phi^{\prime}|^{T_{B}}U(\pi/3,\pi)$
can also detect $\rho_{2}$ as non-AS state in some range of $p$ depending on
the parameter values involved in $|\phi\rangle$. However, we can introduce
nonlinearity to improve the range of detection. For this purpose, we consider
$|\psi\rangle$=$\frac{1}{\sqrt{2}}(|00\rangle+|10\rangle)$ for $F^{1}$-type of
nonlinear witness operator and an orthonormal basis
$\\{|00\rangle,|01\rangle,|02\rangle,|03\rangle,|10\rangle,|11\rangle,|12\rangle,|13\rangle\\}$
for $F^{2}$. Now, using Eq. (II.2), and choosing
$\theta_{1}=2.07345,\;\theta_{2}=2.36710,\;\theta_{3}=1.5128,\;\theta_{4}=1.508,\;\theta_{5}=1.5382,\;\theta_{6}=1.7109,\;\theta_{7}=0.19455$
in $|\phi\rangle$ (where all values are in radians), we observe a clear
improvement over linear witness operators (see Fig. 3).

_Remark 1._ The similar method can also be applied for other values of $b$.

_Remark 2._ For the clear demonstration of the utility of nonlinear witness
operator, we choose a set of $\\{\theta_{i}\\}$. Different values of
$\\{\theta_{i}\\}$ lead to a qualitatively similar result.

Figure 4: The mixing parameter, $p$, with respect to $b$ in the state,
$\rho_{3}$ in Eq. (16). Here we choose $U(\pi/18,5\pi/6)$ and $U_{1}$ as given
in Eqs. (18) and (24). All other specifications are same as in Fig. 2.

####  IV.0.2 Detecting non-absolute PPT states in $3\otimes 3$

Let us move to an example of a class of two-qutrit state. This example is
different than the examples considered before since in $3\otimes 3$, the set
of absolutely separable states and absolutely PPT states are different and we
concentrate on the detection of non-absolutely PPT states. In this purpose,
let us consider a state

$\rho_{3}=p\rho^{\prime}_{b}+\dfrac{1-p}{9}I_{9\times 9},$ (16)

where $\rho_{b}$ horo'99  is given by

$\rho^{\prime}_{b}=\dfrac{2}{7}|\tilde{\psi}\rangle\langle\tilde{\psi}|+\dfrac{b}{7}\sigma_{+}+\dfrac{5-b}{7}\sigma_{-}.$
(17)

Here
$|\tilde{\psi}\rangle=\dfrac{1}{\sqrt{3}}(|00\rangle+|11\rangle+|22\rangle)$,
$\sigma_{+}=\dfrac{1}{3}(|01\rangle\langle 01|+|12\rangle\langle
12|+|20\rangle\langle 20|)$, and $\sigma_{-}=\dfrac{1}{3}(|10\rangle\langle
10|+|21\rangle\langle 21|+|02\rangle\langle 02|)$. The state
$\rho^{\prime}_{b}$ is PPT for $1\leq{b}\leq 4$, and so we confine in this
range of $b$. By using inequalities, (2) and (3), we provide the range of $b$
and $p$ for which the state $\rho_{3}$ is absolutely PPT (see the red circles
in Fig. 4 for the boundary of absolutely PPT and non-absolutely PPT states).
Like Eq. (15), the unitary operator in this case reads as

$\displaystyle U(\phi_{1},\phi_{2})$ $\displaystyle=$
$\displaystyle\cos\phi_{1}[\sigma_{x}\otimes\sigma_{y}\otimes\sigma_{z}]$ (18)
$\displaystyle+$
$\displaystyle\sin\phi_{1}\sin\phi_{2}[\sigma_{y}\otimes\sigma_{z}\otimes\sigma_{x}]$
$\displaystyle+$
$\displaystyle\sin\phi_{1}\cos\phi_{2}[\sigma_{z}\otimes\sigma_{x}\otimes\sigma_{y}]\oplus[1],$

where $[1]$ is a $1\times 1$ matrix with entry $1$. A unitary operator of this
class, specifically, $U(\pi/18,5\pi/6)$, is capable to make non-absolutely PPT
states to NPPT for some region in the $(b,p)$-plane as depicted by squares in
Fig. 4. Notice that we require different unitary operator if the non-
absolutely PPT states belonging to the region between circles and squares have
to make NPPT (as also seen in Fig. 2).

Like in the previous example, we fix $b=1.5$. Let us first construct a linear
witness operator,
$W=U(\pi/18,5\pi/6)^{\dagger}|\phi\rangle\langle\phi|^{T_{B}}U(\pi/18,5\pi/6)$,
which detects the state $\rho_{3}$ as non-absolutely PPT for $0.6<p\leq 1$. In
this case,
$|\phi\rangle=(-0.4476-0.004054i,\;-0.0103-0.009966i,\;-0.001158+0.3953i,\;0.02944-0.04832i,\;0.0003527-0.001285i,\;-0.05052+0.01027i,\;0.000449-0.3918i,\;-0.0478-0.03061i,\;0.6933)^{T}$.

To portray the power of nonlinear witness operator, let us choose
$|\phi^{\prime}\rangle=(p_{1}+ip_{2},\;-0.0103-0.009966i,\;-0.001158+0.3953i,\;0.02944-0.04832i,\;0.0003527-0.001285i,\;-0.05052+0.01027i,\;0.000449-0.3918i,\;-0.0478-0.03061i,\;p_{3}+ip_{4})^{T}$,
with $p_{1}=-0.564882,\;p_{2}=0.471498,\;p_{3}=0.373546,\;p_{4}=0.0$, and for
$F^{1}$ and $F^{2}$, $|\psi\rangle=\frac{1+i}{\sqrt{2}}|22\rangle$ and the
orthonormal basis as
$\frac{1+i}{\sqrt{2}}\\{|00\rangle,|01\rangle,|02\rangle,|10\rangle,|11\rangle,|12\rangle,|20\rangle,|21\rangle,|22\rangle\\}$
respectively. Notice that $|\phi^{\prime}\rangle$ is almost same as
$|\phi\rangle$ in the linear witness operator except the first and the last
entries. With these parameter values, we again identify a range of the noise
value, $p$, for which the performance of $F^{2}$ is better than $F^{1}$ and
the linear witness operator (see Fig. 5).

Figure 5: Plots of nonlinear and linear witness operators, $\langle W\rangle$,
against the noise parameter $p$ for $\rho_{3}$ in Eq. (16). Here $b=1.5$. The
choices of $\\{p_{i}\\}$s and basis required for witness operators are
mentioned in the text. All other specifications are same as in Fig. 1.

##  V DETECTION LOOPHOLE IN NONLINEAR WITNESS operator for non-absolutely
separable (PPT) states

Upto now, all the results derived in Secs. III and IV are under the assumption
that the apparatuses used in the identification process are perfect. We will
now investigate the response of imperfect detectors like no clicking of
detectors (lost events) or additional events on nonlinear witnesses dl'07 ;
ks'19 . We find critical efficiency of detectors so as to avoid the fake
detection of resources.

Towards fulfilling the aim, we decompose a witness operator in local operator
basis such that $W=C_{0}I+\sum_{\alpha}C_{\alpha}{S}_{\alpha}$, where I stands
for the identity operator, $S_{\alpha}$ represents a term from expansion of
$W$ in local operator basis and $C_{\alpha}$ is the corresponding expansion
coefficient dl'07 . Using this decomposition, it can be shown that the
condition for detecting a NPPT state using linear witness operator
experimentally when detector does not work ideally is given by

$\langle{W}_{\phi}\rangle_{m}<C_{0}(1-\frac{1}{\eta_{-}})$ (19)

where $\langle{W}_{\phi}\rangle_{m}$ is experimentally measured value of the
witness operator and $\eta_{-}$ is the lost event efficiency of detector
(which is supplied). We here assume that the additional event efficiency
vanishes. On the other hand, in case of non-linear witness operator, $F^{1}$,
the relation between the measured and the true values modifies as ks'19

$\begin{split}\langle{W}_{\phi}\rangle_{m}<C_{0}(1-\frac{1}{\eta_{-}})+\frac{\eta_{-}}{S(\psi)}(\langle{H}\rangle_{m}^{2}+{K}_{H}^{2}-2\langle{H}\rangle_{m}K_{H})\\\
+\frac{\eta_{-}}{S(\psi)}(\langle{A}\rangle_{m}^{2}+{K}_{A}^{2}-2\langle{A}\rangle_{m}K_{A}),\end{split}$
(20)

where $\langle{H}\rangle_{m}$ and $\langle{A}\rangle_{m}$ are experimentally
measured value of Hermitian and anti-Hermitian part of $X^{T_{B}}$
respectively. Moreover $K_{H}$ and $K_{A}$ are given by

$K_{H}=C_{0H}(1-\frac{1}{\eta_{-}})\text{ and
}K_{A}=C_{0A}(1-\frac{1}{\eta_{-}}),$ (21)

with $C_{0H}$ and $C_{0A}$ being the coefficients of the decomposition of $H$
and $A$ in local operator basis respectively, corresponding to the Identity
matrix.

It is possible to obtain $\langle{W}_{\phi}\rangle_{m}$ for all the examples
that we considered so far. For example, when $\eta_{-}$ is supplied, for the
generalized two-qubits Werner state,
$\langle{W}_{\phi}\rangle_{m}<\frac{1}{4}(1-\frac{1}{\eta_{-}})$ with linear
witness operator while for nonlinear case,
$\langle{W}_{\phi}\rangle_{m}<\frac{1}{4}(1-\frac{1}{\eta_{-}})+\eta_{-}\\{(\langle{H}\rangle_{m}-K_{H})^{2}+\langle{A}\rangle_{m}^{2}\\}$
with $K_{H}=\frac{1}{4}\cos\theta(1-\frac{1}{\eta_{-}})$ and $K_{A}=0$. So,
for two-qubits, the loophole cannot be closed for our chosen linear witness
(with $\theta=0.352$ ) if $\eta_{-}<0.43582$ while in case of nonlinear
witness operator, for a given $X_{nl}=0.4$, the upper bound of $\eta_{-}$ can
further be lowered as $0.3881$ where we assume that $|\psi\rangle$ is
orthogonal to $|\phi\rangle$ in $F^{1}$.

Considering a qubit-qudit state $\rho_{2}$ given in Eq. (13), let us
demonstrate the advantageous role of nonlinear witness operators towards
defeating the inefficiency in detectors. In this case, the loophole of a
linear witness can be closed when
$\langle{W}_{\phi}\rangle_{m}<\frac{1}{8}(1-\frac{1}{\eta_{-}})$, while the
similar condition in presence of nonlinearity reads as
$\langle{W}_{\phi}\rangle_{m}<\frac{1}{8}(1-\frac{1}{\eta_{-}})+\eta_{-}\\{(\langle{H}\rangle_{m}-K_{H})^{2}+\langle{A}\rangle_{m}^{2}\\}$,
where
$K_{H}=\frac{1}{8\sqrt{2}}(1-\frac{1}{\eta_{-}})\sin\theta_{1}\sin\theta_{2}\sin\theta_{3}(\cos\theta_{4}+\sin\theta_{4}\sin\theta_{5}\sin\theta_{6}\sin\theta_{7})$
and $K_{A}=0$.

Let us consider a special case when $|\psi\rangle$ is orthogonal to
$|\phi\rangle$, thereby leading to the vanishing $K_{H}$ and $K_{A}$. For a
fixed $\eta_{-}$, the lost events inefficnecy can be overcome when

$\langle{W}_{\phi}\rangle_{m}<\frac{1}{8}(1-\frac{1}{\eta_{-}})+\frac{\eta_{-}}{S(\psi)}\\{\langle{H}\rangle_{m}^{2}+\langle{A}\rangle_{m}^{2}\\}\equiv
W^{up}.$ (22)

The above equation can be rewritten if we consider
$(\langle{H}\rangle_{m}^{2}+\langle{A}\rangle_{m}^{2})=X_{nl}^{2}$. In Fig. 6,
the upper bound of $\langle{W}_{\phi}\rangle_{m}$ is shown with respect to
$X_{nl}$ for different values of $\eta_{-}$ which indicates that the chances
of detecting non-absolutely separable states (with respect to the detector
efficiency) increase with the increasing value of nonlinearity, $X_{nl}^{2}$.
To visualize it, let us consider the situation when
$\langle{W}_{\phi}\rangle_{m}$ vanishes and the nonlinear term is taken to be
$X_{nl}=0.2$. The detection of non-AS state is then possible when
$\eta_{-}\geq 0.424$, while for $X_{nl}=0.6$, $\eta_{-}\geq 0.275$. Therefore,
the increase of nonlinearity, i.e., $X_{nl}$ enhances the possibility of
identifying non-AS states even in presence of a relatively inferior detector.
If we compare the values with the linear witness operators, we can also show
that nonlinearity in the witness operator helps to overcome the detection
loopholes.

The similar detection inefficiency can also be overcome when the task is to
identify the range of $p$ in $\rho_{3}$ representing the non-absolute PPT
states. In this case, $K_{H}$ and $K_{A}$ can be evaluated to be
$\frac{1}{9\sqrt{2}}(p_{3}+p_{4})(1-\frac{1}{\eta_{-}})$ and
$\frac{1}{9\sqrt{2}}(p_{3}-p_{4})(1-\frac{1}{\eta_{-}})$ respectively and
hence the condition in (20) can also be obtained for a fixed $\eta_{-}$.

Figure 6: Inefficiency in detectors. The upper bound of $\langle
W_{\phi}\rangle_{m}$ in the inequality (22), $W^{up}$, against
$X_{nl}=(\langle{H}\rangle_{m}^{2}+\langle{A}\rangle_{m}^{2})$ for $\rho_{2}$.
We plot for different values of $\eta_{-}$. $|\phi\rangle$ chosen for
computing $F_{1}$ and $U$ are same as in Fig. 3. The behavior of $W^{up}$
indicates that and with the increase of lost events inefficiency, $\eta_{-}$,
nonlinearity helps to obtain the high measured witness value. Both the axes
are dimensionless.

##  VI Conclusion

Quantum information processing tasks can only successfully be realized if the
resource states required for that particular job are prepared and identified
in an efficient manner. Among several available resources, entanglement shared
between multiple parties has become one of the important ingredients in most
of the quantum protocols discovered to date. Interestingly, however, it was
found that there exists a set of unentangled states which can be made
entangled by global unitary transformation while the rest of the states
remains useless, known as absolutely separable states.

Therefore, in the development of quantum technologies, determining non-
absolutely separable states in laboratories can be a significant issue. Among
several identification methods developed in the theory of entanglement, the
most experimental-friendly one is the linear witness operators although, for a
given state, the general method of obtaining an optimal witness operator is
still not known. On the other hand, nonlinear witness operators are shown to
be good alternatives with respect to the detection of useful resources.

The present work develops nonlinear witness operators for detecting non-
absolutely separable as well as non-absolutely positive partial transposed
(PPT) states in presence of an ideal and inefficient detector. We explicitly
constructed nonlinear witness operators for the class of two-qubit states and
admixtures of bound entangled states with white noise in different dimensions.
Specifically, we showed that in presence of noise, nonlinear witness operators
can be more efficient to recognize non-absolutely PPT states compared to their
linear counterparts. Moreover, we found that such advantages are more
pronounced when the detectors are inefficient.

## Acknowledgement

We acknowledge the support from Interdisciplinary Cyber Physical Systems
(ICPS) program of the Department of Science and Technology (DST), India, Grant
No.: DST/ICPS/QuST/Theme- 1/2019/23. SM acknowledges Ministry of Science and
Technology in Taiwan (Grant no. 110-2811-M-006 -501).

## Appendix: different unitary lead to different witness operator

It is to be noted that different unitary transformations make a PPT state to a
NPPT one in a different parameter-range. For example, motivated by the unitary
transformation in two-qubits, we write a global unitary operator in $2\otimes
4$ as

$U_{1}=\dfrac{1}{\sqrt{2}}\begin{pmatrix}1&0&0&0&0&0&0&1\\\
0&\sqrt{2}&0&0&0&0&0&0\\\ 0&0&\sqrt{2}&0&0&0&0&0\\\ 0&0&0&\sqrt{2}&0&0&0&0\\\
0&0&0&0&\sqrt{2}&0&0&0\\\ 0&0&0&0&0&\sqrt{2}&0&0\\\ 0&0&0&0&0&0&\sqrt{2}&0\\\
-1&0&0&0&0&0&0&1\end{pmatrix}$ (23)

while for $3\otimes 3$, a similar unitary transformation takes the form as

$U_{1}=\dfrac{1}{\sqrt{2}}\begin{pmatrix}1&0&0&0&0&0&0&0&1\\\
0&\sqrt{2}&0&0&0&0&0&0&0\\\ 0&0&\sqrt{2}&0&0&0&0&0&0\\\
0&0&0&\sqrt{2}&0&0&0&0&0\\\ 0&0&0&0&\sqrt{2}&0&0&0&0\\\
0&0&0&0&0&\sqrt{2}&0&0&0\\\ 0&0&0&0&0&0&\sqrt{2}&0&0\\\
0&0&0&0&0&0&0&\sqrt{2}&0\\\ -1&0&0&0&0&0&0&0&1\end{pmatrix}.$ (24)

Figs. 2 and 4, triangles symbolize the boundary between PPT states which can
be made NPPT via these unitary operators and the PPT states which remain PPT
even after their applications. Clearly, the above unitary operators are weaker
than $U$ in Eqs. (15) and (18).

## References

  * (1) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Rev. Mod. Phys. 81, 865 (2009).
  * (2) M. A. Nielsen and I. L. Chuang. Quantum Computation and Quantum Information. Cambridge University Press, 2000.
  * (3) A. K. Ekert, Phys. Rev. Lett. 67, 661 (1991).
  * (4) C. H. Bennett and S. 3. Wiesner, Phys. Rev. Lett. 69, 2881 (1992).
  * (5) C. H. Bennett, G. Brassard, C. Crepeau, R. Jozsa, A. Peres, and W. K. Wootters, Phys. Rev. Lett. 70, 1895 (1993).
  * (6) V. Giovannetti, S. Lloyd, and L. Maccone, 412, 417 (2001).
  * (7) A. K. Pati, Phys.Rev. A 63 (2000) 014320; H.-K. Lo, Phys. Rev. A 62, 012313 (2000).
  * (8) R. Raussendorf and H. J. Briegel, Phys. Rev. Lett. 86, 5188 (2001).
  * (9) O. Guhne and G. Toth, Phys. Rep. 474, 1 (2009).
  * (10) S. Das, T. Chanda, M. Lewenstein, A. Sanpera, A. Sen(De) and U. Sen, _The separability versus entanglement problem_ , in _Quantum Information: From Foundations to Quantum Technology Applications_ , second edition, eds. D. Bruß and G. Leuchs (Wiley, Weinheim, 2019), arXiv:1701.02187 [quant-ph].
  * (11) L. Gurvits, Proceedings of the thirty-fifth annual ACM symposium on Theory of computing, Eds. L. L. Larmore and M. X. Goemans, 10 (2003).
  * (12) S. Gharibian, Quant. Inf. and Comp. 10, 343 (2010).
  * (13) R. F. Werner, Phys. Rev. A 40, 4277 (1989).
  * (14) A. Peres, Phys. Rev. Lett. 77, 1413 (1996).
  * (15) M. Horodecki, P. Horodecki, and R. Horodecki, Phys. Lett. A 223, 1 (1996).
  * (16) B. M. Terhal, Phys. Lett. A 271, 319 (2000).
  * (17) O. Guhne, P. Hyllus, D. Bruß, A. Ekert, M. Lewenstein, C. Macchiavello and A. Sanpera, Phys. Rev. A 66, 062305 (2002).
  * (18) E. Knill, Separability from spectrum, published electronically at http://qig.itp.uni-hannover.de/qiproblems/15 (2003).
  * (19) M. Kus and K. Zyczkowski, Phys. Rev. A 63, 032307 (2001).
  * (20) F. Verstraete, K. Audenaert and B. D Moor, Phys. Rev. A 64, 012316 (2001).
  * (21) R. Hildebrand, Phys. Rev. A 76, 052325 (2007).
  * (22) P. B. Slater, J. Geom. Phys. 59, 17 (2009).
  * (23) N. Johnston, Phys. Rev. A 88, 062330 (2013).
  * (24) N. Ganguly, J. Chatterjee, and A. S. Majumdar, Phys. Rev. A 89, 052304 (2014).
  * (25) S. Arunachalam, N. Johnston, and V. Russo, Quant. Inf. Comput. 15, 0694 (2015).
  * (26) M. A. Jivulescu, N. Lupa, and I. Nechita, Quant. Inf. Comput. 15, 1165 (2015).
  * (27) O. Guhne, and N. Lutkenhaus, Phys. Rev. Lett. 96, 170502 (2006).
  * (28) O. Gühne, and N. Lutkenhaus, Journal of Physics: Conference Series 67, 012004 (2007).
  * (29) C.-J. Zhang, Y.-S. Zhang, S. Zhang, and G.-C. Guo, Phys. Rev. A 76, 012334 (2007).
  * (30) T. Moroder, O. Gühne, and N. Lütkenhaus, Phys. Rev. A 78, 032326 (2008).
  * (31) M. Kotowski, M. Kotowski, and M. Kus, Phys. Rev. A 81, 062318 (2010).
  * (32) J. M. Arrazola, O. Gittsovich, and N. Lütkenhaus, Phys. Rev. A 85 062327 (2012).
  * (33) J.-Y. Wu, H. Kampermann, D. Bruß, C. Klöckl, and M. Huber, Phys. Rev. A 86, 022319 (2012).
  * (34) K. Lemr, K. Bartkiewicz, and A. Cernoch, Phys. Rev. A 94, 052334 (2016).
  * (35) K. Sen, S. Das, and U. Sen, Phys. Rev. A 100, 062333 (2019).
  * (36) S.-Q. Shen, J.-M. Liang, M. Li, J. Yu, and S.-M. Fei, Phys. Rev. A 101, 012312 (2020).
  * (37) I. Bengtsson and K. Zyczkowski, Geometry of Quantum States (Cambridge University Press, Cambridge, 2006).
  * (38) M. A. Jivulescu, N. Lupa, I. Nechita, and D. Reeb, Linear Algebra Appl. 469, 276 (2015).
  * (39) S. Haldar, S. Mal, and A. Sen (De), (to be appeared in PRA).
  * (40) R. B. Holmes, Geometric Functional Analysis and its Applications, (Springer-Verlag,Berlin, 1975).
  * (41) M. Lewenstein, B. Kraus, J. I. Cirac, and P. Horodecki, Phys. Rev. A 62, 052310 (2000).
  * (42) P. Horodecki, Phys. Lett. A 232, 333 (1997).
  * (43) P. Horodecki, M. Horodecki, and R. Horodecki, Phys. Rev. Lett. 82, 1056 (1999).
  * (44) P. Skwara, H. Kampermann, M. Kleinmann, and D. Bruß, Phys. Rev. A 76, 012312 (2007).
  * (45) P. Pearle, Phys. Rev. D 2, 1418 (1970); N. Gisin and B. Gisin, Phys. Lett. A 260, 323 (1999); S. Massar, S. Pironio, J. Roland, and B. Gisin, Phys. Rev. A 66, 052112 (2002); A. Cabello, J.-A. Larsson, Phys. Rev. Lett. 98, 220402 (2007); T. Vertesi, S. Pironio, and N. Brunner, Phys. Rev. Lett. 104, 060401 (2010); C. Branciard, Phys. Rev. A 83, 032123 (2011); J. Szangolies, H. Kampermann, and D. Bruß, Phys. Rev. Lett. 118, 260401 (2017).
  * (46) K. Sen, C. Srivastava, S. Mal, A. Sen(De), and U. Sen, Phys. Rev. A 103, 032415 (2021).
  * (47) $d_{1}\otimes d_{2}$ represents the dimensions of the complex Hilbert space, $d_{1}d_{2}$ where $d_{1}$ and $d_{2}$ are the dimensions of subsystems $A$ and $B$ respectively.
  * (48) M. Jezek, J. Fiurasek, and Z. Hradil, Phys. Rev. A 68, 012305 (2003); A. Miranowicz, K. Bartkiewicz, J. Perina, Jr., M. Koashi, N. Imoto, and F. Nori, Phys. Rev. A 90, 062123 (2014); K. Bartkiewicz, A. Černoch, K. Lemr, and A. Miranowicz, Sci. Rep. 6, 19610 (2016).