collisions. In practice, this requires significant work that lies outside the
scope of the present manuscript. In section 5, we have found that a simple
reweighting of a vacuum-generated distribution of $g\to c\bar{c}$ by a factor
$\left(1+P_{g\to c\bar{c}}^{\text{med}}/P_{g\to
c\bar{c}}^{\text{vac}}\right)(z,k_{c}^{2},E_{g})$ can be used to estimate the
medium-modified $g\to c\bar{c}$ branching probability. In the absence of a
full medium-modified parton shower including $P_{g\to c\bar{c}}^{\text{med}}$,
this reweighting prescription will allow us to explore some experimental
signatures of enhanced $c\bar{c}$-radiation in a companion paper
Companion:2022 .

We finally note that the last decade has seen significant theoretical efforts
to improve over the formal accuracy of the original BDMPS-Z formulation. A
first understanding of medium-induced destructive interference in gluon
radiation from more than one parton was reached in studying the radiation
pattern of a QCD antenna in dense matter Armesto:2011ir ; Mehtar-Tani:2012mfa
. Radiative corrections were understood to lead to double logarithmic
enhancements of $\hat{q}$ with in-medium path-length that can be absorbed in a
renormalization of $\hat{q}$ Liou:2013qya ; Blaizot:2014bha ; Wu:2014nca . In
the future, these results may be understood as particular limiting cases of a
ground-breaking body of works by Arnold, Gorda and Iqbal that treats medium-
induced QCD bremsstrahlung beyond leading order in $\alpha_{s}$ and that
analyzes in particular the problem of overlapping formation times in
sequential bremsstrahlung gluons Arnold:2015qya ; Arnold:2020uzm . In parallel
to these theoretical developments, much work aimed at extending formulations
of medium-induced gluon radiation to multi-parton final states via the
development of dedicated Monte Carlo simulation tools Zapp:2008af ;
Zapp:2012ak ; Zapp:2013vla ; Armesto:2009fj ; Schenke:2009gb ; Caucal:2018ofz
; Putschke:2019yrg . Our main result, (6), is the first calculation of the
medium modified splitting to heavy and light quarks, thus completing the set
of BDMPS-Z splitting functions. This will contribute to the further refinement
of medium-modified parton showers. If the present calculation of the medium-
modified $g\to c\bar{c}$ splitting could be extended to $g\to c\bar{c}g$, this
would complement the recent developments of Armesto:2011ir ; Mehtar-
Tani:2012mfa ; Liou:2013qya ; Blaizot:2014bha ; Wu:2014nca ; Arnold:2015qya ;
Arnold:2020uzm .

###### Acknowledgements.

We thank Alexander Huss and Pier Francesco Monni for useful discussions on
parton showers. MA acknowledges support through H2020-MSCA-IF-2019 ExHolo
898223.

##  Appendix A Derivation of the splitting function (6) at $N=1$ opacity
expansion

In this appendix, we provide details of the derivation of (20). Expansion of
the integrand of (6) to first order in $n\,\sigma$ involves two contributions.
The first arises from expansion of the absorption factor
$\exp\left[-\tfrac{1}{2}\int_{\bar{y}}^{\infty}d\xi\,n(\xi)\,\sigma({\bf
r}_{\text{out}})\right]$, and it places a scattering centre at late times
$\xi>\bar{y}$ at which the $c\bar{c}$-pair is _fully formed (ff)_. The second
arises from expansion of the propagator ${\cal K}$, and it places a single
scattering centre at times $y<\xi<\bar{y}$ at which the $c\bar{c}$-pair is
_not fully formed (nff)_ (in the sense that it appears in the amplitude but
not yet in the complex conjugate). In the following, we evaluate these
contributions separately.

The _fully formed_ $N=1$ opacity contribution to (6) reads

$\displaystyle\left(\frac{1}{Q^{2}}\,P_{g\to
c\,\bar{c}}\right)^{\text{med}}_{\text{ff}}=-\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}\,\int_{0}^{\infty}dy\int_{y}^{\infty}d\bar{y}\,\int_{\bar{y}}^{\infty}d\xi\,n(\xi)\exp\left[i\frac{m_{\text{c}}^{2}}{2\mu}(y-\bar{y})-\epsilon|y|-\epsilon|\bar{y}|\right]\,$
$\displaystyle\int d{\bf r}_{\text{out}}\,\sigma({\bf
r}_{\text{out}})e^{-i\,{\bf k}_{\text{c}}.{\bf
r}_{\text{out}}}\left[\left(m_{\text{c}}^{2}+\frac{\partial}{\partial{\bf
r}_{\text{in}}}.\frac{\partial}{\partial{\bf
r}_{\text{out}}}\right)\frac{z^{2}+(1-z)^{2}}{z(1-z)}+2m_{c}^{2}\right]\,{\mathcal{K}}_{0}\left({\bf
r}_{\text{in}}=0,y;{\bf r}_{\text{out}},\bar{y}\right)$
$\displaystyle=\int\frac{d{\bf q}}{(2\pi)^{2}}\,|a({\bf
q})|^{2}\int_{0}^{\infty}d\xi
n(\xi)\Bigg{\\{}(1-\cos[\Gamma_{1}\xi])\left(\frac{1}{Q^{2}_{1}}\,P_{g\to
c\,\bar{c}}\right)^{\text{vac}}$ $\displaystyle\hskip
142.26378pt-(1-\cos[\Gamma_{0}\xi])\left(\frac{1}{Q^{2}_{0}}\,P_{g\to
c\,\bar{c}}\right)^{\text{vac}}\Bigg{\\}}\,.$ (104)

Here, to get to (104), one inserts (7) for the dipole cross section and one
does the ${\bf r}_{\text{out}}$ integration according to (11). Noting that

$\int_{0}^{\infty}dy\int_{y}^{\infty}d\bar{y}\,\int_{\bar{y}}^{\infty}d\xi=\int_{0}^{\infty}d\xi\int_{0}^{\xi}dy\int_{y}^{\xi}d\bar{y}\,,$
(105)

we can perform the longitudinal integrals for density profile $n(\xi)$ with
finite support

$\displaystyle\lim_{\epsilon\to
0}\frac{1}{4E_{g}^{2}}\int_{0}^{\xi}dy\int_{y}^{\xi}d\bar{y}\cos[\Gamma_{0}(\bar{y}-y)]e^{-\epsilon(y+\bar{y})}=\frac{1-\cos[\Gamma_{0}\xi]}{Q^{4}}\,.$
(106)

Analogously, one integrates over the phase factor set by the shifted
transverse energy $\Gamma_{1}$ in (21).

The _not fully formed_ $N=1$ opacity contribution to (6) is found with the
help of (18),

$\displaystyle\left(\frac{1}{Q^{2}}\,P_{g\to
c\,\bar{c}}\right)^{\text{med}}_{\text{nff}}=-\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}\,\int_{0}^{\infty}dy\int_{y}^{\infty}d\bar{y}\,\int^{\bar{y}}_{y}d\xi\,n(\xi)\exp\left[i\frac{m_{\text{c}}^{2}}{2\mu}(y-\bar{y})-\epsilon|y|-\epsilon|\bar{y}|\right]\,$
$\displaystyle\qquad\qquad\times\int d{\bf r}_{\text{out}}e^{-i\,{\bf
k}_{\text{c}}.{\bf
r}_{\text{out}}}\left[\left(m_{\text{c}}^{2}+\frac{\partial}{\partial{\bf
r}_{\text{in}}}.\frac{\partial}{\partial{\bf
r}_{\text{out}}}\right)\frac{z^{2}+(1-z)^{2}}{z(1-z)}+2m_{c}^{2}\right]\,$
$\displaystyle\qquad\qquad\times\int d{\bm{\uprho}}\,{\cal K}_{0}\left({\bf
r}_{\text{in}},y;{\bm{\uprho}},\xi|\mu\right)\,\sigma({\bm{\uprho}})\,{\cal
K}_{0}\left({\bm{\uprho}},\xi;{\bf r}_{\text{out}},\bar{y}|\mu\right)$
$\displaystyle=\int\frac{d{\bf q}}{(2\pi)^{2}}\,|a({\bf
q})|^{2}\int_{0}^{\infty}d\xi
n(\xi)\Bigg{\\{}(1-\cos[\Gamma_{0}\xi])\left(\frac{1}{Q^{2}_{0}}\,P_{g\to
c\,\bar{c}}\right)^{\text{vac}}$
$\displaystyle\qquad\qquad-(1-\cos[\Gamma_{1}\xi])\frac{1}{Q^{2}Q_{1}^{2}}\left[\left(m_{\text{c}}^{2}+{\bf
k}_{c}\cdot({\bf k}_{c}+{\bf
q}_{c})\right)\frac{z^{2}+(1-z)^{2}}{z(1-z)}+2m_{c}^{2}\right]\Bigg{\\}}\,.$
(107)

As the scattering centre is placed in this contribution at a position $\xi$
with $y<\xi<\bar{y}$, the derivative $\frac{\partial}{\partial{\bf
r}_{\text{out}}}$ that acts at time $\bar{y}$ picks up a transverse momentum
${\bf k}_{\text{c}}$, while the derivative $\frac{\partial}{\partial{\bf
r}_{\text{in}}}$ at earlier time $y$ picks up a shifted transverse momentum
${\bf k}_{\text{c}}+{\bf q}$. Compared to the fully formed contribution (104),
this complicates the transverse momentum dependence slightly.

The full $N=1$ opacity contribution to (6) is the sum of (104) and (107). In
this sum, contributions proportional to the interference factor
$(1-\cos[\Gamma_{0}\xi])$ cancel, and we are left with an integrand
proportional to $(1-\cos[\Gamma_{1}\xi])$. To reorganize this integrand, we
use

$\displaystyle\frac{-1}{Q^{2}Q_{1}^{2}}$
$\displaystyle=\frac{1}{2}\left(\frac{1}{Q_{1}^{2}}-\frac{1}{Q^{2}}\right)^{2}-\frac{1}{2}\frac{1}{Q_{1}^{4}}-\frac{1}{2}\frac{1}{Q^{4}}$
(108) $\displaystyle\frac{-{\bf k}_{c}\cdot\left({\bf k}_{c}+{\bf
q}\right)}{Q^{2}Q_{1}^{2}}$ $\displaystyle=\frac{1}{2}\left(\frac{{\bf
k}_{c}+{\bf q}}{Q_{1}^{2}}-\frac{{\bf
k}_{c}}{Q^{2}}\right)^{2}-\frac{1}{2}\frac{({\bf k}_{c}+{\bf
q})^{2}}{Q_{1}^{4}}-\frac{1}{2}\frac{{\bf k}_{c}^{2}}{Q^{4}}.$ (109)

With the help of (108) and (109), the sum of (104) and (107) takes the form

$\displaystyle\left(\frac{1}{Q^{2}}\,P_{g\to c\,\bar{c}}\right)^{\text{med}}=$
$\displaystyle\frac{1}{2}\int\frac{d{\bf q}}{(2\pi)^{2}}\,|a({\bf
q})|^{2}\int_{0}^{\infty}d\xi
n(\xi)\left(1-\cos\left[\Gamma_{1}\xi\right]\right)$
$\displaystyle\times\Bigg{[}$ $\displaystyle\left(\frac{1}{Q^{2}}\,P_{g\to
c\,\bar{c}}\right)^{\text{vac}}_{{\bf k}_{c}\to{\bf k}_{c}+{\bf
q}}-\left(\frac{1}{Q^{2}}\,P_{g\to c\,\bar{c}}\right)^{\text{vac}}$
$\displaystyle+\frac{m_{c}^{2}}{z(1-z)}\left(\frac{1}{Q_{1}^{2}}-\frac{1}{Q^{2}}\right)^{2}+\frac{z^{2}+(1-z)^{2}}{z(1-z)}\left(\frac{{\bf
k}_{c}+{\bf q}}{Q_{1}^{2}}-\frac{{\bf k}_{c}}{Q^{2}}\right)^{2}\Bigg{]}.$
(110)

This completes the derivation of (25). For static medium with
$n(\xi)=n_{0}\theta(\xi-L)\theta(\xi)$ the $\xi$ integral is trivial and we
obtain (20).

##  Appendix B Derivation of the splitting functions (54) and (55) in multiple
soft scattering approximation for static medium

To derive (54) and (55), we focus on two contributions of the full splitting
function (6) defined in terms of the integrals in (51): $I_{4}$ and $I_{5}$.
In the saddle-point approximation (45) and for constant $\hat{q}(\xi)=\hat{q}$
in the region $0<\xi<L$ the term $I_{4}$ receives $\hat{q}$ contributions from
both the absorption factor and the path integral of the two dimensional
harmonic oscillator (47). For $I_{5}$ the absorption factor is unity, but the
propagator is a convolution of the harmonic oscillator with the free
propagator. Our goal is to perform the remaining transverse and one
longitudinal integral.

###  B.1 $I_{4}$

We first substitute the propagator expression for the harmonic oscillator (47)
in (6). Then we take the $\partial_{{\bf r}_{\text{in}}}$ derivative,
integrate ${\bf r}_{\text{out}}$ by parts and set ${\bf r}_{\text{in}}=0$. We
obtain

$\displaystyle I_{4}$
$\displaystyle=2\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}\,\int_{0}^{L}dy\int_{y}^{L}d\bar{y}\,\exp\left[i\frac{m_{\text{c}}^{2}}{2\mu}(y-\bar{y})\right]\,\int
d{\bf r}_{\text{out}}e^{-\frac{1}{4}\hat{q}(L-\bar{y}){\bf
r}_{\text{out}}^{2}-i\,{\bf k}_{\text{c}}.{\bf r}_{\text{out}}}$
$\displaystyle\times\left[\mu\Omega\frac{{\bf r}_{\text{out}}{\bf
k}_{c}-i\frac{1}{2}\hat{q}(L-\bar{y}){\bf
r}_{\text{out}}^{2}}{\sin\Omega(\bar{y}-y)}\frac{z^{2}+(1-z)^{2}}{z(1-z)}+\frac{m_{c}^{2}}{z(1-z)}\right]\frac{\mu\Omega}{2\pi
i\sin\Omega(\bar{y}-y)}e^{-\frac{-i2\mu\Omega\cot\Omega(\bar{y}-y)}{4}{\bf
r}_{\text{out}}^{2}}$ (111)

Noticing that the trigonometric functions depend only on the difference of
longitudinal coordinates, we make the change of variables to $x=\bar{y}-y$ and
$u=L-\bar{y}$. The integration bounds are changed using the following
equivalence of inequalities

$0<y<L\,,\quad y<\bar{y}<L\,,\quad\leftrightarrow\quad 0<x<L\,,\quad
0<u<L-x\,.$ (112)

Expressing the transverse integral in polar coordinates, we write

$\displaystyle I_{4}$
$\displaystyle=2\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}\,\int_{0}^{L}dx\int_{0}^{L-x}du\,\exp\left[-i\frac{m_{\text{c}}^{2}}{2\mu}x\right]\,\int\frac{1}{2}d{r}^{2}_{\text{out}}d\phi
e^{-\frac{1}{4}\hat{q}u{r}_{\text{out}}^{2}-i\,{k}_{\text{c}}{r}_{\text{out}}\cos\phi}$
$\displaystyle\times\left[\mu\Omega\frac{{r}_{\text{out}}{k}_{c}\cos\phi-i\frac{1}{2}\hat{q}u{r}_{\text{out}}^{2}}{\sin\Omega
x}\frac{z^{2}+(1-z)^{2}}{z(1-z)}+\frac{m_{c}^{2}}{z(1-z)}\right]\frac{\mu\Omega}{2\pi
i\sin\Omega x}e^{-\frac{-i2\mu\Omega\cot\Omega x}{4}{r}_{\text{out}}^{2}}\,.$
(113)

Here, the $u$-integral is done easily and the azimuthal integral over $\phi$
leads to Bessel functions. Using
$r_{\text{out}}J_{1}(k_{c}r_{\text{out}})=-\frac{\partial}{\partial
k_{c}}J_{0}(k_{c}r_{\text{out}})$, we find

$\displaystyle I_{4}$
$\displaystyle=2\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}\,\int_{0}^{L}dx\,\exp\left[-i\frac{m_{\text{c}}^{2}}{2\mu}x\right]\,\int_{0}^{\infty}{r}_{\text{out}}d{r}_{\text{out}}e^{-\frac{1}{8}\hat{q}(L-x){r}_{\text{out}}^{2}}$
$\displaystyle\times\Bigg{\\{}\frac{\sinh[\frac{1}{8}\hat{q}(L-x){r}_{\text{out}}^{2}]}{\frac{1}{8}\hat{q}{r}_{\text{out}}^{2}}\left[\mu\Omega\frac{i{k}_{c}\partial_{k_{c}}-2i}{\sin\Omega
x}\frac{z^{2}+(1-z)^{2}}{z(1-z)}+\frac{m_{c}^{2}}{z(1-z)}\right]$
$\displaystyle+i2\mu\Omega\frac{(L-x)e^{-\frac{1}{8}\hat{q}(L-x){r}_{\text{out}}^{2}}}{\sin\Omega
x}\frac{z^{2}+(1-z)^{2}}{z(1-z)}\Bigg{\\}}\frac{\mu\Omega
J_{0}(k_{c}r_{\text{out}})}{i\sin\Omega x}e^{-\frac{-i2\mu\Omega\cot\Omega
x}{4}{r}_{\text{out}}^{2}}\,.$ (114)

The remaining radial integral can be done, using the identities

$\displaystyle\int_{0}^{\infty}dte^{-\beta
t}J_{0}(k_{c}\sqrt{t})=\frac{1}{\beta}e^{-\frac{k_{c}^{2}}{4\beta}}\,,$ (115)
$\displaystyle\int_{0}^{\infty}dte^{-\beta t}\frac{\sinh[\gamma
t]}{t}J_{0}(k_{c}\sqrt{t})=\frac{1}{2}\text{Ei}\left[-\frac{k^{2}_{c}}{4\left(\beta-\gamma\right)}\right]-\frac{1}{2}\text{Ei}\left[-\frac{k^{2}_{c}}{4\left(\beta+\gamma\right)}\right]\,,$
(116)

where $4(\beta-\gamma)=-2i\mu\Omega\cot\Omega x$ and
$4(\beta+\gamma)=\hat{q}(L-x)-i2\mu\Omega\cot\Omega x$. Here Ei is the
exponential integral function. We obtain

$\displaystyle I_{4}$
$\displaystyle=\frac{2\mu}{\hat{q}z(1-z)}2\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}\,\int_{0}^{L}dx\frac{-i\Omega}{\sin\Omega
x}e^{-i\frac{m_{\text{c}}^{2}}{2\mu}x}\Bigg{\\{}\frac{-i2\mu\Omega}{\sin\Omega
x}(z^{2}+(1-z)^{2})\Bigg{(}\text{Ei}\left[\frac{k^{2}_{c}\tan\Omega
x}{i2\mu\Omega}\right]$
$\displaystyle-\text{Ei}\left[\frac{k^{2}_{c}\tan\Omega
x}{i2\mu\Omega-\hat{q}(L-x)\tan\Omega
x}\right]-\exp\left[\frac{k^{2}_{c}\tan\Omega
x}{i2\mu\Omega}\right]+\frac{\exp\left[-\frac{k^{2}_{c}\tan\Omega
x}{i2\mu\Omega-\hat{q}(L-x)\tan\Omega
x}\right]}{1-\frac{\hat{q}(L-x)}{i2\mu\Omega}\tan\Omega x}\Bigg{)}$
$\displaystyle+m_{c}^{2}\left(\text{Ei}\left[\frac{k^{2}_{c}\tan\Omega
x}{i2\mu\Omega}\right]-\text{Ei}\left[\frac{k^{2}_{c}\tan\Omega
x}{i2\mu\Omega-\hat{q}(L-x)\tan\Omega x}\right]\right)\Bigg{\\}}\,.$ (117)

Expressing this equation in the dimensionless variables (56), we obtain (54).

###  B.2 $I_{5}$

For the case $0<y<L$ and $L<\bar{y}$ with constant $\hat{q}(\xi)=\hat{q}$ in
the region $0<\xi<L$, the solution for the path integral is given by the
convolution of the harmonic oscillator propagator (47) and the free propagator
(10)

$\displaystyle{\cal K}\left({\bf r}_{\text{in}}=0,y;{\bf
r}_{\text{out}},\bar{y}|\mu\right)=$ $\displaystyle\int d{\bf x}\,{\cal
K}_{\text{osc}}\left({\bf r}_{\text{in}}=0,y;{\bf x},L|\mu\right){\cal
K}_{0}\left({\bf x},L;{\bf r}_{\text{out}},\bar{y}|\mu\right).$ (118)

We substitute this expression in (6) and we integrate ${\bf r}_{\text{out}}$
by parts. We then integrate over ${\bf r}_{\text{out}}$ using (11), we take
the $\partial_{{\bf r}_{\text{in}}}$ derivative and we set ${\bf
r}_{\text{in}}=0$. In this way, we obtain

$\displaystyle I_{5}$
$\displaystyle=2\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}\,\int_{0}^{L}dy\int_{L}^{y_{\infty}}d\bar{y}\,\exp\left[i\frac{m_{\text{c}}^{2}}{2\mu}(y-\bar{y})-\epsilon\bar{y}\right]\,\int
d{\bf x}e^{-i(\bar{y}-L)\frac{{\bf k}_{c}^{2}}{2\mu}-i\,{\bf
k}_{\text{c}}.{\bf x}}$
$\displaystyle\times\frac{\mu}{z(1-z)}\left[m_{\text{c}}^{2}+\mu\Omega\frac{{\bf
x}\cdot{\bf
k}_{c}}{\sin\Omega(L-y)}(z^{2}+(1-z)^{2})\right]\frac{-i\Omega}{2\pi\sin\Omega(L-y)}e^{-\frac{-i2\mu\Omega{\bf
x}^{2}\cot\Omega(L-y)}{4}}\,.$ (119)

Doing the remaining Gaussian transverse integral, we find

$\displaystyle I_{5}$
$\displaystyle=2\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}\,\int_{0}^{L}dy\int_{L}^{y_{\infty}}d\bar{y}\,\exp\left[i\frac{m_{\text{c}}^{2}}{2\mu}(y-\bar{y})-\epsilon\bar{y}\right]\,$
$\displaystyle\times\frac{1}{\cos\Omega(L-y)}e^{-i(\bar{y}-L)\frac{{\bf
k}_{c}^{2}}{2\mu}}\frac{1}{z(1-z)}e^{-\frac{{\bf
k}_{c}^{2}}{-2i\mu\Omega\cot\Omega x}}\left[m_{\text{c}}^{2}+\frac{{\bf
k}_{c}^{2}}{\cos\Omega(L-y)}(z^{2}+(1-z)^{2})\right]\,.$ (120)

We again notice that the trigonometric functions only depend on the difference
$L-y$. Therefore we make the change of variables $x=L-y$ and $u=\bar{y}-L$.
The integration bounds are determined by the equivalence of

$0<y<L,\quad L<\bar{y}<\infty\quad\leftrightarrow\quad 0<x<L,\quad
0<u<\infty.$ (121)

The $u$ integral is simple and yields

$\displaystyle I_{5}$
$\displaystyle=\frac{2\mu}{z(1-z)}2\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}\,\int_{0}^{L}dx\,\frac{-i}{\cos\Omega
x}\exp\left[-i\frac{m_{\text{c}}^{2}}{2\mu}x+\frac{{\bf k}_{c}^{2}\tan\Omega
x}{2i\mu\Omega}\right]$ $\displaystyle\qquad\times\frac{1}{m_{c}^{2}+{\bf
k}_{c}^{2}}\left[m_{\text{c}}^{2}+\frac{{\bf k}_{c}^{2}}{\cos\Omega
x}(z^{2}+(1-z)^{2})\right]\,.$ (122)

Expressing this equation in the dimensionless variables (56), we obtain (55).

##  Appendix C Derivation of the splitting functions (59) and (60) in the
multiple soft scattering approximation for an expanding medium

For a general time dependent quenching parameter $\hat{q}(\xi)$ with finite
support for $0<\xi<L$, the two-dimensional path integral in (8) is given by
Baier:1998yf ; schulman2012techniques 888See also the lecture notes on the
single-dimensional harmonic oscillator by Andreas Wipf at https://www.tpi.uni-
jena.de/~wipf/lectures/pfad/pfad3

${\cal K}\left({\bf r}_{\text{in}}=0,y;{\bf
r}_{\text{out}},\bar{y}|\mu\right)=\frac{\mu}{2\pi
iD(\bar{y},y)}e^{i\frac{\mu}{2D(\bar{y},y)}\left({\bf
r}^{2}_{\text{out}}\frac{\partial}{\partial\bar{y}}D(\bar{y},y)-{\bf
r}^{2}_{\text{in}}\frac{\partial}{\partial{y}}D(\bar{y},y)-2{\bf
r}_{\text{out}}\cdot{\bf r}_{\text{in}}\right)},$ (123)

where the auxiliary function $D(\bar{y},y)$ solves the second order
differential equation

$\frac{\partial^{2}}{\partial^{2}\xi}D(\xi,y)=-\Omega^{2}(\xi)D(\xi,y),\quad
D(y,y)=0,\quad\left.\frac{\partial D(\xi,y)}{\partial\xi}\right|_{\xi=y}=1,$
(124)

with $\Omega^{2}(\xi)=\frac{\hat{q}(\xi)}{2i\mu}$. For constant $\hat{q}$, the
solution is $D(\bar{y},y)=\frac{1}{\Omega}\sin\Omega(\bar{y}-y)$ and we
recover (47). In the following, we repeat the steps of Appendix B for this
more general time-dependent case.

###  C.1 $I_{4}$

We first substitute the propagator expression for generalized harmonic
oscillator (123) in (6). Then we take the $\partial_{{\bf r}_{\text{in}}}$
derivative, integrate ${\bf r}_{\text{out}}$ by parts and set ${\bf
r}_{\text{in}}=0$. We obtain

$\displaystyle I_{4}$
$\displaystyle=2\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}\,\int_{0}^{L}dy\int_{y}^{L}d\bar{y}\,\exp\left[-i\frac{m_{\text{c}}^{2}}{2\mu}(\bar{y}-y)\right]\,\int
d{\bf r}_{\text{out}}e^{-\frac{1}{4}\int_{\bar{y}}^{L}d\xi\,\hat{q}(\xi){\bf
r}^{2}_{\text{out}}}e^{-i\,{\bf k}_{\text{c}}\cdot{\bf r}_{\text{out}}}$
$\displaystyle\qquad\qquad\times\left[\frac{\mu}{D(\bar{y},y)}\left({\bf
r}_{\text{out}}\cdot{\bf
k}_{c}-i\frac{1}{2}\int_{\bar{y}}^{L}d\xi\,\hat{q}(\xi){\bf
r}^{2}_{\text{out}}\right)\frac{z^{2}+(1-z)^{2}}{z(1-z)}+\frac{m_{c}^{2}}{z(1-z)}\right]$
$\displaystyle\qquad\qquad\times\frac{\mu}{2\pi
iD(\bar{y},y)}\exp\left[-\frac{1}{4}\\{-2i\mu\partial_{\bar{y}}\log
D(\bar{y},y)\\}{\bf r}^{2}_{\text{out}}\right]\,.$ (125)

Performing the remaining Gaussian integration yields

$\displaystyle I_{4}$
$\displaystyle=2\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}\,\int_{0}^{L}dy\int_{y}^{L}d\bar{y}\,\exp\left[-i\frac{m_{\text{c}}^{2}}{2\mu}(\bar{y}-y)\right]\,\Bigg{[}\frac{m_{c}^{2}}{z(1-z)}+\frac{z^{2}+(1-z)^{2}}{z(1-z)}$
$\displaystyle\times\left(\frac{\partial_{\bar{y}}D(\bar{y},y){\bf
k}^{2}_{c}}{\\{\partial_{\bar{y}}D(\bar{y},y)-D(\bar{y},y)\int_{\bar{y}}^{L}d\xi\,\Omega^{2}(\xi)\\}^{2}}+\frac{2i\mu\int_{\bar{y}}^{L}d\xi\,\Omega^{2}(\xi)}{\partial_{\bar{y}}D(\bar{y},y)-D(\bar{y},y)\int_{\bar{y}}^{L}d\xi\,\Omega^{2}(\xi)}\right)\Bigg{]}$
$\displaystyle\times\frac{1}{\partial_{\bar{y}}D(\bar{y},y)-D(\bar{y},y)\int_{\bar{y}}^{L}d\xi\,\Omega^{2}(\xi)}\exp\left[-i\frac{{\bf
k}_{c}^{2}}{2\mu}\\{\partial_{\bar{y}}\log
D(\bar{y},y)-\int_{\bar{y}}^{L}d\xi\,\Omega^{2}(\xi)\\}^{-1}\right]\,.$ (126)

Switching to dimensionless variables in (58) we obtain (59).

###  C.2 $I_{5}$

We substitute generalized harmonic oscillator path integral (123) in (6)
(using (118)) and integrate ${\bf r}_{\text{out}}$ by parts. Then we can
perform ${\bf r}_{\text{out}}$ integral using (11). Finally we take
$\partial_{{\bf r}_{\text{in}}}$ derivative and set ${\bf r}_{\text{in}}=0$.
We obtain

$\displaystyle I_{5}=2\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}$
$\displaystyle\,\int_{0}^{L}dy\int_{L}^{y_{\infty}}d\bar{y}\,\exp\left[-i\frac{m_{\text{c}}^{2}}{2\mu}(\bar{y}-y)-\epsilon\bar{y}\right]\,\int
d{\bf x}$ $\displaystyle\times e^{-i\,{\bf k}_{\text{c}}.{\bf
x}-i(\bar{y}-L)\frac{{\bf
k}_{c}^{2}}{2\mu}}\left[\left(m_{\text{c}}^{2}+\frac{\mu{\bf x}\cdot{\bf
k}_{c}}{D(L,y)}\right)\frac{z^{2}+(1-z)^{2}}{z(1-z)}+2m_{c}^{2}\right]$
$\displaystyle\times\frac{-i\mu}{2\pi
D(\bar{y},y)}\exp\left[-\frac{1}{4}(-2i\mu\left.\partial_{\bar{y}}\log
D(\bar{y},y)\right|_{\bar{y}=L}){\bf r}^{2}_{\text{out}}\right]\,.$ (127)

We then perform the radial integral

$\displaystyle I_{5}$
$\displaystyle=2\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}\,\int_{0}^{L}dy\int_{L}^{y_{\infty}}d\bar{y}\,\exp\left[-i\frac{m_{\text{c}}^{2}}{2\mu}(\bar{y}-y)-\epsilon\bar{y}-i\frac{{\bf
k}_{c}^{2}}{2\mu}(\\{\left.\partial_{\bar{y}}\log
D(\bar{y},y)\right|_{\bar{y}=L}\\}^{-1}-L+\bar{y})\right]$
$\displaystyle\qquad\qquad\times\left[\frac{{\bf
k}_{c}^{2}}{\left.\partial_{\bar{y}}D(\bar{y},y)\right|_{\bar{y}=L}}\frac{z^{2}+(1-z)^{2}}{z(1-z)}+\frac{m_{c}^{2}}{z(1-z)}\right]\frac{1}{\left.\partial_{\bar{y}}D(\bar{y},y)\right|_{\bar{y}=L}}$
(128)

In this term, also the $\bar{y}$ integral can be done analytically,

$\displaystyle I_{5}$
$\displaystyle=2\,\mathfrak{Re}\,\frac{1}{4\,E_{g}^{2}}\,\int_{0}^{L}dy\frac{-i2\mu}{m_{c}^{2}+{\bf
k}_{c}^{2}}\,\exp\left[-i\frac{m_{\text{c}}^{2}}{2\mu}(L-y)-\epsilon\bar{y}-i\frac{{\bf
k}_{c}^{2}}{2\mu}\\{\left.\partial_{\bar{y}}\log
D(\bar{y},y)\right|_{\bar{y}=L}\\}^{-1}\right]$
$\displaystyle\qquad\qquad\times\left[\frac{{\bf
k}_{c}^{2}}{\left.\partial_{\bar{y}}D(\bar{y},y)\right|_{\bar{y}=L}}\frac{z^{2}+(1-z)^{2}}{z(1-z)}+\frac{m_{c}^{2}}{z(1-z)}\right]\frac{1}{\left.\partial_{\bar{y}}D(\bar{y},y)\right|_{\bar{y}=L}}\,.$
(129)

Switching to dimensionless variables in (58) we obtain (60).

##  Appendix D Comments on the relation between $N=1$ hard and multiple soft
scattering

Our numerical study of the $N=1$ opacity expansion in section 3.5 was based on
a medium parameterized by Yukawa-type scattering centers (34). We refer to it
as $N=1$ _hard_ since these scattering centers display power-law tails that
one expects to arise in scatterings with sufficiently high (_hard_) momentum
transfer. In a multiple scattering picture, this accounts for Molière-type
Moliere:1948zz ; ParticleDataGroup:2020ssz  rare large-angle scattering. For
such a Yukawa potential, the dipole cross section (7) takes the form

$\sigma({\bf
r}^{2})\propto\left(1-\mu_{D}rK_{1}(\mu_{D}r)\right)\xrightarrow[]{r\mu_{D}\ll
1}\frac{\mu_{D}^{2}}{4}\left(1-2\gamma_{E}+\log\left[\frac{4}{\mu_{D}r}\right]\right){\bf
r}^{2}\,,$ (130)

where $K_{1}$ denotes a modified Bessel function and $\gamma_{E}=0.577....$ is
Euler’s constant. Here, the logarithmic correction
$\log\left[\frac{4}{\mu_{D}r}\right]$ prevents us from relating the prefactor
of ${\bf r}^{2}$ to the quenching parameter used in the saddle point
approximation (45). The problem is that (45) assumes the existence of the
second moment $\int\tfrac{d{\bf q}^{2}}{(2\pi)^{2}}|a({\bf q})|^{2}\,{\bf
q}^{2}$ whereas this moment is logarithmically UV divergent for a Yukawa-type
potential. The term $\log\left[\tfrac{4}{\mu_{D}r}\right]$ is the tell-tale
sign of this logarithmic divergence in configuration space. The saddle point
approximation (45) amounts to cutting off the ${\bf q}^{-4}$ tail of $|a({\bf
q})|^{2}$ at large ${\bf q}$ by regulating the $\log r$-dependence at small
$r$.

In contrast, by construction, a Gaussian ansatz for $|a({\bf q})|^{2}$ is free
of UV problems and it is consistent with the dipole approximation (45) since

$|a({\bf q})|^{2}=\frac{16\pi}{\langle{\bf q}^{2}\rangle_{\text{med}}}e^{-{\bf
q}^{2}/\langle{\bf q}^{2}\rangle_{\text{med}}}\Longrightarrow\sigma({\bf
r}^{2})=4\left(1-e^{-\langle{\bf q}^{2}\rangle_{\text{med}}{\bf
r}^{2}/4}\right)=\langle{\bf q}^{2}\rangle_{\text{med}}{\bf
r}^{2}+\mathcal{O}({\bf r}^{4})\,.$ (131)

We consider this Gaussian ansatz with a normalization such that the opacity
(41) is unity,

$\frac{1}{2}n_{0}\,L\int\frac{d{\bf q}}{(2\pi)}|a({\bf
q})|^{2}\dots\equiv\frac{1}{\langle{\bf
q}^{2}\rangle_{\text{med}}}\int_{0}^{\infty}d(q^{2})\,\exp\left[-\frac{q^{2}}{\langle{\bf
q}^{2}\rangle_{\text{med}}}\right]\dots\quad\,.$ (132)

We now ask to what extent the resulting $\left(\tfrac{1}{Q^{2}}P_{g\to
c\bar{c}}\right)^{\text{med}}$ differs from the medium-modification calculated
for Yukawa-type potentials with hard power-low tails.

(a)

(b)

Figure 15: The medium-modification (37) of the $g\to c\bar{c}$ splitting in
the $N=1$ opacity expansion, plotted as a function of $z$ and final transverse
momentum $k_{\text{c}}^{2}$. Upper and lower panels show results for different
values of the average squared momentum transfer from the medium. Values for
$\left(\tfrac{1}{Q^{2}}\,P_{g\to c\,\bar{c}}\right)^{\rm med}$ are displayed
in units of $\text{GeV}^{-2}$.

Using (132) amounts to replacing in the expression (37) for
$\left(\tfrac{1}{Q^{2}}P_{g\to c\bar{c}}\right)^{\text{med}}$

$c_{\text{opaque}}\frac{1}{2\mu_{D}^{2}}\frac{\tfrac{1}{2}}{\left(\tfrac{1}{2}+\tilde{\bf
q}^{2}\right)^{2}}\quad\longrightarrow\quad\frac{1}{\langle{\bf
q}^{2}\rangle_{\text{med}}}e^{-\tilde{\bf q}^{2}}\,.$ (133)

The dimensionless variables in (37) are then understood as measuring masses
and momenta in units of ${\langle{\bf q}^{2}\rangle_{\text{med}}}$,

$\tilde{m}_{\text{c}}^{2}=\frac{{m}_{\text{c}}^{2}}{\langle{\bf
q}^{2}\rangle_{\text{med}}}\,,\qquad\tilde{\bf p}_{\text{c}}^{2}=\frac{{\bf
p}_{\text{c}}^{2}}{\langle{\bf q}^{2}\rangle_{\text{med}}},\qquad\tilde{\bf
q}^{2}=\frac{{\bf q}^{2}}{\langle{\bf q}^{2}\rangle_{\text{med}}}\,,\qquad
E_{g}\equiv e_{g}\frac{1}{2}{\langle{\bf q}^{2}\rangle_{\text{med}}}L\,.$
(134)

Numerical results for this $N=1$ opacity expansion with Gaussian momentum
distribution (132) are displayed in Fig. 15. Close comparison with the results
for $N=1$ hard in Fig. 3 shows that the hard Yukawa-type tails of single
scattering centers yield only to a very mild if not negligible change in the
$k_{\text{c}}^{2}$-distribution. This can be understood on general grounds:
2/3rd (8/9th) of the Yukawa cross section in (133) lies in the region of soft
momentum transfers $\tilde{q}<1$ ($\tilde{q}<2$). Therefore, also the $N=1$
opacity approximation is dominated by small-angle scattering and differences
to a Gaussian distribution may be expected to be only gradual.

There has been repeated interest in isolating Molière-type large angle
scattering contributions to jet quenching as these may inform us about the
nature of the constituents in the QCD plasma. For a recent summary of the
phenomenological motivation, see Ref. DEramo:2018eoy . That Molière-type
large-angle scatterings can be included in jet quenching formulations via a
logarithmic dependence of the dipole cross section has been known since
Baier:1996kr ; Zakharov:1996fv . Indeed, in describing $e^{+}e^{-}$ QED pair
radiation off extended targets, Yukawa-type potentials of the form (34) had
been used to describe the screening of nuclei in QED matter by their electron
clouds; the observation of a logarithmic term in the corresponding dipole
cross section had been made in this context much earlier Zakharov:1987QED ;
Zakharov:1996fv . Conceptually the same logarithmic dependence appears also in
the configuration-space description of QCD saturation physics via Glauber-
Mueller multiple scattering effects Kovchegov:1996ty ; Jalilian-Marian:1996mkd
. In applying this formalism to the longitudinal dynamics of incoming partons
that undergo multiple scattering in a nuclear target, it was shown that the
rare Molière-type large-angle scattering encoded in an $\log r$-term could
result in a Cronin-type enhancement at a characteristic intermediate
transverse momentum scale Baier:2003hr . In quenched jet fragmentation
patterns, however, parametric reasoning indicates that Molière-type large-
angle contributions can become numerically relevant only in a very narrow
region of the entire phase space Kurkela:2014tla . We are not aware of any
proposal that would allow for a clear identification of such a large-angle
scattering contribution within realistic hadronic distributions. Consistent
with these findings, the comparison between $N=1$ hard and $N=1$ Gaussian
presented here indicates that also the medium-modification of $g\to c\bar{c}$
shows only a very small if not negligible dependence on the presence of a
large-angle scattering contribution.

## References

  * (1) M. Connors, C. Nattrass, R. Reed and S. Salur, _Jet measurements in heavy ion physics_ , _Rev. Mod. Phys._ 90 (2018) 025005 [1705.01974].
  * (2) J.D. Bjorken, _Energy Loss of Energetic Partons in Quark - Gluon Plasma: Possible Extinction of High p(t) Jets in Hadron - Hadron Collisions_ , .
  * (3) X.-N. Wang, M. Gyulassy and M. Plumer, _The LPM effect in QCD and radiative energy loss in a quark gluon plasma_ , _Phys. Rev. D_ 51 (1995) 3436 [hep-ph/9408344].
  * (4) M. Gyulassy and X.-n. Wang, _Multiple collisions and induced gluon Bremsstrahlung in QCD_ , _Nucl. Phys. B_ 420 (1994) 583 [nucl-th/9306003].
  * (5) R. Baier, Y.L. Dokshitzer, A.H. Mueller, S. Peigne and D. Schiff, _Radiative energy loss of high-energy quarks and gluons in a finite volume quark - gluon plasma_ , _Nucl. Phys. B_ 483 (1997) 291 [hep-ph/9607355].
  * (6) R. Baier, Y.L. Dokshitzer, A.H. Mueller, S. Peigne and D. Schiff, _Radiative energy loss and p(T) broadening of high-energy partons in nuclei_ , _Nucl. Phys. B_ 484 (1997) 265 [hep-ph/9608322].
  * (7) B.G. Zakharov, _Fully quantum treatment of the Landau-Pomeranchuk-Migdal effect in QED and QCD_ , _JETP Lett._ 63 (1996) 952 [hep-ph/9607440].
  * (8) B.G. Zakharov, _Radiative energy loss of high-energy quarks in finite size nuclear matter and quark - gluon plasma_ , _JETP Lett._ 65 (1997) 615 [hep-ph/9704255].
  * (9) S. Cao and X.-N. Wang, _Jet quenching and medium response in high-energy heavy-ion collisions: a review_ , _Rept. Prog. Phys._ 84 (2021) 024301 [2002.04028].
  * (10) G.-Y. Qin and X.-N. Wang, _Jet quenching in high-energy heavy-ion collisions_ , _Int. J. Mod. Phys. E_ 24 (2015) 1530014 [1511.00790].
  * (11) R. Baier, Y.L. Dokshitzer, A.H. Mueller and D. Schiff, _Radiative energy loss of high-energy partons traversing an expanding QCD plasma_ , _Phys. Rev. C_ 58 (1998) 1706 [hep-ph/9803473].
  * (12) C.A. Salgado and U.A. Wiedemann, _Calculating quenching weights_ , _Phys. Rev. D_ 68 (2003) 014008 [hep-ph/0302184].
  * (13) U.A. Wiedemann, _Gluon radiation off hard quarks in a nuclear environment: Opacity expansion_ , _Nucl. Phys. B_ 588 (2000) 303 [hep-ph/0005129].
  * (14) U.A. Wiedemann, _Jet quenching versus jet enhancement: A Quantitative study of the BDMPS-Z gluon radiation spectrum_ , _Nucl. Phys. A_ 690 (2001) 731 [hep-ph/0008241].
  * (15) R. Baier, Y.L. Dokshitzer, A.H. Mueller and D. Schiff, _On the angular dependence of the radiative gluon spectrum_ , _Phys. Rev. C_ 64 (2001) 057902 [hep-ph/0105062].
  * (16) B.G. Zakharov, _Transverse spectra of radiation processes in-medium_ , _JETP Lett._ 70 (1999) 176 [hep-ph/9906536].
  * (17) Y.L. Dokshitzer and D.E. Kharzeev, _Heavy quark colorimetry of QCD matter_ , _Phys. Lett. B_ 519 (2001) 199 [hep-ph/0106202].
  * (18) M. Gyulassy, P. Levai and I. Vitev, _Reaction operator approach to nonAbelian energy loss_ , _Nucl. Phys. B_ 594 (2001) 371 [nucl-th/0006010].
  * (19) X.-N. Wang and X.-f. Guo, _Multiple parton scattering in nuclei: Parton energy loss_ , _Nucl. Phys. A_ 696 (2001) 788 [hep-ph/0102230].
  * (20) M. Cacciari, S. Frixione, N. Houdeau, M.L. Mangano, P. Nason and G. Ridolfi, _Theoretical predictions for charm and bottom production at the LHC_ , _JHEP_ 10 (2012) 137 [1205.6344].
  * (21) M.L. Mangano, P. Nason and G. Ridolfi, _Heavy quark correlations in hadron collisions at next-to-leading order_ , _Nucl. Phys. B_ 373 (1992) 295.
  * (22) R.K. Ellis and J.C. Sexton, _Explicit Formulae for Heavy Flavor Production_ , _Nucl. Phys. B_ 282 (1987) 642.
  * (23) A. Andronic et al., _Heavy-flavour and quarkonium production in the LHC era: from proton–proton to heavy-ion collisions_ , _Eur. Phys. J. C_ 76 (2016) 107 [1506.03981].
  * (24) M. Cacciari, G.P. Salam and G. Soyez, _FastJet User Manual_ , _Eur. Phys. J. C_ 72 (2012) 1896 [1111.6097].
  * (25) A.J. Larkoski, S. Marzani, G. Soyez and J. Thaler, _Soft Drop_ , _JHEP_ 05 (2014) 146 [1402.2657].
  * (26) P. Ilten, N.L. Rodd, J. Thaler and M. Williams, _Disentangling Heavy Flavor at Colliders_ , _Phys. Rev. D_ 96 (2017) 054019 [1702.02947].
  * (27) C. Aidala et al., _sPHENIX: An Upgrade Concept from the PHENIX Collaboration_ , 1207.6378.
  * (28) Z. Citron et al., _Report from Working Group 5: Future physics opportunities for high-density QCD at the LHC with heavy-ion and proton beams_ , _CERN Yellow Rep. Monogr._ 7 (2019) 1159 [1812.06772].
  * (29) D. Adamová et al., _A next-generation LHC heavy-ion experiment_ , 1902.01211.
  * (30) M. Attems, J. Brewer, G.M. Innocenti, A. Mazeliauskas, S. Park, W. van der Schee et al., _Medium-enhanced c-cbar production_ , .
  * (31) N. Armesto et al., _Comparison of Jet Quenching Formalisms for a Quark-Gluon Plasma ’Brick’_ , _Phys. Rev. C_ 86 (2012) 064904 [1106.1106].
  * (32) JET collaboration, _Extracting the jet transport coefficient from jet quenching in high-energy heavy-ion collisions_ , _Phys. Rev. C_ 90 (2014) 014909 [1312.5003].
  * (33) H. Liu, K. Rajagopal and U.A. Wiedemann, _Calculating the jet quenching parameter from AdS/CFT_ , _Phys. Rev. Lett._ 97 (2006) 182301 [hep-ph/0605178].
  * (34) J. Casalderrey-Solana and D. Teaney, _Heavy quark diffusion in strongly coupled N=4 Yang-Mills_ , _Phys. Rev. D_ 74 (2006) 085012 [hep-ph/0605199].
  * (35) F. D’Eramo, H. Liu and K. Rajagopal, _Transverse Momentum Broadening and the Jet Quenching Parameter, Redux_ , _Phys. Rev. D_ 84 (2011) 065015 [1006.1367].
  * (36) S. Caron-Huot, _O(g) plasma effects in jet quenching_ , _Phys. Rev. D_ 79 (2009) 065039 [0811.1603].
  * (37) M. Panero, K. Rummukainen and A. Schäfer, _Lattice Study of the Jet Quenching Parameter_ , _Phys. Rev. Lett._ 112 (2014) 162001 [1307.5850].
  * (38) G.D. Moore, S. Schlichting, N. Schlusser and I. Soudi, _Non-perturbative determination of collisional broadening and medium induced radiation in QCD plasmas_ , _JHEP_ 10 (2021) 059 [2105.01679].
  * (39) Y.-T. Chien, A. Emerman, Z.-B. Kang, G. Ovanesyan and I. Vitev, _Jet Quenching from QCD Evolution_ , _Phys. Rev. D_ 93 (2016) 074030 [1509.02936].
  * (40) E. Bianchi, J. Elledge, A. Kumar, A. Majumder, G.-Y. Qin and C. Shen, _The $x$ and $Q^{2}$ dependence of $\hat{q}$, quasi-particles and the JET puzzle_, 1702.00481.
  * (41) C. Andrés, N. Armesto, M. Luzum, C.A. Salgado and P. Zurita, _Energy versus centrality dependence of the jet quenching parameter $\hat{q}$ at RHIC and LHC: a new puzzle?_, _Eur. Phys. J. C_ 76 (2016) 475 [1606.04837].
  * (42) J. Noronha-Hostler, B. Betz, J. Noronha and M. Gyulassy, _Event-by-event hydrodynamics $+$ jet energy loss: A solution to the $R_{AA}\otimes v_{2}$ puzzle_, _Phys. Rev. Lett._ 116 (2016) 252301 [1602.03788].
  * (43) J. Casalderrey-Solana, D.C. Gulhan, J.G. Milhano, D. Pablos and K. Rajagopal, _A Hybrid Strong/Weak Coupling Approach to Jet Quenching_ , _JHEP_ 10 (2014) 019 [1405.3864].
  * (44) D. Zigic, I. Salom, J. Auvinen, M. Djordjevic and M. Djordjevic, _DREENA-B framework: first predictions of $R_{AA}$ and $v_{2}$ within dynamical energy loss formalism in evolving QCD medium_, _Phys. Lett. B_ 791 (2019) 236 [1805.04786].
  * (45) C. Andres, N. Armesto, H. Niemi, R. Paatelainen and C.A. Salgado, _Jet quenching as a probe of the initial stages in heavy-ion collisions_ , _Phys. Lett. B_ 803 (2020) 135318 [1902.03231].
  * (46) A. Huss, A. Kurkela, A. Mazeliauskas, R. Paatelainen, W. van der Schee and U.A. Wiedemann, _Predicting parton energy loss in small collision systems_ , _Phys. Rev. C_ 103 (2021) 054903 [2007.13758].
  * (47) JETSCAPE collaboration, _Determining the jet transport coefficient $\hat{q}$ from inclusive hadron suppression measurements using Bayesian parameter estimation_, _Phys. Rev. C_ 104 (2021) 024905 [2102.11337].
  * (48) S. Platzer and M. Sjodahl, _The Sudakov Veto Algorithm Reloaded_ , _Eur. Phys. J. Plus_ 127 (2012) 26 [1108.6180].
  * (49) L. Lönnblad, _Fooling Around with the Sudakov Veto Algorithm_ , _Eur. Phys. J. C_ 73 (2013) 2350 [1211.7204].
  * (50) R. Kleiss and R. Verheyen, _Competing Sudakov Veto Algorithms_ , _Eur. Phys. J. C_ 76 (2016) 359 [1605.09246].
  * (51) R.K. Ellis, W.J. Stirling and B.R. Webber, _QCD and collider physics_ , vol. 8, Cambridge University Press (2, 2011), 10.1017/CBO9780511628788.
  * (52) Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * (53) A.H. Mueller and P. Nason, _Heavy particle content in QCD jets_ , _Phys. Lett. B_ 157 (1985) 226.
  * (54) A.H. Mueller and P. Nason, _HEAVY PARTICLE CONTENT IN QCD JETS_ , _Nucl. Phys. B_ 266 (1986) 265.
  * (55) P. Koch, B. Muller and J. Rafelski, _Strangeness in Relativistic Heavy Ion Collisions_ , _Phys. Rept._ 142 (1986) 167.
  * (56) A. Kurkela and A. Mazeliauskas, _Chemical Equilibration in Hadronic Collisions_ , _Phys. Rev. Lett._ 122 (2019) 142301 [1811.03040].
  * (57) A. Kurkela and A. Mazeliauskas, _Chemical equilibration in weakly coupled QCD_ , _Phys. Rev. D_ 99 (2019) 054018 [1811.03068].
  * (58) H.E. Haber, _Useful relations among the generators in the defining and adjoint representations of SU(N)_ , _SciPost Phys. Lect. Notes_ 21 (2021) 1 [1912.13302].
  * (59) L. Apolinário, J.G. Milhano, G.P. Salam and C.A. Salgado, _Probing the time structure of the quark-gluon plasma with top quarks_ , _Phys. Rev. Lett._ 120 (2018) 232301 [1711.03105].
  * (60) N. Armesto, H. Ma, Y. Mehtar-Tani, C.A. Salgado and K. Tywoniuk, _Coherence effects and broadening in medium-induced QCD radiation off a massive $q{\bar{q}}$ antenna_, _JHEP_ 01 (2012) 109 [1110.4343].
  * (61) Y. Mehtar-Tani, C.A. Salgado and K. Tywoniuk, _The Radiation pattern of a QCD antenna in a dense medium_ , _JHEP_ 10 (2012) 197 [1205.5739].
  * (62) T. Liou, A.H. Mueller and B. Wu, _Radiative $p_{\bot}$-broadening of high-energy quarks and gluons in QCD matter_, _Nucl. Phys. A_ 916 (2013) 102 [1304.7677].
  * (63) J.-P. Blaizot and Y. Mehtar-Tani, _Renormalization of the jet-quenching parameter_ , _Nucl. Phys. A_ 929 (2014) 202 [1403.2323].
  * (64) B. Wu, _Radiative energy loss and radiative $p_{\bot}$-broadening of high-energy partons in QCD matter_, _JHEP_ 12 (2014) 081 [1408.5459].
  * (65) P. Arnold and S. Iqbal, _The LPM effect in sequential bremsstrahlung_ , _JHEP_ 04 (2015) 070 [1501.04964].
  * (66) P. Arnold, T. Gorda and S. Iqbal, _The LPM effect in sequential bremsstrahlung: nearly complete results for QCD_ , _JHEP_ 11 (2020) 053 [2007.15018].
  * (67) K. Zapp, J. Stachel and U.A. Wiedemann, _A Local Monte Carlo implementation of the non-abelian Landau-Pomerantschuk-Migdal effect_ , _Phys. Rev. Lett._ 103 (2009) 152302 [0812.3888].
  * (68) K.C. Zapp, F. Krauss and U.A. Wiedemann, _A perturbative framework for jet quenching_ , _JHEP_ 03 (2013) 080 [1212.1599].
  * (69) K.C. Zapp, _JEWEL 2.0.0: directions for use_ , _Eur. Phys. J. C_ 74 (2014) 2762 [1311.0048].
  * (70) N. Armesto, L. Cunqueiro and C.A. Salgado, _Q-PYTHIA: A Medium-modified implementation of final state radiation_ , _Eur. Phys. J. C_ 63 (2009) 679 [0907.1014].
  * (71) B. Schenke, C. Gale and S. Jeon, _MARTINI: An Event generator for relativistic heavy-ion collisions_ , _Phys. Rev. C_ 80 (2009) 054913 [0909.2037].
  * (72) P. Caucal, E. Iancu, A.H. Mueller and G. Soyez, _A new pQCD based Monte Carlo event generator for jets in the quark-gluon plasma_ , _PoS_ HardProbes2018 (2019) 028 [1812.05393].
  * (73) J.H. Putschke et al., _The JETSCAPE framework_ , 1903.07706.
  * (74) L. Schulman, _Techniques and Applications of Path Integration_ , Dover Books on Physics, Dover Publications (2012).
  * (75) G. Moliere, _Theory of the scattering of fast charged particles. 2. Repeated and multiple scattering_ , _Z. Naturforsch. A_ 3 (1948) 78.
  * (76) F. D’Eramo, K. Rajagopal and Y. Yin, _Molière scattering in quark-gluon plasma: finding point-like scatterers in a liquid_ , _JHEP_ 01 (2019) 172 [1808.03250].
  * (77) B.G. Zakharov _Sov. J. Nucl. Phys._ 46 (1987) 92.
  * (78) Y.V. Kovchegov, _NonAbelian Weizsacker-Williams field and a two-dimensional effective color charge density for a very large nucleus_ , _Phys. Rev. D_ 54 (1996) 5463 [hep-ph/9605446].
  * (79) J. Jalilian-Marian, A. Kovner, L.D. McLerran and H. Weigert, _The Intrinsic glue distribution at very small x_ , _Phys. Rev. D_ 55 (1997) 5414 [hep-ph/9606337].
  * (80) R. Baier, A. Kovner and U.A. Wiedemann, _Saturation and parton level Cronin effect: Enhancement versus suppression of gluon production in p-A and A-A collisions_ , _Phys. Rev. D_ 68 (2003) 054009 [hep-ph/0305265].
  * (81) A. Kurkela and U.A. Wiedemann, _Picturing perturbative parton cascades in QCD matter_ , _Phys. Lett. B_ 740 (2015) 172 [1407.0293].


Entanglement in a two-qubit Heisenberg XXX model with x-components of
Dzyaloshinskii-Moriya and Kaplan-Shekhtman-Entin-Wohlman-Aharony interactions

Rachid<EMAIL_ADDRESS>, El Bouâzzaoui
<EMAIL_ADDRESS>, Mohammed El
<EMAIL_ADDRESS>, Abdelhadi
<EMAIL_ADDRESS>and Abdellatif
<EMAIL_ADDRESS>
1Team of Theoretical Physics and High Energy, Department of Physics, Faculty
of Sciences, Ibn Zohr University, Agadir, Morocco,

PO Box 8106, Agadir, Morocco

2Team of Theoretical Physics, Laboratory L.P.M.C., Department of Physics,
Faculty of Sciences, Chouaib Doukkali University, El Jadida, Morocco,

PO Box 20, 24000 El Jadida, Morocco

3ISPS2I Laboratory, National Higher School of Arts and Crafts (ENSAM), Hassan
II University, 20670, Casablanca, Morocco

This article explores the concurrence in a two-qubit Heisenberg XXX model with
Dzyaloshinskii-Moriya (DM) and Kaplan-Shekhtman-Entin-Wohlman-Aharony (KSEA)
interactions. The concurrence expression was developed using the physical
variables connected with the chosen system. Our results indicate that the
temperature, the spin coupling constant, the x-components of the DM and KSEA
interactions may all play a role in determining the degree of intricacy
between states. Additionally, these findings imply that the separability of
states is obtained at high-temperature domains or by switching the spin
coupling. In contrast, the entanglement of states may be achieved at low
temperatures or by using high values of the x-components of the DM and KSEA
parameters. Furthermore, the DM and KSEA interactions have an identical effect
on the concurrence behaviors at high temperatures.

PACS numbers: 03.67-a, 03.65.Ud, 75.10Jm

Keywords: entanglement,concurrence, Heisenberg chain, Dzyaloshinskii-Moriya
coupling, Kaplan-Shekhtman-Entin-Wohlman-Aharony coupling, thermal state.

##  1 Introduction

Because of its nonlocality [1, 2, 3], entanglements have become an important
tool in quantum information processing [4, 5], including quantum computing [6,
7], quantum cloning [8] and teleportation [9]. Using the Heisenberg model as a
basic system, entangled states may be generated and controlled. This model has
been used to simulate a wide range of physical phenomena, including quantum
dots [10], nuclear spins [11], optical lattices [12], and superconductors
[13]. Entanglement in a Heisenberg system with spin is more important because
of several magnetic properties of high-spin quantum systems, such as the
Haldane gap [14, 15]. It has been shown that the Heisenberg model in the two-
site system has a thermal entanglement for a given spin, such as the
spin-$1/2$ chain [16], spin-$1$ chain [17], and $(1/2,3/2)$ mixed-spin chain
[18].

Many quantum entanglement measurements have been established. The most
important of which is that of Bennett et al., who developed quantitative
measures of quantum entanglement in bipartite systems, first for pure states
[19] and then for hybrid states [20]. These measures utilize entropy and do
not use the extremization technique to get the optimum working ensemble.
Because extremization is notoriously tricky to deal with theoretically, the
authors could only obtain a precise formula for entanglement formation in
Bell-diagonal states [21]. Wootters then proved Hill and Wootters’ theory,
which offers a clear prescript for assessing the entanglement of any two-qubit
system [22]. The Wootters formula requires the solution of a fourth-degree
algebraic problem, requiring tedious Ferrari formulas. Thus, concise
formulations for multiple quantum states are required for practical reasons.

Researchers in quantum information have recently begun to investigate the
Dzyaloshinskii-Moriya (DM) interaction, which was discovered by Dzyaloshinskii
and Moriya and named after them. The scientific community has shown a
significant interest in the DM interaction, which is thought to be one of the
most effective regulating characteristics of quantum entanglement. Condensed
matter systems with correlated components have recently received much
interest. As a result, recent work has focused on quantifying and
characterizing quantum correlations in Heisenberg spin-chain models with DM
interaction [23, 24, 25]. The results show that the DM interaction may control
quantum correlations. The non-classical thermal correlations of the two-qubit
Heisenberg chain subjected to a transverse magnetic field and DM interaction
were also investigated using the local quantum uncertainty (LQU) [26].
Pourkarimi investigated the influence of the DM interaction on quantum discord
in the Heisenberg XY spin chain model [27]. His findings indicate that quantum
discord is more persistent than quantum entanglement in the presence of the DM
interaction. G. L. Kamta et al. explore the entanglement of a two-qubit
anisotropic Heisenberg XY chain in thermal equilibrium at temperature $T$ in
an external magnetic field $B$. They produce entanglement for any finite T
through the combined influences of anisotropic interactions and a magnetic
field B by adjusting its strength [28]. Sun et al. demonstrated that
increasing the DM interaction value increases the quantum discord of the two-
qubit Heisenberg XYZ spin chain model [29]. N. Canosa et al. examine the
entanglement of general mixed states of a two-qubit Heisenberg XYZ chain in a
magnetic field and its detection through different criteria. The exact
separability and weaker conditions implied by the disorder and the von Neumann
entropic criteria are analyzed [30].

In the fundamental study, [31, 32], it was discovered that there is symmetric
helical coupling associated with the DM interaction. It is ignored since the
symmetric helical interaction is significantly weaker than the DM interaction.
Kaplan discovers this symmetric helical interaction between local spins in the
single-band Hubbard model with spin-orbit couplings (SOC) [33]. Following
that, Shekhtman, Entin-Wohlman, and Aharony demonstrate that the mild
ferromagnetism of La2CuO4 may be explained by this non-zero symmetric helical
interaction [34]. As a result, the symmetric helical interaction was named the
Kaplan-Shekhtman-Entin-Wohlman-Aharony (KSEA) interaction. Thermal quantum
entanglement and discord have been investigated in a two-qubit XYZ Heisenberg
model with DM, and KSEA interactions in an inhomogeneous magnetic field [35].
The two-qubit Heisenberg XYZ model with the interactions of DM and KSEA in
thermal equilibrium was studied by A. V. Fedorova et al. [39]. They derived
analytical formulas for the LQU, and local quantum Fisher information (LQFI)
uses available expressions for the entropic quantum discord. By using the
correlation functions and chiral order factors, H. Fu et al. investigated the
impact of the KSEA interaction on the ground state features of three kinds of
spin chains in a transverse field by using the correlation functions and
chiral order factors [37].

Low-dimensional magnetism physics has recently focused on finding novel
magnetic compounds and refining their properties to fulfil emerging technology
demands. Knowing that integrable models of magnetism in low-dimensional
materials are based on the homogeneous Heisenberg XXX model and its
anisotropic version XXZ [38]. Motivated by the cited works and the utility of
Heisenberg XXX models, we investigated the thermal entanglement in a two-qubit
homogeneous Heisenberg XXX model with x-component of DM and KSEA interactions.

The organization of this paper is as follows. In section 2, in detail, we
present the theoretical model of our system subjected to DM and KSEA
interactions. In addition, we determine the ground states and their
entanglement at absolute zero temperature. Section 3 will be devoted to the
thermal quantum entanglement and concurrence expression with some limiting
cases. In section 4, numerical studies will be performed to highlight the
system behavior. Finally, the paper is closed with a findings overview and our
perspectives on the studied system.

##  2 Theoretical model

We consider a two-qubits Heisenberg XXX model with DM and KSEA interactions of
a one-half isotropic Heisenberg XXX model. The Hamiltonian of the system can
be written as:

$\mathcal{H}=\mathcal{H}_{H}+\mathcal{H}_{DM}+\mathcal{H}_{KSEA}$ (1)

where the first term refers to Heisenberg exchange couplings, the second one
denotes the DM interaction, and the third one represents the KSEA interaction.
Generally, the DM Hamiltonian $\mathcal{H}_{DM}$ [31, 32] can be expressed as
follows

$\displaystyle\mathcal{H}_{DM}$ $\displaystyle=$
$\displaystyle\bm{D}.\left(\bm{\sigma}_{1}\times\bm{\sigma}_{2}\right)$
$\displaystyle=$ $\displaystyle
D_{x}\left(\sigma_{1}^{y}\sigma_{2}^{z}-\sigma_{1}^{z}\sigma_{2}^{y}\right)+D_{y}\left(\sigma_{1}^{x}\sigma_{2}^{z}-\sigma_{1}^{z}\sigma_{2}^{x}\right)+D_{z}\left(\sigma_{1}^{x}\sigma_{2}^{y}-\sigma_{1}^{y}\sigma_{2}^{x}\right)$

where $\bm{D}=(D_{x},D_{y},D_{z})$ reflects the DM interaction vector and
$\sigma^{x,y,z}$ represents the Pauli spin matrices. However, the KSEA
Hamiltonian [33, 34] may be expressed as

$\displaystyle\mathcal{H}_{KSEA}$ $\displaystyle=$
$\displaystyle\bm{\sigma}_{1}.\bm{\Gamma}.\bm{\sigma}_{2}$ (3)
$\displaystyle=$
$\displaystyle(\sigma_{1}^{x},\sigma_{1}^{y},\sigma_{1}^{z})\left(\begin{array}[]{ccc}0&\Gamma_{z}&\Gamma_{y}\\\
\Gamma_{z}&0&\Gamma_{x}\\\ \Gamma_{y}&\Gamma_{x}&0\\\
\end{array}\right)\left(\begin{array}[]{c}\sigma_{2}^{x}\\\ \sigma_{2}^{y}\\\
\sigma_{2}^{z}\\\ \end{array}\right)$ $\displaystyle=$
$\displaystyle\Gamma_{x}\left(\sigma_{1}^{y}\sigma_{2}^{z}+\sigma_{1}^{z}\sigma_{2}^{y}\right)+\Gamma_{y}\left(\sigma_{1}^{x}\sigma_{2}^{z}+\sigma_{1}^{z}\sigma_{2}^{x}\right)+\Gamma_{z}\left(\sigma_{1}^{x}\sigma_{2}^{y}+\sigma_{1}^{y}\sigma_{2}^{x}\right)$

where $\bm{\Gamma}$ denotes the traceless symmetric tensor. Thus, both the DM
and KSEA interactions are generally composed of three components.  In the
present work, we thought to investigate the effect of x-component interactions
on thermal entanglements. To do this, we restrict ourselves to the case where
the DM and KSEA interactions occur along the x-axis. Then, our Hamiltonian can
be reads as

$\mathcal{H}=J\left(\sigma_{1}^{x}\sigma_{2}^{x}+\sigma_{1}^{y}\sigma_{2}^{y}+\sigma_{1}^{z}\sigma_{2}^{z}\right)+D_{x}\left(\sigma_{n}^{y}\sigma_{n+1}^{z}-\sigma_{n}^{z}\sigma_{n+1}^{y}\right)+\Gamma_{x}\left(\sigma_{n}^{y}\sigma_{n+1}^{z}+\sigma_{n}^{z}\sigma_{n+1}^{y}\right)$
(11)

where $J$ represents the real coupling constant for the spin interaction.
Antiferromagnetic chains are when the value of $J>0$, while ferromagnetic
chains are those in which the value of $J<0$. It should be noted that $D_{x}$
signifies the x-component of the DM interaction, while $\Gamma_{x}$ denotes
the x-component of the KSEA interaction. The Hamiltonian (11) can be expressed
in the standard computing basis $|00>$, $|01>$, $|10>$, $|11>$ as

$\mathcal{H}=\left(\begin{array}[]{cccc}J&-i\Gamma_{x}+iD_{x}&-i\Gamma_{x}-iD_{x}&0\\\
i\Gamma_{x}-iD_{x}&-J&2J&i\Gamma_{x}+iD_{x}\\\
i\Gamma_{x}+iD_{x}&2J&-J&i\Gamma_{x}-iD_{x}\\\
0&-i\Gamma_{x}-iD_{x}&-i\Gamma_{x}+iD_{x}&J\\\ \end{array}\right)$ (12)

The eigenvalue equation’s solution yields to the eigenvalues

$\displaystyle\epsilon_{1,2}$ $\displaystyle=$ $\displaystyle J\pm
2\Gamma_{x}$ (13) $\displaystyle\epsilon_{3,4}$ $\displaystyle=$
$\displaystyle-J\pm 2\eta$ (14)

where $\eta=\sqrt{D_{x}^{2}+J^{2}}$, and the related eigenvectors

$\displaystyle|\varphi_{1}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}|00\rangle-\frac{i}{2}|01\rangle-\frac{i}{2}|10\rangle+\frac{1}{2}|11\rangle$
(15) $\displaystyle|\varphi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}|00\rangle+\frac{i}{2}|01\rangle+\frac{i}{2}|10\rangle+\frac{1}{2}|11\rangle$
$\displaystyle|\varphi_{3}\rangle$ $\displaystyle=$
$\displaystyle-\frac{1}{\sqrt{2}}\sin(\theta_{1})|00\rangle-\frac{i}{\sqrt{2}}\cos(\theta_{1})|01\rangle+\frac{i}{\sqrt{2}}\cos(\theta_{1})|10\rangle+\frac{1}{\sqrt{2}}\sin(\theta_{1})|11\rangle$
$\displaystyle|\varphi_{4}\rangle$ $\displaystyle=$
$\displaystyle-\frac{1}{\sqrt{2}}\sin(\theta_{2})|00\rangle+\frac{i}{\sqrt{2}}\cos(\theta_{2})|01\rangle-\frac{i}{\sqrt{2}}\cos(\theta_{2})|10\rangle+\frac{1}{\sqrt{2}}\sin(\theta_{2})|11\rangle$

where $\theta_{1,2}$ are defined by

$\displaystyle\theta_{1,2}$ $\displaystyle=$
$\displaystyle\arctan\left(\frac{D_{x}}{\eta\mp J}\right)$ (16)

First, we address the entanglement of the system’s ground state at absolute
zero temperature, which is necessary before we can explore thermal quantum
entanglements. It is reasonable to suppose that the ground state entanglement
occurs in both the antiferromagnetic and ferromagnetic instances since the
energies are the exchange coupling $J$ functions.

In order to highlight the dependency of the ground state energy on KSEA
coupling constant $\Gamma_{x}$, we explore three significant cases. The first
one is $\Gamma_{x}>0$. Using equations (13) and (14), the ground state
energies can be written as

$\displaystyle\epsilon_{2}$ $\displaystyle=$ $\displaystyle
J-2\Gamma_{x},\qquad\text{if}\quad\Gamma_{x}>J+\eta$ (17)
$\displaystyle\epsilon_{4}$ $\displaystyle=$
$\displaystyle-J-2\eta,\qquad\text{if}\quad\Gamma_{x}<J+\eta$

Consequently, for $\Gamma_{x}>J+\eta$, the ground state it will be the
entangled state $|\varphi_{2}\rangle$, regardless of the sign of $J$, and when
$\Gamma_{x}<J+\eta$, the ground state is the entangled state
$|\varphi_{4}\rangle$ independently of the sign of $J$. Because
$|\varphi_{2}\rangle$ and $|\varphi_{4}\rangle$ are maximally entangled at low
temperatures, thus conducting to the maximum concurrence of $\mathcal{C}=1$ is
achieved.

In the second case, independently of the sign of $J$ for $\Gamma_{x}<0$, the
ground state energies are given by

$\displaystyle\epsilon_{1}$ $\displaystyle=$ $\displaystyle
J+2\Gamma_{x},\qquad\text{if}\quad\Gamma_{x}<-J-\eta$ (18)
$\displaystyle\epsilon_{4}$ $\displaystyle=$
$\displaystyle-J-2\eta,\qquad\text{if}\quad\Gamma_{x}>-J-\eta$

where $|\varphi_{1}\rangle$ and $|\varphi_{4}\rangle$ are the corresponding
maximally entangled states leading to the maximum concurrence $\mathcal{C}=1$.

However, in the last case where $\Gamma_{x}=0$, the ground state energies are
dependent on the sign of $J$. In addition, if $J<0$, those energies are given
by

$\displaystyle\epsilon_{1}$ $\displaystyle=$
$\displaystyle\epsilon_{2}=J,\qquad\text{if}\quad J<-\eta$ (19)
$\displaystyle\epsilon_{4}$ $\displaystyle=$
$\displaystyle-J-2\eta,\qquad\text{if}\quad J>-\eta.$

In contrast, for $J>0$, the ground state energy is given by

$\displaystyle\epsilon_{4}$ $\displaystyle=$
$\displaystyle-J-2\eta,\qquad\text{for any value of }\eta$ (20)

For the antiferromagnetic chains, the significant ground state is $|\varphi
4\rangle$. It was also the ground state for ferromagnetic chains when $J$
fulfilled the condition $J+\eta>0$, leading to maximal concurrence
$\mathcal{C}=1$. To complete this investigation, we illustrate the three cases
discussed above in the phase diagram shown in Fig. 1.

Figure 1: (Color online) Phase diagram at $T=0$ of a two qubits XXX Heisenberg
model. (a) $\Gamma_{x}\neq 0$ for $J=1$, (b) $\Gamma_{x}=0$.

##  3 Thermal quantum entanglements

After determining our system’s spectrum, using the density matrix $\rho(T)$ to
depict the system’s state on thermal equilibrium at a given temperature $T$.
The expression of $\rho(T)$ is presented by

$\rho(T)={1\over\mathbb{Z}}e^{-\beta\mathcal{H}}$ (21)

where

$\mathbb{Z}={\rm Tr}e^{-\beta\mathcal{H}}$ (22)

Where $\mathbb{Z}$ is the canonical ensemble partition function and
$\beta=1/k_{B}T$ is the inverse thermodynamic temperature, where $k_{B}$ is
the Boltzmann’s constant, which is treated as unity in the following for
simplicity. That will be succeeded by using the spectral decomposition of the
Hamiltonian (12), which allows the thermal density matrix $\rho(T)$ to be
represented as

$\rho(T)={1\over\mathbb{Z}}\sum_{l=1}^{4}e^{-\beta\epsilon_{l}}|\phi_{l}\rangle\langle\phi_{l}|$
(23)

The density matrix of the system, as discussed before in thermal equilibrium,
may be represented in the normal computational basis by inserting (13) and
(15) into equation (23) and obtaining

$\rho(T)={1\over\mathbb{Z}}\left(\begin{array}[]{cccc}a&i\mu&i\nu&c\\\
-i\mu&b&d&-i\nu\\\ -i\nu&d&b&-i\mu\\\ c&i\nu&i\mu&a\\\ \end{array}\right)$
(24)

the elements matrix is represented by the equations

$\displaystyle a$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(2e^{-\beta\epsilon_{3}}\sin^{2}\left(\theta_{1}\right)+2e^{-\beta\epsilon_{4}}\sin^{2}\left(\theta_{2}\right)+e^{-\beta\epsilon_{1}}+e^{-\beta\epsilon_{2}}\right)$
(25) $\displaystyle b$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(2e^{-\beta\epsilon_{3}}\cos^{2}\left(\theta_{1}\right)+2e^{-\beta\epsilon_{4}}\cos^{2}\left(\theta_{2}\right)+e^{-\beta\epsilon_{1}}+e^{-\beta\epsilon_{2}}\right)$
$\displaystyle c$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(-2e^{-\beta\epsilon_{3}}\sin^{2}\left(\theta_{1}\right)-2e^{-\beta\epsilon_{4}}\sin^{2}\left(\theta_{2}\right)+e^{-\beta\epsilon_{1}}+e^{-\beta\epsilon_{2}}\right)$
$\displaystyle d$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(-2e^{-\beta\epsilon_{3}}\cos^{2}\left(\theta_{1}\right)-2e^{-\beta\epsilon_{4}}\cos^{2}\left(\theta_{2}\right)+e^{-\beta\epsilon_{1}}+e^{-\beta\epsilon_{2}}\right)$
$\displaystyle\mu$ $\displaystyle=$
$\displaystyle-\frac{1}{4}\left(-e^{-\beta\epsilon_{3}}\sin\left(2\theta_{1}\right)+e^{-\beta\epsilon_{4}}\sin\left(2\theta_{2}\right)+e^{-\beta\epsilon_{1}}-e^{-\beta\epsilon_{2}}\right)$
$\displaystyle\nu$ $\displaystyle=$
$\displaystyle-\frac{1}{4}\left(e^{-\beta\epsilon_{3}}\sin\left(2\theta_{1}\right)-e^{-\beta\epsilon_{4}}\sin\left(2\theta_{2}\right)+e^{-\beta\epsilon_{1}}-e^{-\beta\epsilon_{2}}\right)$

Consequently, the partition function is clearly described by

$\mathbb{Z}=2e^{\beta J}\cosh\left(2\beta\eta\right)+2e^{-\beta
J}\cosh(2\beta\Gamma_{x})$ (26)

This matrix $\rho$ is examined in order to find out how much entanglement is
linked with it.

$R=\rho S\rho^{\ast}S$ (27)

where $\rho$ has a complex conjugate in the form of $\rho^{\ast}$, and $S$ is
determined by the formula

$S=\sigma^{y}\otimes\sigma^{y}$ (28)

where $\sigma^{y}$ denotes the Pauli matrix, and the $R$ matrix may be
obtained by a simple calculation by

$R=\left(\begin{array}[]{cccc}R_{11}&R_{12}&R_{13}&R_{14}\\\
R_{12}^{\ast}&R_{22}&R_{23}&R_{13}^{\ast}\\\
R_{13}^{\ast}&R_{23}&R_{22}&R_{12}^{\ast}\\\ R_{14}&R_{13}&R_{12}&R_{11}\\\
\end{array}\right)$ (29)

where the $R$ components matrix are denoted by

$\displaystyle R_{11}$ $\displaystyle=$ $\displaystyle
a^{2}+c^{2}+\mu^{2}+\nu^{2}$ (30) $\displaystyle R_{22}$ $\displaystyle=$
$\displaystyle b^{2}+d^{2}+\mu^{2}+\nu^{2}$ $\displaystyle R_{14}$
$\displaystyle=$ $\displaystyle 2ac+2\mu\nu$ $\displaystyle R_{12}$
$\displaystyle=$ $\displaystyle i(\mu(a+b)+\nu(c+d))$ $\displaystyle R_{13}$
$\displaystyle=$ $\displaystyle i(\nu(a+b)+\mu(c+d))$ $\displaystyle R_{23}$
$\displaystyle=$ $\displaystyle 2bd+2\mu\nu$

To grasp the meaning of $R$, consider that ${\rm Tr}R$, which ranges from $0$
to $1$, measures the degree of equality between $\rho$ and $\rho^{\ast}$,
which indicates how closely $\rho$ approximates a mixture of generalized Bell
states. Additionally, the eigenvalues of $R$ are invariant under local unitary
transformations of the individual qubits, which qualifies them to be included
in a formula for entanglement, as entanglement must also be invariant under
such transformations. Using (29) and (30), one can check quickly that the
square roots of the eigenvalues of the matrix $R$ are given by

$\displaystyle\lambda_{1}$ $\displaystyle=$
$\displaystyle\left(\frac{(a-c)^{2}+(b-d)^{2}+2(\mu-\nu)^{2}-(a+b-c-d)\sqrt{(a-b-c+d)^{2}+4(\mu-\nu)^{2}}}{2\mathbb{Z}^{2}}\right)^{1\over
2}$ (31) $\displaystyle\lambda_{2}$ $\displaystyle=$
$\displaystyle\left(\frac{(a-c)^{2}+(b-d)^{2}+2(\mu-\nu)^{2}+(a+b-c-d)\sqrt{(a-b-c+d)^{2}+4(\mu-\nu)^{2}}}{2\mathbb{Z}^{2}}\right)^{1\over
2}$ $\displaystyle\lambda_{3}$ $\displaystyle=$
$\displaystyle\left(\frac{(a+c)^{2}+(b+d)^{2}+2(\mu+\nu)^{2}-(a+b+c+d)\sqrt{(a-b+c-d)^{2}+4(\mu+\nu)^{2}}}{2\mathbb{Z}^{2}}\right)^{1\over
2}$ $\displaystyle\lambda_{4}$ $\displaystyle=$
$\displaystyle\left(\frac{(a+c)^{2}+(b+d)^{2}+2(\mu+\nu)^{2}+(a+b+c+d)\sqrt{(a-b+c-d)^{2}+4(\mu+\nu)^{2}}}{2\mathbb{Z}^{2}}\right)^{1\over
2}$

When attempting to quantify the amount of entanglement associated with $\rho$,
we take into consideration the concurrence [21, 22], which may be described as
follows:

$\mathcal{C}=\max\left[0,2\max\left(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}\right)-\sum_{i=1}^{4}\lambda_{i}\right]$
(32)

The expression mentioned above is significant in determining the extent of
entanglement between states; however, the concurrence values are included in
the range from zero, which represents an unentangled state, to one, which
signifies a maximally entangled state. After deriving the concurrence
expression, which is implicitly reliant on the coupling constant $J$, the
x-components of the DM and KSEA interactions, and the temperature $T$.
Consequently, we now have all of the components necessary to examine the
behavior of the proposed system concerning the previously specified
quantities. Before beginning the numerical section, it is needed to explore a
few limiting cases to simplify the expression of concurrence and understand
the behavior of our system.

The first limiting case concerns the high-temperature regime. To make things
easier, we will use the same formulas for the various expressions in the
previous section, assuming that $J=\Gamma_{x}=D_{x}=1$. The eigenvalues of the
matrix $R$ are represented in the following manner:

$\displaystyle\lambda_{1}$ $\displaystyle\simeq$
$\displaystyle\left(\frac{1}{8}\left(1-2\sqrt{2}\right)\beta+\frac{1}{16}\right)^{1\over
2}$ (33) $\displaystyle\lambda_{2}$ $\displaystyle\simeq$
$\displaystyle\left(\frac{1}{8}\left(2\sqrt{2}+1\right)\beta+\frac{1}{16}\right)^{1\over
2}$ $\displaystyle\lambda_{3}$ $\displaystyle\simeq$
$\displaystyle\left(\frac{1}{16}-\frac{3\beta}{8}\right)^{1\over 2}$ (34)
$\displaystyle\lambda_{4}$ $\displaystyle\simeq$
$\displaystyle\left(\frac{\beta}{8}+\frac{1}{16}\right)^{1\over 2}$

Following a comparison of the eigenvalues of $R$ at high temperature, the
expression of concurrence, in this case, maybe expressed by the relationship

$\mathcal{C}\simeq\max\left(0,\frac{1}{2}\left(\left(2\sqrt{2}+1\right)\beta-1\right)\right)$
(35)

When $T$ approaches a substantial value, the expression of concurrence reveals
that $\mathcal{C}\rightarrow 0$ makes the system states separable and not
entangled at high temperatures.

Strong exchange coupling $J$ is the second limiting case where we consider
that $J\gg\Gamma_{x}$ and $J\gg D_{x}$. As a result, the quantities
$\theta_{1}$ and $\theta_{2}$ have the values $\theta_{1}\simeq\pi/2$ and
$\theta_{2}\simeq D_{x}/2J$, respectively, and the eigenvalues of the matrix
$R$ are dependent on the signature of $J$. For an antiferromagnetic system
($J>0$), the eigenvalues of the matrix $R$ are represented by

$\displaystyle\lambda_{1}$ $\displaystyle\simeq$
$\displaystyle\frac{1}{\left(2\cosh(2\beta\Gamma_{x})+e^{4\beta J}+1\right)}$
(36) $\displaystyle\lambda_{2}$ $\displaystyle\simeq$
$\displaystyle\frac{e^{4\beta J}}{\left(2\cosh(2\beta\Gamma_{x})+e^{4\beta
J}+1\right)}$ $\displaystyle\lambda_{3}$ $\displaystyle\simeq$
$\displaystyle\frac{1}{\left(2\cosh(2\beta\Gamma_{x})+e^{4\beta J}+1\right)}$
$\displaystyle\lambda_{4}$ $\displaystyle\simeq$
$\displaystyle\frac{1}{\left(2\cosh(2\beta\Gamma_{x})+e^{4\beta J}+1\right)}$

and consequently, the expression of concurrence is given by

$\mathcal{C}\simeq\max\left(0,\frac{e^{4\beta
J}-3}{\left(2\cosh(2\beta\Gamma_{x})+e^{4\beta J}+1\right)}\right)$ (37)

If we look at the previous concurrence, it is easily possible to see that when
$J\rightarrow+\infty$ , the value of the concurrence is
$\mathcal{C}\rightarrow 1$, which indicates that the states are maximally
entangled. For $J\rightarrow 0$, it is clear that the concurrence tends to
$0$, which suggests that the state’s system has become more separable. For a
ferromagnetic system ($J<0$), the expression of concurrence is expressed by

$\mathcal{C}\simeq\max\left(0,-\frac{e^{4\beta
J}+1}{\left(2\cosh(2\beta\Gamma_{x})+e^{4\beta J}+1\right)}\right)$ (38)

The previous concurrence can have a value of $\mathcal{C}\rightarrow 1$ when
$J\rightarrow-\infty$, indicating that the states are most entangled. When
$J\rightarrow 0$, it is evident that the concurrence goes to $0$, meaning that
the state’s system has become more separable.

Finally, the last limiting case is studying strong DM and KSEA interactions.
In this perspective, we assume that $D_{x}\gg J$ and $\Gamma_{x}\gg J$, which
implies that the quantities $\theta_{1}$ and $\theta_{2}$ take on the values
of $\theta_{1}=\theta_{2}\simeq\pi/4$, and as a result, the eigenvalues of the
matrix $R$ are expressed in the following form

$\displaystyle\lambda_{1}$ $\displaystyle\simeq$ $\displaystyle\frac{e^{-\beta
D_{x}}}{\left(2\left(\cosh\left(2\beta
D_{x}\right)+\cosh\left(2\beta\Gamma_{x}\right)\right)\right)^{1\over 2}}$
(39) $\displaystyle\lambda_{2}$ $\displaystyle\simeq$
$\displaystyle\frac{e^{\beta D_{x}}}{\left(2\left(\cosh\left(2\beta
D_{x}\right)+\cosh\left(2\beta\Gamma_{x}\right)\right)\right)^{1\over 2}}$
$\displaystyle\lambda_{3}$ $\displaystyle\simeq$
$\displaystyle\frac{e^{-\beta\Gamma_{x}}}{\left(2\left(\cosh\left(2\beta
D_{x}\right)+\cosh\left(2\beta\Gamma_{x}\right)\right)\right)^{1\over 2}}$
$\displaystyle\lambda_{4}$ $\displaystyle\simeq$
$\displaystyle\frac{e^{\beta\Gamma_{x}}}{\left(2\left(\cosh\left(2\beta
D_{x}\right)+\cosh\left(2\beta\Gamma_{x}\right)\right)\right)^{1\over 2}}$

It is clear that the (39) depends on $\Gamma_{x}$ and $D_{x}$. Therefore, the
expression of concurrence is also dependent on the $\Gamma_{x}$ and $D_{x}$.
For $\Gamma_{x}>D_{x}$, the concurrence is expressed by

$\mathcal{C}\simeq\max\left(0,\frac{\sqrt{2}\left(\sinh\left(\beta\Gamma_{x}\right)-\cosh\left(\beta
D_{x}\right)\right)}{\left(\cosh\left(2\beta
D_{x}\right)+\cosh\left(2\beta\Gamma_{x}\right)\right)^{1\over 2}}\right)$
(40)

The last equation shows that the concurrence tends toward asymptotic value
$\mathcal{C}\simeq 1$ for strong KSEA interaction parameters $\Gamma_{x}$,
which means that the states are maximally entangled in this case. For
$\Gamma_{x}\rightarrow 0$, the concurrence tends to $0$, which indicates that
the states are separated for this set of limits. In the second case, for
$\Gamma_{x}<D_{x}$, the concurrence takes the form

$\mathcal{C}\simeq\max\left(0,\frac{\sqrt{2}\left(\sinh\left(2\beta
D_{x}\right)-\cosh\left(2\beta\Gamma_{x}\right)\right)}{\left(\cosh\left(2\beta
D_{x}\right)+\cosh\left(2\beta\Gamma_{x}\right)\right)^{1\over 2}}\right)$
(41)

When $D_{x}\rightarrow+\infty$, the concurrence takes the value
$\mathcal{C}\simeq 1$, which indicates that the states are entangled.
Inversely, when $D_{x}\rightarrow 0$, the concurrence approaches $0$,
reflecting states’ separability for this specific condition.

To illustrate the suggested model’s overall performance, we will devote the
following section to a numerical evaluation of the concurrence $\mathcal{C}$,
described below. Following that, we supplied a few plots in terms of the
parameters of the considered system, such as temperature, spin coupling
exchange, DM, and KSEA interactions, and we will have further discussions to
conclude this.

##  4 Numerical results

This section of the paper will quantitatively investigate several features of
entanglement in a two-qubit Heisenberg XXX chain with x-components of DM and
KSEA interactions. First, we will examine concurrence as a function of
temperature $T$ and the constant coupling exchange for the spin interaction
$J$ by setting the x-component of DM and KSEA interactions to the value
$D_{x}=\Gamma_{x}=1$. Second, we will set the coupling constants for the spin
coupling exchange and the x-component of the KSEA interaction to
$\Gamma_{x}=J=1$ and plot the concurrence $\mathcal{C}$ as a function of the
x-component of the DM interaction $D_{x}$. Finally, we investigate the
concurrence’s behavior regarding the x-component of the KSEA interaction with
fixed values $J=D_{x}=1$. For simplicity, the values $k_{B}=\hbar=1$ shall be
assumed.

Figure 2: (Color online) (a) The concurrence versus logarithmic temperature
for different $J$. (b) The concurrence versus of constant coupling for
different $T$ ($\Gamma_{x}=D_{x}=1$).

In Fig. 2, we display the concurrence $\mathcal{C}$ in terms of logarithmic
temperature $Log(T)$ for various coupling constants for the spin interaction
$J=$ 0.4, 0.6 and 0.8 (Fig. 2.a). And versus coupling constant for the spin
interaction $J$ for different values of temperature $T=$ 0.4, 0.6 and 0.8
(Fig. 2.b). From Fig. 2.a, one can observe threshold temperatures $T_{c}$
beyond which the entanglement disappears. The critical temperature $T_{c}$ is
also impacted by the parameter $J$, such that $T_{c}$ and $J$ grow in
proportion. Something more significant is that in the low-temperature and
precisely in the interval $0<Log(T)<0.3$, although the value of $J$ rises, the
entanglement stays equal to $1$, this means that in the vicinity of $T=0$, the
temperature is significant in front of the coupling. In addition, with the
increase in temperature, we can achieve long-lasting entanglement by
increasing the coupling between the spin (the blue curve).In Fig. 2.b, for
$J=0$, we have a degeneracy for the ground state seen in the line delimiting
separate ground states in Fig. 1.b. In this case, the ground state is a
statistical mixture of the two states $|\varphi_{2}\rangle$ and
$|\varphi_{4}\rangle$ which has zero concurrences. For small values of $J$,
the entanglement of the ferromagnetic and antiferromagnetic chains increases
in the same way. According to Fig. 2.b, we can see that the concurrence
$\mathcal{C}$ for an antiferromagnetic system tends towards one rapidly for
large values of $J$. In contrast, for a ferromagnetic chain and large values
of $|J|$, the concurrence tends towards one, but, in comparison to
antiferromagnetic systems, this progress is being made slowly. In conclusion,
our system’s states are more entangled at low temperatures or for large values
of $|J|$; but, at high temperatures or for $J=0$, the system’s states are more
separable.

Figure 3: (Color online) The concurrence versus x-component of DM interaction
for different $T$ ($J=1$ and $\Gamma_{x}=1$).

In Fig. 3, we plot the concurrence $\mathcal{C}$ in terms of x-component
$D_{x}$ at low temperature for $T=$ 1, 2 and 3 (a) and at high temperature for
$T=$ 10, 20 and 30 (b). From Fig. 3, the first observation is that the
concurrence is symmetrical for $D_{x}=0$. The concurrence tends to a fixed
value $\mathcal{C}=1$ for the high absolute value of the x-component of DM
interaction even as the temperature rises. Moreover, for $D_{x}=0$, we notice
that the $\mathcal{C}$ presents a minimum at low temperature, which disappears
with the decrease in temperature and vice versa, confirmed by the previous
study in section (3). In addition, we note that the entanglement increases
symmetrically, and with the decrease in temperature, we may get a short time
entanglement. Otherwise, the concurrence shows an interval T-dependent, where
$\mathcal{C}$ stays zero at high temperatures near $D_{x}=0$. Hence, an
increase in temperature indicates a widening of this interval. Then we can
conclude that the state’s system is influenced by the DM interaction, in which
the large values of $D_{x}$ make the system more entangled. At high
temperatures for small values of $D_{x}$, the system becomes less entangled,
and the system’s states become completely separable.

Figure 4: (Color online) The concurrence versus x-component of KSEA
interaction for different $T$ ($J=D_{x}=1$).

In Fig. 4, we display the concurrence $\mathcal{C}$ in terms of the
x-component of the KSEA interaction $\Gamma_{x}$ at low temperature for $T=$
1, 2 and 3 (a) and high temperature for $T=$ 10, 20, and 30 (b). In Fig. 4.b,
the first finding is that for $\Gamma_{x}=0$, the concurrence is symmetrical.
For the large absolute value of the x-component of the KSEA interaction at the
low temperature shown in Fig. 4.a, the concurrence tends to a constant value
of $\mathcal{C}=1$ even as the temperature increases, i.e., the system’s
ground state for $\Gamma_{x}>0$ is a maximally entangled state $\varphi_{2}$
(the yellow area in Fig. 1.a), whereas for $\Gamma_{x}<0$, the ground state is
$\varphi_{1}$ (the green area in Fig. 1.a), which is likewise a maximally
entangled state, implying that $\mathcal{C}=1$. Furthermore, for
$2<\Gamma_{x}<3$ ($-3<\Gamma_{x}<-2$), the ground state of the system is an
equal mixture of the doublet states $\varphi_{2}$ and $\varphi_{4}$
($\varphi_{1}$ and $\varphi_{4}$), seen in the lines delimiting distinct
ground states in Fig. 1.a, thus involving to $\mathcal{C}=0$. Additionally, we
see that entanglement rises symmetrically, and thus with a drop in
temperature, we may encounter a brief temporal entanglement. Otherwise, the
concurrence displays an interval in which $\mathcal{C}$ remains zero at high
temperatures near $\Gamma_{x}=0$, and the width of this interval is dependent
on the temperature. Thus a temperature rise suggests a lengthening of this
interval. In conclusion, the KSEA interaction influences the state’s system.
Large values of $\Gamma_{x}$ make the system more entangled, while small
values of $\Gamma_{x}$ make the system less entangled at high temperatures,
implying that the system’s states become separable. Furthermore, at high
temperatures, the DM and KSEA interactions have the same impact on the
concurrence ( Fig. 3.b and Fig. 4.b.)

##  5 Conclusion and perspectives

Concurrence, a measure of entanglement, is investigated in a two-qubit
Heisenberg XXX chain with x-components of DM and KSEA interactions. The
Hamiltonian model is given, and through mathematical calculations, the
eigenstates entanglement has been determined, and the thermal state at a
finite temperature is explicitly derived. The entanglement of the ground state
at zero temperature limits and the associated phase diagram has been
discussed. Afterwards, we obtained the concurrence expression, which depends
on the spin’s coupling exchange constant $J$, the DM and KSEA interactions’
x-components, and the temperature $T$. Subsequently, we have studied various
limits. Indeed, we analyzed the case at high temperature, the strong and weak
spin coupling exchange J compared to the couplings of DM and KSEA. The
behaviors of the concurrence-measured entanglements for our study have been
investigated numerically. In this paper, we have concluded that temperature
$T$, the x-components of the DM and KSEA interactions, may all play a role in
determining the degree of intricacy between the states to a greater or lesser
extent. Moreover, it is possible to infer from these results that the
separability of the states can be obtained in the high-temperature domains,
$J=0$, or at high temperature for small values of the constants of DM or KSEA
interactions. Again, the entanglement of the states can be obtained for large
values x-components of the DM and KSEA interactions or low temperature.
Moreover, DM and KSEA interactions similarly affect concurrence behaviors at
high temperatures.

Still some interesting questions to be addressed. Can we use the studied
system to investigate entanglement’s dynamic behaviors and describe the basic
features of the quantum entanglement at a finite time [40]? A related question
arose, what about other correlation measurements to study? These issues and
associated questions are under consideration.

## Acknowledgment

Our sincere gratitude goes to Mohamed Monkad, head of the Laboratory for
Physics of Condensed Matter (LPMC) in the Faculty of Sciences at Chouaïb
Doukkali University, for his essential assistance.

Funding: Not applicable.

Data availability: Not applicable.

Code availability: Not applicable.

## Declarations

Conflict of interest: The authors declare that they have no conflict of
interest.

## References

  * [1] E. Schrödinger, Proc. Camb. Phil. Soc. 31, 555 (1935).
  * [2] A. Einstein, B. Podolsky and N. Rosen, Physical Review. 47, 777 (1935).
  * [3] J. S. Bell, Physics, 1, 195–200 (1964).
  * [4] C. H. Bennett et al., Phys. Rev. Lett. 70, 1895 (1993).
  * [5] M. Murao et al. Phys. Rev. A, 59, 156 (1999).
  * [6] L. K. Grover, Phys. Rev. Lett. 80, 4329 (1998).
  * [7] L. K. Grover, Phys. Rev. Lett. 79, 325 (1997).
  * [8] G. D. Chiara et al., Phys. Rev. A 70, 062308 (2004).
  * [9] Y. Yeo, Phys. Rev. A 68, 022316 (2003).
  * [10] B. Trauzettel et al., Nat. Phys. 3, 192 (2007).
  * [11] B. E. Kane, Nature 393, 133 (1998).
  * [12] A. Sørensen and K. Mølmer, Phys. Rev. Lett. 83, 2274 (1999).
  * [13] M. Nishiyama et al., Phys. Rev. Lett. 98, 047002 (2007).
  * [14] R. M. Morra et al., Phys. Rev. B 38, 543 (1988).
  * [15] X. Wang , S. Qin and L. Yu, Phys. Rev. B 60, 14529 (1999).
  * [16] L. Zhou , H. S. Song et al., Phys. Rev. A 68, 024301 (2003).
  * [17] W. Song, Chin. Phys. B 18, 3251 (2009).
  * [18] K. T. Guo et al., J. Phys. A 43, 505301 (2010).
  * [19] C. H. Bennett, H. J. Bernstein, S. Popescu and B. Schumacher, Phys. Rev. A 53, 2046 (1996).
  * [20] C. H. Bennett, D. P. DiVincenzo, J.A. Smolin and W. K. Wootters, Phys. Rev. A 54, 3824 (1996).
  * [21] W. K. Wootters, Phys. Rev. Lett. 80, 2245 (1998).
  * [22] S. Hill and W. K. Wootters, Phys. Rev. Lett. 78, 5022 (1997).
  * [23] R. Jafari and A. Langari, Int. J. Quantum Inf. 9, 1057 (2011).
  * [24] R. Jafari, M. Kargarian and A. Langari, M. Siahatgar, Phys. Rev. B 78, 214414 (2008).
  * [25] E. Mehran, S. Mahdavifar and R. Jafari, Phys. Rev. A 89, 049903 (2014).
  * [26] N. Habiballah, Y. Khedif and M. Daoud, Eur. Phys. J. D 72, 154 (2018).
  * [27] M. R. Pourkarimi, Int. J. Theor. Phys. 57, 1158 (2018).
  * [28] G. L. Kamta and A. F. Starace, Phys. Rev. Lett. 88, 107901 (2002).
  * [29] Y. Sun Y, X. Ma , J. Guo, Quantum Inf. Process. 19, 1 (2020).
  * [30] N. Canosa and R. Rossignoli, Phys. Rev. A 69, 052306 (2004).
  * [31] T. Moriya, Phys. Rev. 120, 91 (1960).
  * [32] T. Moriya, Phys. Rev. Lett. 4, 228 (1960).
  * [33] T. A. Kaplan, Z. Phys. B, 49, 313 (1983).
  * [34] L. Shekhtman, O. Entin-Wohlman, A. Aharony, Phys. Rev. Lett. 69, 836 (1992).
  * [35] M. A Yurischev, Quantum Inf. Process. 19, 336 (2020).
  * [36] A. V. Fedorova and M. A. Yurischev, Quantum Inf. Process. 20, 169 (2021).
  * [37] H. Fu, M. Zhong and P. Tong, Europ. Phys. J. B 93, 80 (2020).
  * [38] M. Motamedifar, Phys. A Stat. Mech. Appl. 483, 280 (2017).
  * [39] A. V. Fedorova and M. A. Yurischev, Quantum Inf. Process. 20, 169 (2021).
  * [40] C. A. G. Gutiérrez and R. R. Ancheyta, Int. J. Quantum Inf. 14, 1650031 (2016).