Collective modes in anisotropic systems

Margaret E. Carrington
Department of Physics and Astronomy, Brandon University,
Brandon, Manitoba R7A 6A9, Canada
Winnipeg Institute for Theoretical Physics, Winnipeg, Manitoba, Canada

Bailey M. Forster
Department of  Physics and Astronomy, Brandon University,
Brandon, Manitoba R7A 6A9, Canada

Sofiya Makar
Department of Mathematics, Brandon University,
Brandon, Manitoba R7A 6A9, Canada

We study  collective modes in anisotropic plasmas of quarks and gluons using a quasi-particle picture and a hard loop approximation.
We use a general class of anisotropic distribution functions, and we consider chirally asymmetric systems.
We introduce a complete tensor basis to decompose the gluon polarization tensor into a set of nine scalar functions.
We derive and solve the corresponding dispersion equations.
Imaginary modes are particularly important because of their potential influence on plasma dynamics.
We explore in detail their dependence on the chiral chemical potential and the parameters that characterise the anisotropy of the system.
We show that our generalized distributions produce dispersion relations that are much richer in structure than those obtained with a simple one parameter deformation of an isotropic distribution.
In addition, the size and domain of the imaginary solutions are enhanced, relative to those obtained with a one parameter deformation.
Finally, we show that the influence of even a very small chiral chemical potential is significantly magnified when anisotropy is present.

§ INTRODUCTION

We consider a relativistic plasma of charged particles that is not necessarily in equilibrium, but is in the regime where a quasi-particle description is valid.
We use an effective phase-space description in terms of distribution functions that are not assumed thermal or isotropic.
A convenient parametrization for the distribution function of a spheroidially anisotropic system was introduced in [1, 2].
The idea is to introduce one parameter that effectively stretches or squeezes the spherically symmetric isotropic distribution in one direction.
This parametrization is interesting in the context of heavy ion collisions where one expects the parton distribution that is produced by the collision to be squeezed along the beam axis.
Dispersion relations are obtained from the poles of the dressed propagators of the partonic quasi-particles.
One finds that the spectrum of collective excitations is much richer than the spectrum of a system in thermal equilibrium.
It is especially interesting that there are imaginary modes, which are associated with instabilities and could play an important role in the thermalization of the plasma [3, 4, 5, 6, 7, 8].
Anisotropic plasmas have been studied intensively in the context of quark-gluon plasmas, which thermalize earlier than expected, for reasons that are as yet unclear.
Collective modes have been calculated for both gluons [1, 2, 9] and quarks [10].
An ellipsoidal generalization was developed in [11, 12].
These distributions have been used to study various properties of quark-gluon plasma (QGP) including
transport coefficients [13, 14, 15, 16],
heavy quark bound states [17, 18, 19, 20],
and photon and dilepton production [10, 21, 22, 24, 23].

Isotropic plasmas also support imaginary collective modes, if there is a chiral imbalance [25].
Chiral plasmas with spheroidal anisotropy were studied in Ref. [26].
Chiral systems are of interest in a wide variety of different situations in particle physics, nuclear physics, condensed matter physics, and cosmology.
One important example is the chiral magnetic effect (CME) which is the production of a parity violating current in a plasma, in the presence of a magnetic field and an asymmetry between left and right handed fermions [27].

A chirally asymmetric plasma can be described  using a chiral chemical potential defined as $\mu_5 \equiv (\mu_R-\mu_L)/2$, where $\mu_L$ and $\mu_R$
are the chemical potentials of the left and right handed fermions.

The CME is particularly interesting in the context of heavy-ion collisions because it has been argued that the quark gluon plasma that is produced in non-central collisions may contain regions where $\mu_5$ is locally finite [28, 30, 31, 32, 33].

Instabilities in chiral plasmas have also been studied in electroweak theory at large lepton chemical potential [34] and in the context of the early universe at $T\gg\mu_5$ [35].

Collective modes can be studied using either kinetic theory, or effective field theories.
The equivalence of the kinetic theory method and a hard thermal loop (HTL) effective theory approach, for isotropic and chirally symmetric systems, was shown in [36, 37, 38].
For anisotropic chirally symmetric systems, the connection between kinetic theory and a generalization of the HTL effective theory which is called a hard loop (HL) effective theory has also been established [39, 40].
Isotropic chiral kinetic theories have been developed in [41, 42, 43, 44, 45], and applied to a variety of physical problems, see for example [46, 47, 48].

In this paper we use anisotropic distributions that are more general than those used in  previous works, and we study the spectrum of collective modes in anisotropic plasmas with non-zero chiral chemical potential.
The method can be applied generally to either a QED plasma of ultrarelativistic electrons and
positrons or a QGP.

We use natural units where $\hbar = c =1$. The indices $i,j,k = 1, 2, 3$ and $\mu, \nu = 0, 1, 2, 3$ label, respectively, the Cartesian spatial coordinates and those of Minkowski space. Our metric is mostly minus $g_{\mu\nu} = (1,-1,-1,-1)_{\rm diag}$.
We use capital letters for four-vectors so that, for example, $P^2=p_0^2-\vec p\cdot\vec p =  p_0^2-p^2$. We define the unit vector $\hat p = \vec p/|\vec p|$ and we will also use $\hat p_0 = p_0/p$. Real solutions to dispersion equations will be denoted $\omega(\vec p)$ and imaginary solutions are written $i\gamma(\vec p)$.

§ ISOTROPIC FORMALISM

In vacuum the photon polarization tensor can be written in terms of one scalar function using a transverse projection operator as $\Pi^{\mu\nu}(p_0,\vec p) = (P^2 g^{\mu\nu}-P^\mu P^\nu)\,\Pi(p_0,\vec p)$.

At finite temperature the rest frame of the heat bath, which we define with the four-vector $n^\mu = (1,0,0,0)$, breaks Lorentz invariance.
Time-like axial gauge (TAG) is  particularly useful at finite temperature, because the gauge condition is imposed in the heat bath rest frame.
In TAG we only have to consider the spatial components of the propagator, for which the Dyson equation has the form
D^ij -1(p_0,p⃗) = D_0^ij -1(p_0,p⃗) -Π^ij(p_0,p⃗) = P^2δ^ij + p^2 p̂^i p̂^j  -Π^ij(p_0,p⃗) .

In an equilibrated chirally symmetric plasma, the polarization tensor is symmetric and has two independent components, one transverse and the other longitudinal with respect to the three-vector $\vec p$.

An isotropic system with finite chiral chemical potential was first considered in [25], and we review the formalism presented in that paper below.
In the presence of finite $\mu_5$ the polarization tensor has an additional asymmetric component and can be written
Π^ij(p_0,p⃗) = (δ^ij-p̂^ip̂^j) Π_T(p_0,p⃗) + p̂^ip̂^j Π_L(p_0,p⃗) + iϵ^ijmp̂^m Π_A(p_0,p⃗)  .
The propagator, obtained from inverting the Dyson equation, is
D^ij(p_0,p⃗)    =    (δ^ij-p̂^ip̂^j) (P^2-Π_T(p_0,p⃗))/(P^2-Π_T(p_0,p⃗))^2 - Π_A^2(p_0,p⃗) + p̂^ip̂^j1/p_0^2-Π_L(p_0,p⃗)

+    iϵ^ijmp̂^m Π_A(p_0,p⃗)/(P^2-Π_T(p_0,p⃗))^2 - Π_A^2(p_0,p⃗) .

The dispersion relations are obtained from the dispersion equations which give the poles of the retarded propagator. From equation (<ref>) we have that the dispersion equations for an isotropic chiral plasma are
p_0^2-Π_L(p_0,p⃗) = 0

P^2-(Π_T(p_0,p⃗) + Π_A(p_0,p⃗) ) = 0

P^2- (Π_T(p_0,p⃗) - Π_A(p_0,p⃗) ) = 0  .
Throughout this paper we refer always to the retarded polarization tensor, and we will not include any subscripts to indicate this.

Our notation for the equilibrium distribution function is (see equation (<ref>))
n(k) =n^+(k) =  1/e^β(k-μ)+1   and   n̅(k) = n^-(-k) =  1/e^β(k+ μ)+1  .
The distributions for right/left handed particle/anti-particles are written $n_R(k)$, $\bar n_R(k)$, $n_L(k)$ and $\bar n_L(k)$, where, for example,  $n_R(k) = n(k)\big|_{\mu=\mu_R}$.

We calculate the 1-loop photon polarization tensor in the HTL approximation (see Appendix <ref> for details).
The contribution from right handed fermions is
Π^ij_R(p_0,p⃗) = Π^ij_Reven(p_0,p⃗) + Π^ij_Rodd(p_0,p⃗)

Π^ij_Reven(p_0,p⃗) = 2g^2∫d^3 k/(2π)^3n_R(k)+n̅_R(k)/k
(δ^ij+v^i p^j + p^i v^j/P·V +iϵ- P^2 v^i v^j/(P·V+iϵ)^2)

Π^ij_Rodd(p_0,p⃗) =  i g^2 P^2 ϵ^ijm ∫d^3 k/(2π)^3n_R(k) - n̅_R(k)/k^2
(p_0 v^m - p^m)/(P·V +iϵ)^2
and the corresponding expression for left handed fermions is obtained from the transformation $\mu_R\to \mu_L$.
We use $P\cdot V = p_0-\vec p\cdot\vec v$ and
$\vec k/\sqrt{k^2+m^2} \approx \hat k \equiv \vec v$ since massless fermions are consistent with the hard loop approximation.
It is straightforward to modify the HTL integrals to the case of a QCD plasma by changing the QED coupling constant to the QCD one, and using [49]
n(k)+n̅(k) →N_f/2(n_q(k)+n̅_q(k)) + N_c n_g(k)

n(k)-n̅(k) →N_f/2(n_q(k)- n̅_q(k))
where the subscripts $q$ and $g$ refer to quark (or anti-quark) and gluon distribution functions.

Equations (<ref>, <ref>) can be rewritten in a form that is sometimes more useful by integrating by parts.
A straightforward calculation gives
Π^ij_Reven(p_0,p⃗)    =     - 2 g^2 ∫d^3 k/(2π)^3 ∂(n_R(k)+n̅_R(k))/∂k^m  v^i (δ^jm + v^j p^m/P·V )

Π^ij_Rodd(p_0,p⃗)    =     - i g^2 ϵ^ijm ∫d^3k/(2π)^3 1/k ∂/∂k^l(n_R(k) - n̅_R(k))
(p_0 δ^lm + (p_0 v^m-p^m)p^l/P·V
) .

From the definitions in equations (<ref>, <ref>) one finds
Π^ij(p_0,p⃗)    =    1/2(Π^ij_R(p_0,p⃗) + Π^ij_L(p_0,p⃗))

=    1/2(Π^ij_Reven + Π^ij_Leven) + 1/2(Π^ij_Rodd + Π^ij_Lodd)

=    1/2(Π^ij_Reven + Π^ij_Reven|_μ_R →μ_L) + 1/2(Π^ij_Rodd + Π^ij_Rodd|_μ_R →-μ_L)  .

The functions $\Pi_T(p_0,\vec p)$, $\Pi_L(p_0,\vec p)$ and $\Pi_A(p_0,\vec p)$ in equations (<ref>, <ref>) are calculated by applying the projection operators

${\cal P}_T^{ij} = (\delta^{ij}-\hat p^i\hat p^j)/2$, ${\cal P}_L^{ij} = \hat p^i\hat p^j$, and ${\cal P}_A^{ij} = -i\epsilon^{ijm}\hat p^m/2$ to

$\Pi^{ij}(p_0,\vec p)$.
The resulting expressions for the three scalar components of the polarization tensor are
Π_T(p_0,p⃗) = m_D^2 p_0^2/2p^2 (1 - P^2/2p_0 p ln(p_0+p+iϵ/p_0- p+iϵ) )

Π_L(p_0,p⃗) = - m_D^2 p_0^2/p^2( 1- p_0/2p ln(p_0+p+iϵ/p_0- p+iϵ) )

Π_A(p_0,p⃗) = -g^2 μ_5 P^2 /2π^2 p (1-p_0/2p ln(p_0+p+iϵ/p_0- p+iϵ) )
where we have used $\mu_5 = (\mu_R-\mu_L)/2$ and defined the Debye mass parameter
m_D^2    =    2g^2 ∫d^3k/(2π)^3 1/k ([n(k)+n̅(k)]_μ_R + [n(k)+n̅(k)]_μ_L )  = g^2 (T^2/3 + 1/2π^2(μ_R^2+μ_L^2)) .

The transverse and longitudinal components of the polarization tensor are the familiar HTL results, and in the chirally symmetric limit $\mu_R = \mu_L$ gives $\mu_5=0$ and therefore $\Pi_A(p_0,\vec p)=0$.
Throughout the rest of this paper we use the shorthand notation $\hat\mu_5 = g^2\mu_5/(2\pi^2)$.
In all of our numerical calculations we set the Debye mass to one, which is equivalent to defining all dimensionful quantities in units of the Debye mass.

The dispersion relations are obtained by solving equations (<ref>-<ref>)
where the functions $\Pi_L(p_0,\vec p)$, $\Pi_T(p_0,\vec p)$ and $\Pi_A(p_0,\vec p)$ are given in equation (<ref>).
There are pure real and pure imaginary solutions that appear in positive and negative pairs.

Equation (<ref>) has two real solutions which are just the usual longitudinal HTL modes, and are written $\pm\omega_L(p)$.
Equation (<ref>) also has two real solutions, which we denote $\pm\omega_+(p)$.
Equation (<ref>) has two real solutions, which we call $\pm\omega_-(p)$, and two extra solutions that appear at values of the wave vector below the critical value $p_{\rm crit}=\hat\mu_5$.
These solutions are pure imaginary and denoted $\pm i\gamma_-(p)$.

The critical wave vector at which the imaginary solutions appear can be found with a Nyquist analysis (see  Appendix <ref> for details).
When $\hat\mu_5=0$ we have $\omega_+(p) = \omega_-(p) = \omega_T(p)$ where $\omega_T(p)$ is the transverse HTL mode.

In figure <ref> we show solutions to equations (<ref> - <ref>) for different values of the parameter $\hat\mu_5$.

Real ($\omega$) and imaginary ($\gamma$) solutions to the dispersion equations (<ref>-<ref>) for different values of $\hat\mu_5 = g^2\mu_5/(2\pi^2)$. The units are defined by setting $m_D=1$. The diagonal gray line in the top figure is the light-cone.

§ ANISOTROPIC FORMALISM

§.§ Distribution function

Next we consider the situation that the momentum distribution is not isotropic.
A simple way to account for momentum anisotropy is to replace the distribution function in equations (<ref>, <ref>) using $n(k) \to n(\vec k) = C_\xi \, n\big(k H_\xi(\vec v)\big)$ where $\vec v = \vec k/k$, and similarly for $\bar n(k)$. The subscript $\xi$ indicates dependence on a set of anisotropy parameters that can be used to construct a distribution that is deformed relative to the isotropic one. The factor $C_\xi$ is a normalization and can be defined in different ways depending on the calculation being done. Our choice is explained at the end of this section.

We can construct a completely general expression for the function $H_\xi(\vec v)$ as a sum of terms that are products of anisotropy parameters and dot products of the vector $\vec v$ with two perpendicular unit vectors, which we will write $\hat n_3$ and $\hat n_1$.
We are primarily interested in distributions of partons produced in heavy-ion collisions. In this context, we will take $\hat n_3$ to define the beam axis and $\hat n_1$ to give the direction of transverse anisotropy.
We restrict to functions that satisfy the condition $H_\xi(\vec v) = H_\xi(-\vec v)$ and use an expression of the form
H^2_ξ(v⃗)     =    ξ_0
+ξ_2(n⃗_1·v⃗)^2+ξ_9 (n⃗_3·v⃗)^2
+ξ_6 (n⃗_1·v⃗)(n⃗_3·v⃗)

+ξ_4 (n⃗_1·v⃗)^4
+ ξ_8 (n⃗_1·v⃗)^3 (n⃗_3·v⃗)
+ξ_11(n⃗_1·v⃗)^2 (n⃗_3·v⃗)^2
+ξ_13 (n⃗_1·v⃗)(n⃗_3·v⃗)^3
+ξ_14 (n⃗_3·v⃗)^4  . When $\xi_0=1$ and $\xi_{i\ne 0}=0$, we have $H^2_\xi(\vec v)=1$,  and the distribution is isotropic.
For an arbitrary choice of the anisotropy parameters $\xi_i$, the isotropic distribution is expanded in the direction of $\vec v$ if $H_\xi(\vec v)<1$, and contracted if $H_\xi(\vec v)>1$.
The values of the anisotropy parameters must be chosen so that $H^2_\xi(\vec v)$ is positive for all orientations of the vector $\vec v$,
which is equivalent to the requirement that $H_\xi(\vec v)$, and therefore the argument of the distribution function, is real and positive.
In section <ref> we will show that our expression (<ref>) is equivalent to including spherical harmonics up to second order.

Equation (<ref>) includes the parametrizations used in several previous calculations.
In the original work of Ref. [1] the authors used $\xi_0=1$ and $\xi_9\in(-1,\infty)$ with  all other parameters set to zero.
The chosen value of $\xi_9$ allows one to consider either a slightly prolate momentum distribution ($-1<\xi_9<0$), which is elongated along the beam axis, or an arbitrarily oblate distribution ($\xi_9>0$), which  is squeezed in the direction of the beam axis.
In the calculations of Ref. [11, 12], where ellipsoidal asymmetry was included for the first time,
the authors used $\xi_0=1$ and various values of $\xi_2$ and $\xi_9$ with all other parameters zero.

In Ref. [9] the authors considered a one parameter deformation using $\xi_0=1+\sigma$
and $\xi_9=-\sigma$ with $\sigma\in(-1,\infty)$. The resulting distribution is slightly oblate for $-1<\sigma<0$ and prolate to any degree for $\sigma>0$.

Our definition (<ref>) can be extended by including terms with higher order products of the scalar products $(\hat n_3\cdot\vec v)$ and $(\hat n_3\cdot\vec v)$.
The only restriction we impose, besides positivity, is the condition $H_\xi(\vec v) = H_\xi(-\vec v)$. This means there are no terms where the sum of the exponents is odd,  like for example $(\hat n_3\cdot\vec v)(\hat n_1\cdot\vec v)^2$.
The reason we require that $H_\xi(\vec v)$ is even under the transformation $\vec v \to -\vec v$ is as follows.
We said above that when general distributions are considered the polarization tensor is obtained from equations (<ref>, <ref>) by making the replacement $n(k) \to n(\vec k) = C_\xi \, n\big(k H_\xi(\vec v)\big)$. In fact, there is an extra contribution to the anisotropic version of equation (<ref>) of the form
Π^ij_extra(p_0,p⃗) ∼g^2∫d^3 k/(2π)^3 (n_R(k⃗)+n̅_R(k⃗))
(v^i v^j/p_0 - p⃗·v⃗ +iϵ-v^i v^j/p_0 + p⃗·v⃗ +iϵ )
that integrates to zero for any distribution that is even under the transformation $\vec k \to -\vec k$.
If the condition $H_\xi(\vec v) = H_\xi(-\vec v)$ is not satisfied, equation (<ref>) would produce a non-zero contribution to the polarization tensor that dominates over the HL terms, and is not present in the result obtained from semiclassical kinetic theory (which is equivalent to the HL expression).

We will calculate the anisotropic polarization tensor using equations (<ref>, <ref>) with the replacement $n(k) \to C_\xi \,n\left(k H_\xi(\vec v)\right)$ and similarly for $\bar n(k)$.
We define
$\tilde k = k H_\xi(\vec v)$
and perform a straightforward change of variable so that the integral over $\vec k$ is written as the product of an integral over the new variable $\tilde k$ and a factored integral involving the angular variables.
Equations (<ref>, <ref>) become
Π^ij_Reven(p_0,p⃗) =   (m_D^2)_R    C_ξ ∫dΩ/4π 1/H_ξ^4(v⃗) v^i M^l(Ω)     (δ^jl + p^l v^j/P·V + iϵ)

Π^ij_Rodd(p_0,p⃗) =   i g^2 μ_R/2π^2 ϵ^ijm p   C_ξ  ∫dΩ/4π 1/H_ξ^3(v⃗) M^l(Ω)
(p_0 δ^lm + (p_0 v^m-p^m)p^l/P·V + iϵ

where we have defined
M^l(Ω) = (1/2k ∂(k^2 H_ξ^2(v⃗))/∂k^l)

and used $(m^2_D)_R = g^2(T^2/3+\mu_R^2/\pi^2)$.
The full expression for $\Pi^{ij}$ is obtained using equation (<ref>), as in the isotropic case.

We close this section by explaining how we define the function of the anisotropy parameters which we call $C_\xi$, which can be thought of as a normalization of the anisotropic distribution function.
This factor is determined by enforcing the condition
∫d^3k/(2π)^3 n_μ=0(k)/k ≡C_ξ∫d^3k/(2π)^31/k   n_μ=0(k H_ξ(v⃗))
from which we obtain
C_ξ^-1 = ∫dΩ/4π1/H^2_ξ(v⃗) .
This choice of normalization gives that when the chemical potentials are zero the mass parameter in equation (<ref>) is independent of the anisotropy parameters.
At zero chemical potential the entire spectrum of collective modes depends on this one scale (since we have assumed massless partons - see under equation (<ref>)), and it is therefore natural to adopt a normalization that leaves the Debye mass invariant under a change of the anisotropy parameters.

We also note that using equations (<ref>, <ref>, <ref>) it is easy to see that the even part of the polarization tensor is unchanged when the complete set of anisotropy parameter are scaled uniformly: $\xi_i \to \Lambda \xi_i$ where $\Lambda$ is any positive constant. This means that when we compare the spectrum of collective modes obtained from different choices of anisotropy parameters, we are comparing the effects of a specific deformation, and not the influence of an overall change of scale.

§.§ Projection operators

To decompose our general expression for the $3\times 3$ polarization tensor we need a complete set
of nine projection operators. This can be done most conveniently using three normalized unit vectors that are related to the plasmon momentum $\vec p$ and the two anisotropy vectors $\hat n_1$ and $\hat n_3$.
From these three-vectors we construct three ortho-normal vectors which we call $(\hat p, ~n_f,~m_F)$. We define
p̂ = p⃗/p

n_f = ñ_f/√(ñ_f·ñ_f)   with   ñ_f = n_3 - (n_3·p̂) p̂

m_F=m̃_F/√(m̃_F·m̃_F)   with   m̃_F = m̃_f - (n_f·m̃_f) n_f   and   m̃_f = n_1 - (n_1 ·p̂) p̂ .
Using these unit vectors the projection operators are defined as
P_1^ĳ = m_F^i m_F^j  ,   P_2^ĳ = p̂^i p̂^j  ,   P_3^ĳ = n_f^i n_f^j

P_4^ĳ = p̂^i n_f^j + n_f^i p̂^j  ,   P_5^ĳ = p̂^i m_F^j + m_F^i p̂^j  ,   P_6^ĳ =  n_f^i m_F^j + m_F^i n_f^j

P_7^ĳ =  n_f^i p̂^j - p̂^i n_f^j  ,  P_8^ĳ =  n_f^i m_F^j -  m_F^i n_f^j  ,  P_9^ĳ =  p̂^i m_F^j - m_F^i p̂^j  .

We note that $\hat p^i \hat p^j + n_f^i n_f^j + m_F^i m_F^j = \delta^{ij}$ and therefore we can write
Π^ij        = π_1(δ^ij-p̂^i p̂^j) +π_2 P_2^ij + (π_3 - π_1) P_3^ij+π_4 P_4^ij+π_5 P_5^ij+π_6 P_6^ij

≡α(δ^ij-p̂^i p̂^j) +βP_2^ij + γP_3^ij + δP_4^ij + λP_5^ij+τP_6^ij  .

For future reference we write down the operators that project out each component of the polarization tensor:
P_α^ij=P_1^ij  ,
P_β^ij=P_2^ij  ,

P_δ^ij=1/2P_4^ij  ,
P_λ^ij=1/2P_5^ij  ,
P_τ^ij = 1/2 P_6^ij  .
They satisfy conditions of the form $P^{ij}_m P^{ji}_n = {\cal T}_{mn}$ where the 9$\times$9 tensor ${\cal T}$ is

2 P_1     0     0     0     P_5     P_6     0     P_8     P_9

0     2 P_2     0     P_4     P_5     0     P_7     0     P_9

0     0     2 P_3     P_4     0     P_6     P_7     P_8     0

0     P_4     P_4     2 (P_2+P_3)     P_6     P_5     0     P_9     P_8

P_5     P_5     0     P_6     2 (P_1+P_2)     P_4     P_8     P_7     0

P_6     0     P_6     P_5     P_4     2 (P_1+P_3)     -P_9     0     -P_7

0     P_7     P_7     0     P_8     -P_9     -2 (P_2+P_3)     -P_5     P_6

P_8     0     P_8     P_9     P_7     0     -P_5     -2 (P_1+P_3)     -P_4

P_9     P_9     0     P_8     0     -P_7     P_6     -P_4     -2 (P_1+P_2)

) .

The projectors in equation (<ref>) form a complete basis, and therefore the polarization tensor can be decomposed as
Π^ij = π_1 P_1^ij+π_2 P_2^ij+π_3 P_3^ij+π_4 P_4^ij+π_5 P_5^ij+π_6 P_6^ij+i π_7 P_7^ij+i π_8 P_8^ij+ i π_9 P_9^ij

where we have omitted the functional arguments to shorten the notation.
A decomposition that is related to the first six of our projection operators was defined in Ref. [50].
The last three projection operators are anti-symmetric in their indices, and the corresponding dressing functions can only be non-zero if the chiral chemical potential is non-zero.
We invert the Dyson equation (<ref>) using (<ref>) and obtain the propagator
D D^ij    =    [ (p_0^2-π_2) (P^2-π_3)-π_4^2-π_7^2 ] P_1^ij
+ [ (P^2-π_1) (P^2-π_3)-π_6^2-π_8^2 ] P_2^ij

+    [ (p_0^2-π_2) (P^2-π_1)-π_5^2-π_9^2 ] P_3^ij
+ [ (P^2-π_1)π_4 +π_5 π_6+π_8 π_9 ] P_4^ij

+    [(P^2-π_3) π_5 +π_4 π_6+π_7 π_8  ] P_5^ij
+ [ (p_0^2-π_2) π_6 +π_4 π_5-π_7 π_9 ] P_6^ij

+    i [ (P^2-π_1) π_7 +π_5 π_8-π_6 π_9  ] P_7^ij
+ i [ (p_0^2-π_2) π_8 +π_5 π_7+π_4 π_9  ] P_8^ij

+    i [ (P^2-π_3)  π_9 -π_6 π_7+π_4 π_8 ] P_9^ij
where we have defined
D    =    (p_0^2-π_2) (P^2-π_1) (P^2-π_3)
-(P^2-π_1) (π_4^2+π_7^2)
-(p_0^2-π_2) (π_6^2+π_8^2)

- (P^2-π_3) (π_5^2+π_9^2)
-2(π_4 π_5 π_6-π_7 π_9 π_6+π_5 π_7 π_8+π_4 π_8π_9) .
The anisotropic dispersion relations are the solutions of the dispersion equation, ${\cal D}=0$.

§.§ Coordinate system

To make further process we must define a  coordinate sytem. The vector along the beam axis can be chosen as $\hat n_3' = (0,0,1)$ and the vector that defines the direction of asymmetry in the transverse plane can then be chosen as
$\hat n_1' = (1,0,0)$. The external momentum vector then has the general form
$\vec p^{\;\prime} = p(\sin\theta \cos\phi, \sin\theta \sin\phi, \cos\theta)$.
We note that using these coordinates it is easy to see that the parametrization in equation (<ref>) is equivalent to an expansion in spherical harmonics $Y_{lm}(\theta,\phi)$ with $l\in(0,4)$ and $m\in(-l,l)$.

We can define coordinates that are more convenient for calculational purposes by performing a rotation [12].
We rotate counter-clockwise about the $z$-axis with angle $\theta_z = \pi/2-\phi$ and then counter-clockwise about the $x$-axis with $\theta_x=\theta$.
The result is that the basis vectors become
n̂_1 =  (sin(ϕ), cos(θ)cos(ϕ), sin(θ)cos(ϕ))

n̂_3 = (0, -sin(θ), cos(θ))

p̂ = p⃗ /p= (0,0,1)  .
Any scalar or pseudo-scalar quantity will give the same result when calculated with either the original or the rotated coordinate system.
This can be easily be seen from the fact that the scalar products $(\hat p\cdot\hat n_1)$, $(\hat p\cdot\hat n_3)$ and $(\hat n_1\cdot\hat n_3)$, and the pseudo-scalar $\epsilon^{ijk}\hat n_1^i \hat n_3^j \hat p^k$, are invariant. Likewise the dispersion relations are the same in either coordinate system.
Another useful thing about the coordinates in equation (<ref>) is that the isotropic limit has a simple form.
It is easy to show using equation (<ref>) that
$\delta^{ij}-\hat p^i\hat p^j = P_1^{ij}+P_3^{ij}$, $\hat p^i\hat p^j = P_2^{ij}$ and $\epsilon^{ijm}p^m = P^{ij}_8$.
In the isotropic limit there are only three different non-zero components in equation (<ref>) which are
$\pi_2=\Pi_L$ and
$\pi_8 = \Pi_A$.
Using these results one can see that equations (<ref>, <ref>) reduce to (<ref>).

The variables in the angular integrations in equations (<ref>, <ref>) are independent variables, and will be written in component form as
v⃗ = (sin(θ')cos(ϕ'),sin(θ')sin(ϕ'),cos(θ')) .
We will use the notation $x=\cos(\theta)$ and $x'=\cos(\theta')$.

§.§ Anisotropy parameters

The anisotropy parameters $\xi_9$ and $\xi_2$ produce uniform stretching or squeezing of the isotropic distribution along the beam axis (in the case of $\xi_9$), and along the vector $\hat n_1$ in the transverse plane (for $\xi_2$).
We have introduced additional anisotropy parameters that can be used to produce distribution functions that are deformed relative to the isotropic distribution in more general ways, and produce dispersion relations with more structure.
As discussed in section <ref>, we consider only distributions that satisfy the condition $f(\vec k) = f(-\vec k)$.
This restricts the form of the terms in equation (<ref>).
Each term is chosen so that $e_1+e_3$ is even, where the exponent of the factor $(\vec n_1\cdot\vec v)$ is denoted $e_1$ and the exponent of the factor $(\vec n_3\cdot\vec v)$ is called $e_3$.
There are two different possibilities: $e_1$ and $e_3$ individually even, or $e_1$ and $e_3$ individually odd.

The elliptically anisotropic distribution considered in [12] depends only on the two parameters $\xi_9$ and $\xi_2$, both of which have $e_1$ and $e_3$ individually even.

Any distribution for which $e_1$ and $e_3$ are individually even for all terms in (<ref>) can be mapped onto an elliptically anisotropic distribution with angularly dependent parameters.

We consider for example
the distribution constructed with the set of anisotropy parameters $\xi_0=1$, $\xi_4=20$, $\xi_9=10$, $\xi_{11}=20$, $\xi_{14}=60$
(distribution 6 in table <ref>).
This distribution is equivalent to an elliptically anisotropic distribution with effective coefficients $\xi_2(\Omega)$ and $\xi_9(\Omega)$

as shown in figure <ref>.
The new parameters $\xi_4$, $\xi_{11}$ and $\xi_{14}$ effectively modify the constant value of $\xi_9$ and generate a non-zero $\xi_2$.

Spherical plots of the coefficients $\xi_2(\Omega)$ (left panel) and $\xi_9(\Omega)$ (right panel) that correspond to the mapping of distribution 6 in table <ref> onto an elliptically anisotropic distribution. The values in the figure are obtained with $\beta=1$ and $\mu=0$.

The parameters $\xi_6$, $\xi_8$ and $\xi_{13}$ are coefficients of terms in equation (<ref>) with $e_1$ and $e_3$ individually odd.
These terms produce a distortion of the isotropic distribution that is qualitatively different.
Some examples of the dispersion relations obtained from distributions with non-zero values of these coefficients are presented in section <ref>.

To illustrate more directly the roles of the different anisotropy parameters we look at  plots of  $C_\xi \, n\big(k H_\xi(\vec v)\big)$ with $\beta=1$ and $\mu=0$, in the $(p_x,p_y)$ and $(p_x,p_z)$ planes.
In figure <ref> we show the results  for several choices of anisotropy parameters.
The first panel shows, for reference, the isotropic distribution in the $(p_x,p_y)$ plane, the figure for the $(p_x,p_z)$ plane is identical. In the second panel we show a prolate distribution in the $(p_x,p_z)$ plane, and one sees clearly that the distribution is stretched along the beam axis.
The third panel shows the squeezing of an oblate distribution along the beam axis.
In the fourth and fifth panels we show an elliptically oblate distribution, which is asymmetric in both planes.
In the sixth panel we show an example of the kind of structure that can be obtained with a specific choice of the extra anisotropy parameters we have introduced.

Contour plots of the distribution for different choices of the anisotropy parameters. In each case the parameters that are not specified are set to zero.

§ DISPERSIONS RELATIONS

§.§ Analytic structure

The isotropic and anisotropic dressing functions have the same  analytic structure, which can be determined from equations (<ref>, <ref>), using that equations (<ref>, <ref>) give $P\cdot V = p(\hat p_0-x')$ where $\hat p_0 = p_0/p$.

When $\hat p_0$ is pure real the epsilon prescription is needed to define the integrands at $\hat p_0=x'$.
Shifting the integration variable $\vec v\to -\vec v$ one shows that $\pi_i^*(\hat p_0) = \pi_i(-\hat p_0)$, where the subscript $i$ indicates any component of the polarization tensor. This means that when $\hat p_0$ is real, the real part of each dressing function is even in $\hat p_0$ and the imaginary part is odd. We also note that all dressing functions are pure real for $\hat p_0>1$, since the imaginary part comes from the discontinuity at $\hat p_0=x'$ and $|x'|<1$.

When $\hat p_0$ is not pure real the epsilon prescription is not needed, and we have $\pi_i^*(\hat p_0) = \pi_i(\hat p^*_0)$. When $p_0$ is pure imaginary $\pi_i(\hat p^*_0) = \pi_i(\hat p_0)$, which can be shown by shifting the integration variable, and therefore the dressing functions are pure real and even in $\hat p_0$.

The information above can be summarized as
Re(p̂_0)>1  and  Im(p̂_0)=0       →       π_i(p̂_0)   real  and even

p̂_0  imaginary         →       π_i(p̂_0)   real  and even

Re(p̂_0)<1  and  Im(p̂_0)=0      →
Re( π_i(p̂_0))    even

Im (π_i(p̂_0))    odd . .		It is easy to see that all components of the polarization tensor satisfy
Π(p_0^*) = Π^*(p_0)  for  p_0  complex or pure imaginary

Re[Π(p_0)] = Re[Π(-p_0)]  and  Im[Π(p_0)] = -Im[Π(-p_0)]  for  p_0  pure real .

§.§ Method

The integrand for the polarization tensor is given by equations (<ref>, <ref>, <ref>).
The nine scalar components $\pi_1$-$\pi_9$ are obtained by contracting with the appropriate projection operators, which are easily obtained from equations (<ref>, <ref>). These calculations are straightforward but tedious, and they are therefore done with Mathematica.
The resulting integrals are calculated numerically using Gauss-Legendre quadrature.
For real valued $\hat p_0$, the real part of the dressing functions has a single integrable singularity at $x'=\hat p_0$ that is easily handled with a numerical principal part prescription, and the imaginary part can be obtained from the residue of the pole. For imaginary valued $\hat p_0$ there are no difficulties. After the dressing functions are obtained, the dispersion equation is solved with a standard binary search algorithm.

§.§ Real solutions

For any choice of the anisotropy parameters and the chiral chemical potential $\mu_5$ there are three real solutions.
We show some examples of the dispersion relations in figure <ref>.

Real solutions for distributions 1, 3 and 8 (as defined in table <ref>). In the bottom right panel the lighter lines are the soutions with $\hat\mu_5=0$ and the dark lines are $\hat\mu_5=1$. In all graphs the light gray line shows the light cone.

§.§ Imaginary solutions

We are primarily interested in imaginary solutions, because of their potential to influence plasma dynamics.
We will compare the imaginary modes produced by different choices of the anisotropy parameters, and explore their dependence on these parameters as a function of wave vector.
We have previously used the notation $\gamma(\vec p)$ to indicate the dependence of the dispersion relation on the wave vector, but in this section we will use the more explicit notation $\gamma(p,x,\phi)$ where $x=\cos(\theta)$.
For different sets of anisotropy parameters
we compare the maximum value of the imaginary solution,
the position of this maximum,
and the integral of the solution over its domain, which we write  $\int\gamma\equiv \int_0^\infty dp \int_{-1}^1 dx\int_0^{2\pi} d\phi \,\gamma(p,x,\phi)$.

We define three variables that characterise the degree of oblateness, the transverse eccentricity, and the azimuthal asymmetry of a given distribution:
χ= ⟨k_x^2⟩+ ⟨k_y^2⟩/2⟨k_z^2⟩-1

ε= ⟨k_y^2⟩-⟨k_x^2⟩/⟨k_y^2⟩+ ⟨k_x^2⟩

v_2 = ∫dk   k^2 ∫dϕ  cos(2ϕ)   n(k H_ξ(v⃗))_x=0/∫dk   k^2 ∫dϕ  n(k H_ξ(v⃗))_x=0  .
The angle brackets denote averaging over the momentum phase space, for example $\langle k_x^2\rangle = \int d^3 k k_x^2 C_\xi n(k H_\xi(\vec v))$.

In table <ref> we show some of the main features of the largest imaginary solution of the dispersion relation for several different distributions.

In all cases the anisotropy parameters that are not specified are zero.
The first distribution is the isotropic one. The second and third have spheroidal anisotropy, with parameters that make them prolate and oblate, respectively. The fourth is a distribution with spheroidal anisotropy that is strongly oblate.
The fifth is an elliptically anisotropic distribution.
Distributions 6-9 are examples constructed using the additional anisotropy parameters that we have introduced.
Distributions 6 and 7 involve only terms for which $e_1$ and $e_3$ are even (see the discussion at the beginning of section <ref>).

The second to last column shows that the maximum of the largest mode can be much greater than the result obtained from a distribution with spheroidal anisotropy with a comparable value of the oblateness parameter $\chi$.
In addition, the last column shows that the integral of the solution over its full domain is greatly enhanced, relative to the result obtained from a comparable distribution with spheroid oblateness.
For example, distribution 7 has approximately the same value of $\chi$ as distribution 3, but $\gamma_{\rm max}$ is about twice as big, and the integral of the solution is almost 9 times larger.
Distributions 8 and 9 are examples for which $e_1$ and $e_3$ are individually odd for some of the terms in (<ref>).
The magnitudes of the modes produced by these distributions are similar to what is obtained from distributions without these odd coefficients, but the angular dependence of the solutions is much more complex. At the end of this section we  show several plots to demonstrate this.

Finally, the table shows that the effect of a non-zero chiral chemical potential is always to increase the magnitude and domain of the imaginary solution. The solutions produced with even moderate anisotropy are much larger than those obtained with the same chemical potential in an isotropic system.

#          $\xi_i$         $\hat\mu_5$         $\chi$          $\varepsilon$           $v_2$            $\gamma_{\rm max}(p,x,\phi)$            $\int \gamma$
1     $\xi_0=1$     0     0     0     0     -     -

$0.3$                    $4.7\times10^{-3}$(.19,-,-)     0.017

2     $(\xi_0,\xi_9)=(11,-10)$     0      -0.91     0     0     0.16(0.45,0,-)     0.59

0.3                             0.23(0.68,0,-)     0.86

3    $(\xi_0,\xi_9)=(1,10)$     0      10     0     0     0.19(0.63,1,-)     0.08

0.3                     0.34(1.07,1,-)     0.24

4     $(\xi_0,\xi_9)=(1,500)$     0      500     0     0     0.39(1.57,1,-)     0.32

0.3                     1.60(21.0,1,-)     6.90

5      $(\xi_0,\xi_2,\xi_9)=(1,5,10)$     0     5.42     0.71     0.59     0.20(0.58,1,0)     0.17

0.3                    0.29(0.88,1,0)     0.34

6     $\xi_{0,4,9,11, 14}=(1,20,10,20,60)$      0     17.9     0.85     0.7     0.33(0.78,1,1.49)     0.56

0.3                    0.50(1.23,1,0.06)     1.15

7     $\xi_{0,4,9,11, 14}=(1,40,-2,40,55)$      0     10.1     0.91     0.85     0.32(0.68,0.79,0)     0.70

0.3                    0.43(1.06,0.80,0)     1.27

8     $\xi_{0,2,4,6,8,9,11,13,14}=$     0     8.15     0.90     0.97     0.36(0.88,0.67,$\pi$)     0.61

$(1,8,37,-21,21,8,15,-28,40)$      0.3                    0.56(1.65,0.63,$\pi$)     1.27

9     $\xi_{0,9,13,14}=(1,1,20,20)$     0     1.80     -0.35     0     0.23(0.70,0.85,0)     0.16

0.3                    0.31(0.94,0.86,0)     0.27

Different sets of anisotropy parameters, the corresponding values of the parameters $\xi$, $\varepsilon$ and $v_2$ defined in equation (<ref>), the size and location of the largest imaginary solution, and the integral of the solution over its domain. In each case the anisotropy parameters that are not specified are set to zero.

For a given distribution, the number of imaginary solutions is 0 or 1 or 2, depending on the magnitude and direction of the wave vector.
For each mode there is a critical vector $p_{\rm crit}(x,\phi,\xi_i)$ such that the imaginary mode appears for $p$ less than this critical value.
Analytic results for the critical wave vectors can be obtained in the limit of weak anisotropy using a Nyquist analysis. The method is explained in detail in Appendix <ref>, where we show how to identify the angular variables that produce the largest critical wave vectors for a given set of anisotropy parameters.

To understand the structure of the solutions obtained from a given distribution in more detail, we show several figures below.

In figure <ref> we compare the imaginary solutions obtained from the distributions 3 and 6 for a specific choice of angles. The figure illustrates the concept of the critical wave vector.
For example, in the right panel of figure <ref>, the critical wave vectors for modes 1 and 2
are, respectively, 2.33 and 0.62.
It is true in general that when the domain of the solution increases, the maximum value of the solution also increases, as seen in figure <ref>.

The imaginary solutions of the dispersion equation ${\cal D}=0$ obtained from distributions 3 (left panel) with $\hat\mu_5=0$ and $x=1$,  and the two imaginary solutions obtained from distribution 6 (right panel) with $\hat\mu_5=0$, $x=1$ and $\phi=\pi/2$.

In figure <ref> we show the largest imaginary mode obtained from distribution 6 for different choices of angles, and two different values of the chiral chemical potential.
The figure shows that increasing the chiral chemical potential increases the domain and size of the solution, but not uniformly as a function of the angular variables.

The largest imaginary solution of the dispersion equation ${\cal D}=0$ using the distribution 6, two different chiral chemical potentials, and several different values of the angular variables.

The dependence of the imaginary modes on the angular variables can be seen most clearly using contour plots.
In the top left panel of figure <ref> we show a contour plot of the magnitude of the largest imaginary solution for distribution 3 with $\hat\mu_5=0$, as a function of $p$ and $x$ (there is no $\phi$ dependence). The solution is non-zero in a small region of the phase space where $x$ is close to 1 and $p$ is small.

In the top right panel of figure <ref> we show the largest imaginary solution for distribution 2 with $\hat\mu_5=0$,
which is also $\phi$ independent. The magnitude of the solution is smaller, but the domain is larger.

The bottom two panels show the $\phi$ dependent solutions obtained from distribution 5.

Contour plots of the largest imaginary solution for several distributions with $\hat\mu_5=0$. The top left and top right panels show the $\phi$ independent solution obtained from distributions 3 and 2, respectively. The bottom panels are solutions obtained from distribution 5.

The distributions 6-9 involve anisotropy parameters that correspond to higher spherical harmonics, and have a much richer structure.
In figures <ref> and <ref> we show some contour plots of the largest imaginary solution obtained from these distributions with $\hat\mu_5=0$, for different choices of $x=\cos(\theta)$ and $\phi$.
A non-zero value of the chirality parameter $\hat\mu_5$ produces solutions that are larger and spread over a larger domain (as shown in table <ref>). To show more clearly the effect of $\hat\mu_5$ we plot in figure <ref> and <ref> the difference of the solutions obtained with $\mu_5$=0.3 and $\mu_5=0$ for the two distributions meg8 and meg4. The dispersion relation with $\hat\mu_5=0.3$ is shifted by an amount $\Delta p = 1.25-0.79$ for meg8 and $\Delta p = 1.65-0.88$ for meg4, which corresponds to the difference in the location of the maximum values (see table <ref>), so that the maxima occur at the same values of the coordinates $(p,x,\phi)$.

Contour plots of the largest imaginary solution from the distributions 6 (left two panels) and 7 (right two panels) with $\hat\mu_5=0$.

Contour plots of the largest imaginary solution from the distributions 9 (left two panels) and 8 (right two panels) with $\hat\mu_5=0$.

The effect of the anisotropy parameters which correspond to terms with $e_1$ and $e_3$ individually odd (see the discussion at the beginning of section <ref>) can be seen by calculating $s_2\equiv \int_0^\infty dp \int^1_0 dx \int^\pi_0 d\phi \,  \sin(2\phi) \gamma(p,x,\phi)$.
Any distribution for which $\xi_6=\xi_8=\xi_{13}=0$ will give $s_2=0$, which corresponds physically to symmetry across the reaction plane.
If any of these anisotropy parameters are non-zero, the value of $s_2$ will also be non-zero.
Physically this effect could be produced by a collision at non-zero impact parameter of two ions of different size.
Distribution 6 and 7, for which we should find $s_2=0$, give $s_2 \approx 9\times 10^{-10}$ and $s_2 \approx 7\times 10^{-8}$, which can be taken as a measure of the numerical error in our computation.
Distributions 8 and 9 give, respectively,  $s_2=-0.059$ and $s_2=0.068$.

Finally, to show the effect of $\hat\mu_5\ne 0$, we plot in figure <ref> the solutions obtained with $\hat\mu_5$=0.3 using the distributions 6 and 7. Comparison with figure <ref> shows that the chiral chemical potential increases the size and domain of the imaginary solution, but not uniformly.

Contour plots of the largest imaginary solution for the distributions 6 (left two panels) and 7 (right two panels) with $\hat\mu_5=0.3$.

§ CONCLUSIONS

We have studied plasmons in anisotropic plasmas using a hard loop approach within the quasi-particle regime.
We have used more general distribution functions than have been previously considered, with a technique that could be easily generalized
to account for the anisotropies that are relevant in any given physical situation.

The distribution function we have used depends on 9 anisotropy parameters and is normalized so that at zero chemical potential the Debye mass parameter is the same for all distributions.
In an ultra-relativistic plasma at zero chemical potential this parameter is the only dimensionful scale. It is therefore meaningful to compare the solutions obtained from the dispersion equations for different choices of anisotropy parameters.

Imaginary solutions to the dispersion equation are important because they are associated with plasma instabilities, which have the potential to strongly influence plasma dynamics.
The important effect of imaginary collective modes on physical quantities in anisotropic plasmas, like transport coefficients, production rates, and bound state properties, is well documented in the literature.

We have shown that the more general distribution functions we have introduced can significantly increase imaginary solutions of the dispersion equation in terms of both the magnitude of the mode, and the domain of wave vectors over which it exists.
The distribution that is most commonly used in the literature depends on only one anisotropy parameter, which describes the oblateness of a distribution that is squeezed in one direction.
An example is our distribution 3.
In comparison,  our distribution 7, which has approximately the same value of the parameter $\xi$ that characterises its oblateness (see table <ref>), produces a largest imaginary mode that is almost twice as big, and the integral of the solution over the momentum phase space is approximately 10 times larger.

We have also shown that the effect of the chiral chemical potential is greatly enhanced by anisotropy. This is seen clearly from the results in table <ref> which show that, relative to the isotropic result, the imaginary modes are much greater when
even very moderate anisotropy is introduced.

This work was supported by the Natural Sciences and Engineering Research Council of Canada Discovery  Grant  program.

§ RETARDED POLARIZATION TENSOR

In this Appendix we calculate the one-loop retarded polarization tensor at finite temperature and chemical potential with the real-time formulation of finite temperature field theory, using the Keldysh representation [51, 52, 53].

The electron propagator is a 2$\times$2 matrix of the form
G     =     (
G_rr    G_ra

G_ar     G_aa

= (
G_sym    G_ret

G_adv     0

where the retarded, advanced and symmetric propagators are given by
G_ret(p_0,p)  = (P + m)r(p_0,p)

G_adv(p_0,p) = (P + m) a(p_0,p)

G_sym(p_0,p) = (P + m) f(p_0,p)
r(p_0,p)     ≡    1/P^2-m^2+i0^+ sgn(p_0)

a(p_0,p)    ≡    1/P^2-m^2-i0^+ sgn(p_0)

[2mm] f(p_0,p)    ≡    -2πi
[ (1-2 n^+(p_0) ) Θ(p_0) +  (1-2 n^-(p_0) ) Θ(-p_0) ]
δ(P^2 - m^2)
n^+(p_0) = 1/e^β(p_0-μ)+1    and      n^-(p_0) = 1/e^-β(p_0-μ)+1 .

The self-energy has the form
Π   =    (
Π_rr     Π_ra

Π_ar     Π_aa
= (
0     Π_adv

Π_ret     Π_sym
and the vertex function is a 2$\times$2$\times$2 tensor which can be written
Γ^μ   =    -ie γ^μ(

{Γ_rrr,Γ_rra}     {Γ_rar,Γ_raa}

{Γ_arr,Γ_ara}     {Γ_aar,Γ_aaa}

) = -ie γ^μ(

{0,1}     {1,0}

{1,0}     {0,1}

) .

The contribution to the retarded self-energy from the one-loop diagram is
i Π_ret^μν (p_0,p) = i Π_ar^μν(p_0,p) = (-ie)^2/2∑_ii'jj'∫d^4k/(2π)^4
γ^μΓ_aij   G_jj'(k) G_i'i(k-p)   γ^νΓ_ri'j' .
The sum over Keldysh indices $\{i,i',j,j'\}\in\{r,a\}$ is easily done because $G_{aa}=0$ and a vertex function with an odd number of $a$ indices vanishes.  The result is

Π_ret^μν (p_0,p)
= ie^2/2
Tr[γ^μ(K + m) γ^ν(K-P + m)]
[r(K) f(K-P) + f(K) a(K-P)]

where we have used the short hand notation $r(K)$ to represent $r(k_0,k)$, and similarily for other propagator components. From this point on we will suppress the subscript that indicates the retarded component of the polarization tensor. In addition we consider only the relativistic limit $T\gg m$ and therefore we drop the electron mass in the integrand in equation (<ref>). We also introduce the short-hand notation $dK\equiv d^4k/(2\pi)^4$.

We will consider a chirally asymmetric plasma which is characterised by a difference in the chemical potentials of the right and left handed fermions, which are denoted $\mu_R$ and $\mu_L$. We define

Π_R^μν (p_0,p)
= ie^2/2
Tr[(1+γ^5)γ^μK  γ^ν(K-P)]
[r(K) f(K-P) + f(K) a(K-P)]|_μ=μ_R

Π_L^μν (p_0,p)
= ie^2/2
Tr[(1-γ^5)γ^μK γ^ν(K-P)]
[r(K) f(K-P) + f(K) a(K-P)] |_μ=μ_L

Π^μν(p_0,p) = 1/2(Π_R^μν (p_0,p) +Π_L^μν (p_0,p))   .
The integrals in equations (<ref>) can be rewritten by performing a shift of the 4-vector $K \to K+P$ on the term in square brackets that does not have a factor $f(K)$. All terms in the resulting expression have a factor $f(K)$ and, using equation (<ref>), this factor is proportional to $\delta(K^2)$ in the relativistic limit. The integral over $k_0$ can be done using this delta function.
We work in time-like axial gauge and calculate only the spatial components of the polarization tensor. After doing the $k_0$ integral we obtain the results in equations (<ref>, <ref>).

§ NYQUIST ANALYSIS

The dispersion equation that we have solved numerically is ${\cal D}=0$ where ${\cal D}$ is given in equation (<ref>).
When numerical methods are used to find the solutions of an equation, one must input the range over which the search will be performed, and solutions outside this range will be missed.
A Nyquist analysis can be used to determine the number of solutions that exist so that the search range can be extended as needed.

To explain the idea we discuss a generic equation of the form

f(p_0) = 0
and  define the function

F(p_0) ≡f^'(p_0)/f(p_0) = d/dp_0lnf(p_0)  .
We consider the contour integral

∮_C dp_0/2πi F(p_0)
where the contour is a positively (counter-clockwise) oriented closed loop, which is chosen so that $F(p_0)$ is analytic inside the loop except at isolated points. The integral is equal to the sum of the residues. The residue of $F(p_0)$ at a zero of $f(p_0)$ of order $l$ is $l$, and the residue of $F(p_0)$ at a pole of $f(p_0)$ of order $l$ is $-l$. This gives

∮_C dp_0/2πi F(p_0) = n_Z - n_P
where $n_Z$ and $n_P$ are the numbers of zeros and poles of $f(p_0)$ inside the contour $C$, taking into account the fact that each zero and pole of order $l$ is counted $l$ times. The purpose of the Nyquist analysis is to determine $n_Z$.

The anisotropic dressing functions are obtained numerically from the integrals in equations (<ref>, <ref>).
Using $P\cdot V = p(\hat p_0 - x')$ (see equations (<ref>, <ref>)) one finds that the dressing functions have a cut along  the real axis for $|\hat p_0|<1$, and the dispersion equation therefore has a cut for $|p_0|<p$.
This means that the contour $C$ can be chosen as in Fig. <ref>.

The contour in the plane of complex $\hat p_0$ which is used to compute the number of solutions of the dispersion equations.

The integrals along the lines connecting the circular contour $C_\infty$ to $C_{\rm cut}$ always compensate each other and therefore the contour integral (<ref>) equals

∮_C_∞ dp_0/2πi F(p_0) +
∮_C_cut dp_0/2πi F(p_0) = n_Z - n_P .

The contribution from the big circle is calculated by writing $p_0 = |p_0|e^{i \phi}$ and taking $|p_0|\to\infty$ which gives

∮_C_∞ dp_0/2πi F(p_0) = lim_|p_0| →∞ p_0 F(p_0) ≡n_∞ .

The integral along the cut can be calculated using that $F(p_0)$ is the logarithmic derivative of $f(p_0)$ which gives

∮_C_cut dp_0/2πi F(p_0)
= 1/2πi∮_C_cut  d/dp_0lnf(p_0)
= 1/2πi ( ln f(p_0end) - ln f(p_0start) ) ≡n_W
where $p_{0{\rm start}}$ is the (arbitrarily chosen) starting point of the contour which encloses the cut, and $p_{0{\rm end}}$ is the end point.
The quantity $n_W$ is the number of times that the curve in the plane of complex $f$
travels counter-clockwise around the point $f=0$, and is called a winding number.
Combining the results (<ref>, <ref>), we have that equation (<ref>) gives

n_Z = n_P + n_∞+ n_W  .

A Nyquist analysis of the equation ${\cal D}=0$ can be done analytically in the weakly anisotropic limit.
When $\mu_5=0$ the dispersion equation in the weakly anisotropic limit factors into three separate equations
f^(1)(p_0,p⃗) =  P^2-π_1(p_0,p⃗)=0

f^(2)(p_0,p⃗) = 1/p_0^2(p_0^2-π_2(p_0,p⃗)) = 1-π_2(p_0,p⃗)/p_0^2 =0

f^(3)(p_0,p⃗) = P^2-π_3(p_0,p⃗)=0  .
When $\mu_5\ne 0$ equations (<ref>, <ref>) are replaced by the weakly anisotropic versions of equations (<ref>, <ref>) which are
f^(4)(p_0,p⃗) = P^2-(1/2(π_1(p_0,p⃗)+π_3(p_0,p⃗)) + π_8(p_0,p⃗))=0

f^(5)(p_0,p⃗) = P^2-(1/2(π_1(p_0,p⃗)+π_3(p_0,p⃗)) -π_8(p_0,p⃗))=0  .
In all cases it is easy to see that $n_P=0$.
In equation (<ref>) we have divided out a zero mode, so that $n^{(2)}_\infty = 0$, and for all of the other dispersion equations it is straightforward to see that $n_\infty = 2$.

The non-trivial part of the calculation is the determination of the winding number, which is found by mapping the closed contour $C_{\rm cut}$ in the plane of complex $p_0$ onto a path in the plane of complex $f(p_0)$.
The points $p_{0{\rm start}}$ and $p_{0{\rm end}}$ in equation (<ref>) have the same modulus, but their phases differ by $2\pi$.
The winding number equals the number of times that the curve in the plane of complex $f$
travels counter-clockwise around the point $f=0$ (some detailed examples of these mappings can be found in Ref. [9]).
There are two possibilities, which are shown in figure <ref>.

Mappings of the contour $C$ in the plane of complex $p_0$. In the right panel the corresponding winding number is two, and in the left panel it is zero.

In the right panel, the mapping circles the origin twice, and the winding number is 2, while in the left panel the winding number is zero.
To distinguish between the two possibilities for each dispersion equation, we need to determine the sign of Re[$f(0,p)$].

We discuss first the dispersion equation (<ref>), which is different from the other four.
In the weakly anisotropic limit we have
f^(2)(0,p⃗) = 2+f̃(ξ_i,x,ϕ)
where $\tilde f(\xi_i,x,\phi)$ is a dimensionless function of the angular variables that depends linearly on the anisotropy parameters.
Since we have assumed weak anisotropy we have $f^{(2)}(0,\vec p)$ always positive, which means $n_W=2$, and therefore equation (<ref>) gives  $n_Z=2$.
These two solutions are always present, and numerically we find they are always real. They are anisotropic modifications of the HTL longitudinal modes (denoted $\pm\omega_L(\omega)$).

For all of the remaining dispersion equations the number of solutions  is
n_Z=2+n_W .
For equations (<ref>) and (<ref>) the winding number is determined from the sign of the functions
f^(1)(0,p⃗) = -(p^2- p^2_1crit(ξ_i,x,ϕ))

f^(3)(0,p⃗) = -(p^2- p^2_3crit(ξ_i,x,ϕ)) .
When the chirality parameter $\mu_5$ is non-zero we consider equations (<ref>, <ref>) instead of equations (<ref>, <ref>).
For weak anisotropy the winding number is determined from the sign of the functions
f^(4)(0,p⃗) = -(p +  p_5crit(ξ_i,x,ϕ))

f^(5)(0,p⃗) = -(p - p_5crit(ξ_i,x,ϕ)) .

The results for $p^2_{1{\rm crit}}$, $p^2_{3{\rm crit}}$ and $p_{5{\rm hat}}$ are given in equations (<ref>, <ref>, <ref>).
15 p_1crit^2    =
ξ_8 x √(1-x^2) (4 x^2-1) cos^3(ϕ)
+(ξ_11 (1-2 x^2)^2-5 ξ_2 (x^2-1)) cos^2(ϕ)

+    1/2 x √(1-x^2) cos(ϕ) (-9 ξ_8 cos(2 ϕ)+10 ξ_6+9 ξ_8+6 ξ_13-8 ξ_13 x^2)

-    sin^2(ϕ) (5 ξ_2-3 ξ_4+ξ_11+9 ξ_4 (x^2-1) cos(2 ϕ)+9 ξ_4 x^2-3 ξ_11 x^2)

+    x^2 (5 ξ_9+6 ξ_14-4 ξ_14 x^2)+2 ξ_4 (-2 x^4+x^2+1) cos^4(ϕ)

15 p_3crit^2    =    -3 ξ_4 x^2 sin^2(2 ϕ)+ξ_11 (2 x^2-1) sin^2(ϕ)
+4 ξ_8 x √(1-x^2) (4 x^2-1) cos^3(ϕ)

+    x √(1-x^2) cos(ϕ) (-3 ξ_8 cos(2 ϕ) + 10 ξ_6+3 ξ_8+4 ξ_13 (3-4 x^2))
+5 ξ_9 (2 x^2-1)

+    cos^2(ϕ) (6 ξ_4 sin^2(ϕ)+5 ξ_2 (1-2 x^2)+2 ξ_11
(8 x^4-8 x^2+1))

+    2 ξ_4 (-8 x^4+4 x^2+1) cos^4(ϕ) - 2 ξ_14 (8 x^4-12 x^2+3)

30/μ̂_5 p_5hat    =    30+10 ξ_2+6 ξ_4-5 ξ_9-9 ξ_14+ 10 √(1-x^2) x (3 (ξ_6+ξ_8 x^2)
+ξ_13 (3-4 x^2)) cos(ϕ)

-    5 (x^2-1)
(ξ_4 (x^2-1) cos(4 ϕ)
-2 ξ_8 x √(1-x^2) cos(3 ϕ)

+(3 ξ_2+2 ξ_4+ξ_11+4 ξ_4 x^2-4 ξ_11 x^2) cos(2 ϕ))

-5 x^2 (3 (ξ_2-2 ξ_9+ξ_4 x^2)
+4 ξ_14 (2x^2-3))+ξ_11 (20 x^4-15 x^2+2)  .

In the isotropic case $p^2_{1{\rm crit}}(\xi_i,x,\phi) = p^2_{3{\rm crit}}(\xi_i,x,\phi) = 0$, the right side of equations (<ref>, <ref>) cannot be positive, and therefore the winding number is zero and both dispersion equations have two solutions. Remembering that in the isotropic limit $\pi_1=\pi_3=\Pi_T$, we see that these are the HTL transverse solutions ($\pm\omega_T(p)$).
When non-zero anisotropy parameters are introduced, and when $p$ is small enough, there can be choices of the angles and anisotropy parameters for which $p^2_{1{\rm crit}}(\xi_i,x,\phi)>p^2>0$ and/or $p^2_{3{\rm crit}}(\xi_i,x,\phi)>p^2>0$.
In this case, the winding number is two, which means that extra solutions to equations (<ref>, <ref>) exist.

From equation (<ref>) we have that in the isotropic limit $p_{5{\rm crit}}(\xi_i,x,\phi)=\hat \mu_5$, which means that $n_W^{(4)}=0$ and $n_W^{(5)}=2$ if $p<\hat\mu_5$ and zero otherwise. Thus we have that for $\mu_5\ne 0$ an additional pair of solutions, which turn out to be imaginary, exist even in the isotropic limit (see figure <ref>). The domain of these modes can be larger for certain orientations of the wave vector if the system is anisotropic.

The expressions in equations (<ref>, <ref>, <ref>) are only valid in the weakly anisotropic limit, however, we can use them to obtain information about the orientations of the wave vector $\vec p$ for which large imaginary solutions could be found with specific sets of anisotropy parameters. To obtain more insight, we can plot the coefficient of each anisotropy parameter in equations (<ref> - <ref>). We show the results for three of the coefficients in figure <ref>.

The first panel shows the dependence of the critical wave vectors on $\xi_9$, which is the coefficient of $(n_3\cdot v)^2$ in equation (<ref>). Using our coordinate system this dot product is independent of the azimuthal angle $\phi$,
and therefore we have  that the terms in the critical wave vectors that are proportional to $\xi_9$ are also independent of $\phi$.
The second panel shows the dependence of the critical wave vectors on $\xi_2$ which is the coefficient of $(n_3\cdot v)^2$. Both  the $\xi_9$ and $\xi_2$ terms have been included in previous calculations. In the third panel we show the dependence of the critical wave vectors on $\xi_6$, which is the coefficient of $(n_1\cdot v)(n_3\cdot v)$, and in the last panel we show the dependence on
$\xi_8$, the coefficient of $(n_1\cdot v)^3(n_3\cdot v)$.

Contour plots of the coefficients of the anisotropy parameters for $p_{1{\rm crit}}^2$ (left), $p_{3{\rm crit}}^2$ (center) and $p_{5{\rm crit}}$ (right).  The axes show $x=\cos(\theta)\in(0,1)$ and $\phi\in(0,\pi/2)$.

[1]
P. Romatschke and M. Strickland,

Phys. Rev.  D 68, 036004 (2003).

[2]
P. Romatschke and M. Strickland,

Phys. Rev.  D 70, 116006 (2004).

[3]
St. Mrówczyński,

Phys. Lett. B 314, 118 (1993).

[4]
P. B. Arnold, J. Lenaghan and G. D. Moore,

JHEP 0308, 002 (2003).

[5]
A. Kurkela and G. D. Moore,

JHEP 1112, 044 (2011).

[6]
A. Kurkela and G. D. Moore,

JHEP 1111, 120 (2011).

[7]
A. Ipp, A. Rebhan and M. Strickland,

Phys. Rev. D 84, 056003 (2011).

[8]
M. Attems, A. Rebhan and M. Strickland,

Phys. Rev. D 87, 025010 (2013).

[9]
M. E. Carrington, K. Deja and St. Mrówczyński,

Phys. Rev. C 90, 034913 (2014).

[10]
B. Schenke and M. Strickland,

Phys. Rev. D 76, 025023 (2007).

[11]
B. S. Kasmaei, M. Nopoush and M. Strickland,

Phys. Rev. D 94, 125001 (2016).

[12]
B. S. Kasmaei and M. Strickland,

Phys. Rev. D 97, 054022 (2018).

[13]
R. Baier and Y. Mehtar-Tani,

Phys. Rev. C 78, 064906 (2008).

[14]
P. Romatschke and M. Strickland,

Phys. Rev. D 71, 125008 (2005).

[15]
M. E. Carrington, K. Deja and St. Mrówczyński,

Phys. Rev. C 92, 044914 (2015).

[16]
M. E. Carrington, St. Mrówczyński and B. Schenke,

Phys. Rev. C 95, 024906 (2017).

[17]
A. Dumitru, Y. Guo and M. Strickland,

Phys. Lett. B 662, 37 (2008).

[18]
Y. Burnier, M. Laine and M. Vepsalainen,

Phys. Lett. B 678, 86 (2009).

[19]
A. Dumitru, Y. Guo and M. Strickland,

Phys. Rev. D 79, 114003 (2009).

[20]
K. Boguslavski, B. S. Kasmaei and M. Strickland,

arXiv:2102.12587 [hep-ph].

[21]
M. Martinez and M. Strickland,

Phys. Rev. C 78, 034917 (2008).

[22]
S. Hauksson, S. Jeon and C. Gale,

Phys. Rev. C 103, 064904 (2021).

[23]
B. S. Kasmaei and M. Strickland,

Phys. Rev. D 102, 014037 (2020).

[24]
B. S. Kasmaei and M. Strickland,

Phys. Rev. D 99, 034015 (2019).

[25]
Y. Akamatsu and N. Yamamoto,

Phys. Rev. Lett. 111, 052002 (2013).

[26]
A. Kumar, J. R. Bhatt and P. K. Kaw,

arXiv:1405.2865 [hep-th].

[27]
K. Fukushima, D. E. Kharzeev and H. J. Warringa,

Phys. Rev. D 78, 074033 (2008).

[28]
D. Kharzeev, R. D. Pisarski and M. H. G. Tytgat,

Phys. Rev. Lett. 81, 512 (1998).

[30]
D. Kharzeev, A. Krasnitz and R. Venugopalan,

Phys. Lett. B 545, 298 (2002).

[31]
D. Kharzeev,

Phys. Lett. B 633, 260 (2006).

[32]
D. E. Kharzeev, L. D. McLerran and H. J. Warringa,

Nucl. Phys. A 803, 227 (2008).

[33]
D. E. Kharzeev,

Annals Phys. 325, 205 (2010).

[34]
A.N. Redlich and L.C.R. Wijewardhana, Phys. Rev. Lett. 54, 907 (1985).

[35]
M. Laine,

JHEP 10, 056 (2005).

[36]
J. P. Blaizot and E. Iancu,

Phys. Rev. Lett. 70, 3376 (1993).

[37]
P. F. Kelly, Q. Liu, C. Lucchesi and C. Manuel,

Phys. Rev. Lett. 72, 3461 (1994).

[38]
D. F. Litim and C. Manuel,

Phys. Rept.  364, 451 (2002).

[39]
St. Mrówczyński and M. H. Thoma,

Phys. Rev. D 62, 036011 (2000).

[40]
St. Mrówczyński, A. Rebhan and M. Strickland,

Phys. Rev. D 70, 025004 (2004).

[41]
D. T. Son and N. Yamamoto,

Phys. Rev. Lett. 109, 181602 (2012).

[42]
M. A. Stephanov and Y. Yin,

Phys. Rev. Lett. 109, 162001 (2012).

[43]
D. T. Son and N. Yamamoto,

Phys. Rev. D 87, 085016 (2013).

[44]
C. Manuel and J. M. Torres-Rincon,

Phys. Rev. D 89,  096002 (2014).

[45]
S. Carignano, C. Manuel and J. M. Torres-Rincon,

Phys. Rev. D 98, 076005 (2018).

[46]
C. Manuel and J. M. Torres-Rincon,

Phys. Rev. D 92,  074018 (2015).

[47]
S. Carignano and C. Manuel,

Phys. Rev. D 99,  096022 (2019).

[48]
S. Carignano and C. Manuel,

Phys. Rev. D 103,  116002 (2021).

[49]
M.  Le Bellac,
Thermal Field Theory
(Cambridge University Press, Cambridge, 2000).

[50]
R. Ghosh, B. Karmakar and A. Mukherjee,

Phys. Rev. D 102,  114002 (2020).

[51]
R. L. Kobes, G. W. Semenoff and N. Weiss,

Z. Phys. C 29, 371 (1985).

[52]
M. Le Bellac and H. Mabilat,

Z. Phys. C 75, 137 (1997).

[53]
M. E. Carrington, T. Fugleberg, D. S. Irvine and D. Pickering,

Eur. Phys. J. C 50, 711 (2007).


# Machine Learning-based Analysis of Electronic Properties as Predictors of
Anticholinesterase Activity in Chalcone Derivatives111A footnote for the title

Thiago Buzelli  [  Bruno Ipaves<EMAIL_ADDRESS>[  Wanda Pereira
Almeida  [  Douglas Soares Galvao<EMAIL_ADDRESS>[  Pedro Alves da
Silva Autreto<EMAIL_ADDRESS>[

###### Abstract

In this study, we investigated the correlation between the electronic
properties of anticholinesterase compounds and their biological activity.
While the methodology of such correlation is well-established and has been
effectively utilized in previous studies, we employed a more sophisticated
approach: machine learning. Initially, we focused on a set of $22$ molecules
sharing a common chalcone skeleton and categorized them into two groups based
on their IC50 indices: active and inactive. Utilizing the open-source software
Orca, we conducted calculations to determine the geometries and electronic
structures of these molecules. Over a hundred parameters were collected from
these calculations, serving as the foundation for the features used in machine
learning. These parameters included the Mulliken and Lowdin electronic
populations of each atom within the skeleton, molecular orbital energies, and
Mayer’s free valences. Through our analysis, we developed numerous models and
identified several successful candidates for effectively distinguishing
between the two groups. Notably, the most informative descriptor for this
separation relied solely on electronic populations and orbital energies. By
understanding which computationally calculated properties are most relevant to
specific biological activities, we can significantly enhance the efficiency of
drug development processes, saving both time and resources.

###### keywords:

American Chemical Society, LaTeX

UFABC] Center for Natural and Human Sciences (CCNH) of Federal University of
ABC (UFABC) UFABC] Center for Natural and Human Sciences (CCNH) of Federal
University of ABC (UFABC) UNICAMP] Institute of Chemistry, State University of
Campinas (UNICAMP) UNICAMP] Applied Physics Department and Center for
Computational Engineering and Sciences, State University of Campinas (UNICAMP)
UFABC] Center for Natural and Human Sciences (CCNH) of Federal University of
ABC (UFABC) IR,NMR,UV

##  1 Introduction

Developing novel pharmaceutical compounds is a complex, resource-intensive,
and time-consuming process. While drugs promise to improve life expectancy and
deliver effective disease treatment, they also pose the inherent risk of
adverse reactions and use abuse. The pipeline development of new
pharmacological products generally involves multi-sequential stages, such as:
1. identification and discovery of compounds exhibiting therapeutic activity;
2. rigorous in vitro testing to assess their biological properties; 3.
comprehensive in vivo studies to investigate their dynamics in animal models;
4. human clinical trials 1.

Optimizing the selection of new active compounds through computer simulations
offers significant appeal for cost reduction, reagent optimization, and
accelerated product development. Computational approaches, such as docking,
allow the analyses of molecule-ligand interactions to identify suitable
candidates that fit (energetically and geometrically) into protein pocket
binding sites 2. In principle, selecting structures with a good potential for
a desired pharmaceutical activity through interaction energy values is
effectively possible. Other approaches aim to correlate geometric features,
chemical composition, and electronic structure data with biological activity.

The use of computer-generated models has already achieved remarkable success
in correlating biological activity with molecular electronic features obtained
from electronic structure calculations. For instance, the electronic structure
of acids can be directly linked to the biological activity of their
derivatives 3. Additionally, there is a correlation between spectroscopy data
and the biological activity of steroids, highlighting the significance of the
molecules’ skeletal composition 4. Other studies established a correlation
between polycyclic aromatic hydrocarbons (PAHs) electronic structure and their
carcinogenic potential through the so-called universal quantum-molecular
descriptors 5. Recently, Autreto and co-authors successfully established a
correlation between the electronic indices and biological properties of novel
febrifugine derivatives 6. Finally, the use of the electronic theory of cancer
(ETC) revealed significant correlations between the PAHs’ carcinogenicity and
their energy difference values between the highest occupied molecular orbital
(HOMO) and the second-highest occupied molecular orbital (HOMO-1), as well as
with localized charge density within some molecular regions 7, 8.

Recent advances in machine learning (ML) techniques have been successfully
used for complex problems characterized by large volumes of dynamic,
unstructured, and error-containing data. With the exponential growth of
generating data and rapid technological advancements in software and computer
hardware, we entered the era of big data science 9. Data science combines
mathematical principles, computer science, statistical analyses, and
programming techniques to address multifaceted problems. By exploiting the
computational power of modern computers and advanced complex algorithms, ML
can create a large number of models from a single data set. This capability
allows us to compare and combine models, thereby enhancing the quality of the
analyses in extracting and identifying correlations even with sparse data 9.
It is worth emphasizing that the use of ML extends beyond problems involving
large data sets. In fact, it can be instrumental in elucidating relationships
or even unveiling underlying laws governing complex phenomena where scientific
knowledge remains undiscovered or is too intricate to comprehend through
traditional approaches 10.

In this study, we carried out an investigation on a set of $22$ molecules
categorized as anticholinesterases. These compounds exhibit the ability to
inhibit different forms of cholinesterases and have found applications in
diverse domains such as pesticide development, glaucoma treatment, and
myasthenia management 11. Notably, anticholinesterases have received
considerable attention in the context of Alzheimer’s therapy, with three out
of the four drugs currently used for its treatment belonging to the class of
anticholinesterase agents (namely, galantamine, donepezil, and rivastigmine)
11.

In the present work, We have developed a new approach to identifying active
and inactive molecules by integrating ML techniques with electronic structure
parameters. We have obtained significantly enhanced efficiency in determining
the prediction limits of our models through an extensive analysis of hundreds
of thousands of combinations simultaneously. Based on that, we developed an
optimized computational strategy to identify active and/or inactive molecules
with specific biological activity.

##  2 Methodology

###  2.1 Target molecules

The molecules under investigation exhibit a shared structural framework, as
can be seen from Figures 1 and 2. This characteristic is important as it makes
it easier to obtain structural correlations across the group set. The
biological activity of these molecules is presented based on the IC50 index of
each molecule in Table 1.

To determine the IC50 values, the compounds were diluted in phosphate buffer
(pH 7.4) and DMSO at six different concentrations: 200, 150, 100, 50, 25, and
10 micromolar. These diluted solutions were then tested against
acetylcholinesterase. The acetylcholinesterase enzyme used in the study was
obtained from Sigma-Aldrich Co and was sourced from Electrophorus electricus
(electric eel), Type VI S. The assay was conducted using Ellman’s colorimetric
method. It should be noted that the maximum concentration of DMSO used in the
experiments was 1% 12, 11.

Figure 1: Skeletal structure (chalcone group) of our anticholinesterase group set. Each atom within the skeleton represents a point where multiple electronic parameters are attributed. Table 1:  Categorization of molecules into two distinct groups (active and inactive) according to their IC50 values. ID |  IC50 | Group
---|---|---
1 | 0.039 | Active
2 | 1.471 | Inactive
3 | 1.480 | Inactive
4 | 0.261 | Active
5 | 2.567 | Inactive
6 | 4.813 | Inactive
7 | 0.172 | Inactive
8 | 0.230 | Inactive
9 | 0.170 | Inactive
10 | 0.050 | Active
11 | 0.020 | Active
12 | 0.031 | Active
13 | 0.018 | Active
14 | 0.528 | Inactive
15 | 0.028 | Active
16 | 0.739 | Inactive
17 | 1.023 | Inactive
18 | 2.528 | Inactive
19 | 0.008 | Active
20 | 0.026 | Active
21 | 0.035 | Active
22 | 0.027 | Active

Based on the IC50 values, the group of molecules was divided into two
subgroups: active and inactive, as indicated in Figure 2. The IC50 index
represents the potency of the inhibitor compound in inhibiting 50% of a
biological or biochemical process. A low IC50 value indicates a more efficient
and powerful inhibitory activity of the molecule. We considered molecules with
an IC50 value below 0.1 µM (0.1 x $10^{3}$ mol/m3) as belonging to the active
subgroup, while those with IC50 values above this threshold belonging to the
inactive subgroup. These are typical classification values found in the
literature.

Figure 2: Molecules labeled as active are represented by a red circle on their
right side, while molecules labeled as inactive are represented by a blue one
on their right side.

###  2.2 Calculation of electronic indices

After classifying the molecules into these two distinct subgroups, the
subsequent procedure involved the computation of electronic indices associated
with each element of the molecular skeleton. These indices will serve as our
(classification active/inactive) feature set and form the foundation for
creating correlation models using machine learning algorithms.

The electronic structure of the investigated molecules was obtained using the
Orca computational code 13. We used a semi-empirical approach based on the
Hartree-Fock method 14, 15, 16. The molecular geometries were optimized using
the Quasi-Newton BFGS algorithm, the PM3 semi-empirical method, and a def2-SVP
basis set 17. The SCF convergence criteria in ORCA for geometry optimizations
were set to TightSCF, with an energy change of $1.0\times 10^{-8}$ au,
$\text{TolE}=5\times 10^{-6}$, $\text{TolRMSG}=1\times 10^{-4}$,
$\text{TolMaxG}=3\times 10^{-4}$, $\text{TolRMSD}=2\times 10^{-3}$, and
$\text{TolMaxD}=4\times 10^{-3}$ 13.

###  2.3 Electronic indices - features

Using the optimized geometries of the anticholinesterase molecules (Figure 2),
we calculated the following PM3 electronic indices for each skeleton atom:

  * •

Mulliken and Lowdin orbital charges (s and p).

  * •

Mulliken and Lowdin atomic charges.

  * •

Orbital Energies values: HOMO-3, HOMO-2, HOMO-1, HOMO, GAP, LUMO, LUMO+1,
LUMO+2, LUMO+3.

  * •

Mayer valency.

HOMO, LUMO, and GAP refer to Highest Occupied Molecular Orbital, Lowest
Unoccupied Molecular Orbital, LUMO-HOMO energy value, respectively.

To identify the best descriptors, we used an algorithm called SISSO (Sure
Independence Screening and Sparsifying Operator), which combines two
techniques: SIS (Sure Independence Screening) and OS (Sparsifying Operator).
The SIS technique is related to reducing the dimensionality of the system,
while the OS is the sparsifying operator used in the compressed sensing
approach 18.

The SISSO input file consists of a worksheet that contains a molecule (labeled
example) on each line, collectively forming the data set. Each column
associated with a labeled example represents a feature, and the collection of
these columns composes the feature set. The SISSO output provides a ranking of
maps indicating the most effective separation between the active and inactive
data sets. This separation metric is determined by a pair of features known as
descriptors. These descriptors assign a point to each labeled example,
resulting in the formation of a two-dimensional graph that depends on the
chosen descriptors.

The set of input data, known as the training set, serves as the basis for the
learning algorithm. Prior to implementing machine learning methods, a subset
of the data set, molecules in our case, is typically reserved for testing the
model. This subset is referred to as the test set. In our study, we have
separated molecules $8$ and $9$ from the inactive group and molecules $10$ and
$11$ from the active group, which constitute our test set. The remaining
molecules form our training set, which is used in the machine learning
algorithm to establish a correlation model between electronic indices and
biological activity.

##  3 Results and discussions

We have performed statistical analyses based on the Pearson correlation
coefficient to select the most relevant features (Figure 7 in supplementary
info). We executed a filtering process by excluding parameters with an
absolute value of the Pearson coefficient ($|\rho|$) greater than or equal to
$0.8$. This process resulted in a refined feature set consisting of $35$
attributes. We refer to this curated collection as Pearson’s feature set. Here
are the specific parameters that compose this set:

  * •

Orbital energies: HOMO-3, HOMO-1, HOMO, LUMO+1

  * •

Mulliken population analysis: (M0, M1, M2, M3, M4, M5, M6, M8, M10, M11, M12,
M13, M15). MX indicates the Mulliken charge on atom X; (MS7, MS9, MS10). MYX
Indicates the Mulliken charge in the Y orbital (which can be S or P) of the X
atom.

  * •

Lowdin population analysis: (L2, L4, L6, L7, L8, L10). LX indicates the Lowdin
charge on the X atom; (LS0, LS1, LS5, LS7, LS11, LS13, LS14, LS15). LYX
indicates the Lowdin charge on the Y orbital of the X atom.

  * •

Mayer’s free valence analysis: All attributes pertaining to this analysis were
removed.

For each feature set, we generated a large number of models. However, on
average, only 20 of them exhibited a correlation of 100% accuracy with the
biological activity of the molecules, as determined by the electronic indices.
This correlation allowed us to effectively classify the molecules into two
distinct groups: active and inactive. Consequently, we successfully mitigated
underfitting problems, although the limited data size still posed a risk of
overfitting.

To assess the performance of these models, we carried out a validation
procedure using an independent test set. As a result, we narrowed down the
selection to four models for the Pearson feature set. These models, except for
orbital charges, are consistent with the discussions previously presented in
the literature using other approaches 19. The following parameter models
achieved a perfect 100% accuracy in separating active and inactive molecules:

  1. 1.

$\left[\frac{\textnormal{M}10/\textnormal{M}4}{\textnormal{M}1-\textnormal{M}3}\right],\left[\frac{\textnormal{M}10/\textnormal{L}4}{\textnormal{M}1-\textnormal{M}3}\right]$

  2. 2.

$\left[\frac{\textnormal{HOMO}_{-1}-\textnormal{HOMO}}{\textnormal{M}1-\textnormal{M}3}\right],\left[\frac{\textnormal{L}10/\textnormal{L}2}{\textnormal{M}1-\textnormal{M}3}\right]$

  3. 3.

$\left[\frac{\textnormal{HOMO}_{-1}-\textnormal{HOMO}}{\textnormal{M}1-\textnormal{M}3}\right],\left[\frac{\textnormal{M}10/\textnormal{M}4}{\textnormal{M}1-\textnormal{M}3}\right]$

  4. 4.

$\left[\frac{\textnormal{HOMO}_{-1}-\textnormal{HOMO}}{\textnormal{M}1-\textnormal{M}3}\right],\left[\frac{\textnormal{M}10/\textnormal{L}4}{\textnormal{M}1-\textnormal{M}3}\right]$

It is worth noting that among the four models, there is a recurring presence
of factors such as (M1-M3) and ($\Delta H$ = HOMO-1 \- HOMO), indicating their
significance in the separation of the molecular groups. It is important to
emphasize that the term $\Delta H$ appears in numerous studies related to
structure-activity correlations.

Among these models, models $1$, $2$, and $4$ incorporate a combination of
Mulliken and Lowdin charges. Despite their shared purpose of describing the
same physical quantities, these two charge schemes have distinct mathematical
foundations. Consequently, we consider only model 3 as an appropriate approach
for correlating biological activity with electronic indices. Taking into
account the aforementioned observations, we can conclude that descriptor 3
yields the most favorable outcome among all the descriptors used in the
machine learning methodology.

Figure 3 illustrates a two-dimensional map of the descriptors, which rely on
Mulliken’s charge description and the differences in orbital energy values
($\Delta H$). This model is specifically tailored for the training set, with
four molecules chosen for the test set. Notably, there is a well-defined
separation between the active and inactive groups, with no observed
intersections.

Figure 3: (a) A two-dimensional map of machine learning results using
descriptor 3 and the training set, and (b) a selected zoomed part of the map.
The final descriptors are plotted on the axes, and it is evident that the
groups do not overlap in the border region. The system metric indicates a
negative overlap (separation) of -3.20.

Figure 4 presents the ML model, based on the training set, expanded with
externally calculated descriptor values for the test set molecules. These
additional test set data points slightly alter the boundary definitions as
they exhibit a closer alignment with their respective groups. Remarkably, it
is worth emphasizing that the horizontal descriptor alone is capable of
achieving 100% accuracy in classifying the groups into active and inactive
categories.

Figure 4: (a) A two-dimensional map of machine learning results, including an
externally calculated test set, and (b) a selected zoomed region of the map.
The figure demonstrates that the model effectively describes the groups, even
though the molecules only slightly redefine the boundaries.

In Figure 5, we present the map of descriptors obtained through our
reimplementation of machine learning using the combined dataset, consisting of
both the training set and the test set. Notably, the fundamental
characteristics of the model remain unchanged, but we observe a slight change
in the boundary definitions. More specifically, when considering the model
with all entries, the intersection of the groups shifts from -3.20 to -3.05.
This shift indicates that there is no significant overfitting, suggesting that
our model is not excessively reliant on specific data points. Additionally, we
find that the descriptor (HOMO-1-HOMO)/(M1-M3) continues to effectively be
able to distinguish between the two distinct groups of molecules: active and
inactive. This is evident from the horizontal division observed in the inset
of Figure 5.

Figure 5: (a) A two-dimensional map of machine learning descriptors, including
the test set molecules (8, 9, 10, and 11), calculated separately, and (b) a
selected zoomed region of the map. The dashed line clearly indicates a group
separation, primarily observed with the horizontal axis descriptor. The
molecules previously discarded are color-highlighted.

The physical interpretation of the factor $\Delta H$ in the context of
structure-activity relationships still remains somewhat elusive, despite its
frequent use in studies involving the IC50 index and other biological
quantifiers 19, 20, 21, 22, 5. The HOMO energy, commonly associated with
ionization potential and nucleophilic regions indicative of enhanced
reactivity, is significantly important. This is further supported by the
presence of a charge difference within the descriptor (M1-M3), precisely
located in highly reactive regions, in general, associated with reactive bond
sites.

Notably, the ML horizontal descriptor alone can with which accuracy classify
the active/inactive molecules. This validates that a simple rule is able to
establish the relationship between molecular structure and biological
activity. To further validate this observation, we build a Boolean decision
tree (Figure 6). By applying the following two rules, we achieved a 100%
accuracy in the separation of inactive and active groups:

  * •

Condition 1: The charge difference between atoms 1 and 3 (M1 and M3) is less
than -0.002.

  * •

Condition 2: The energy value of the penultimate occupied molecular orbital
(HOMO-1) is higher than -9.32.

According to condition 1, we were able to successfully differentiate all
molecules except for numbers 10, 19, and 22, which were mistakenly categorized
as inactive. However, applying the second rule, we rectified this
misclassification by reassigning these molecules to the appropriate active
group.

Figure 6: Simplified Boolean tree diagram. The diagram illustrates that two
conditions are adequate to entirely segregate the two classes of molecules.

##  4 Conclusions

In summary, the difference in charges (M1-M3) observed in both descriptors
indicates the relevance of the reactivity in the aromatic region between atoms
1 and 3 in influencing the biological behavior of these molecules. Moreover,
the horizontal descriptor proves adequate in determining group behavior, while
the difference between (HOMO-1-HOMO) energy values was also identified as a
crucial parameter for our model. Although we cannot assert with absolute
certainty that our model is definitive, owing to limitations in sample size,
the results remain consistent with other studies in the field, further
underscoring the pivotal role of $\Delta H$ in structure-activity
correlations. Our study presents a streamlined methodology for optimizing the
search process for molecules with specific biological activity, exemplified
here for the chalcone skeleton. By using common skeletal structures to design
new molecules (for instance, incorporating different chemical groups), we can,
in principle, predict the activity/inactivity (based on the IC50 index) of
these novel molecules through the parameter features identified by the ML
models. Notably, these predictions rely solely on electronic structure
factors, utilizing the [(HOMO${-1}$-HOMO)/(M1-M3)] descriptor calculated via
the PM3 method based. The methodology proposed here is completely general and
can used for other classes of compounds. Studies along these lines are in
progress.

The authors would like to thank CNPq (Grants 308428/2022-6 and 150595/2023-9).
The authors would also like to thank CCM-UFABC for the computational
resources.

Figure 7 show the statistical analysis. In this figure, each pixel corresponds
to the Pearson coefficient (CP), which represents the correlation between two
features.

Figure 7: Pearson coefficient (CP). Each pixel represents the correlation
between two features.

## References

  * Lombardino and Lowe III 2004 Lombardino, J. G.; Lowe III, J. A. The role of the medicinal chemist in drug discovery—then and now. _Nature Reviews Drug Discovery_ 2004, _3_ , 853–862
  * Azam and Abbasi 2013 Azam, S. S.; Abbasi, S. W. Molecular docking studies for the identification of novel melatoninergic inhibitors for acetylserotonin-O-methyltransferase using different docking routines. _Theoretical Biology and Medical Modelling_ 2013, _10_ , 1–16
  * Novak and Kovač 2010 Novak, I.; Kovač, B. Electronic structure and biological activity: Barbiturates vs. thiobarbiturates. _Chemical Physics Letters_ 2010, _493_ , 242–244
  * Novak and Kovač 1999 Novak, I.; Kovač, B. Electronic structure and biological activity of steroids. _Biophysical Chemistry_ 1999, _78_ , 233–240
  * Braga and Galvão 2003 Braga, S. F.; Galvão, D. S. A Structure- Activity Study of Taxol, Taxotere, and Derivatives Using the Electronic Indices Methodology (EIM). _Journal of chemical information and computer sciences_ 2003, _43_ , 699–706
  * Autreto and Lavarda 2008 Autreto, P. A. d. S.; Lavarda, F. C. Febrifugine derivative antimalarial activity: quantum mechanical predictors. _Revista do Instituto de Medicina Tropical de São Paulo_ 2008, _50_ , 21–24
  * Barone et al. 1996 Barone, P.; Camilo Jr, A.; Galvão, D. Theoretical approach to identify carcinogenic activity of polycyclic aromatic hydrocarbons. _Physical review letters_ 1996, _77_ , 1186
  * Braga et al. 1999 Braga, R.; Barone, P.; Galvao, D. Identifying carcinogenic activity of methylated polycyclic aromatic hydrocarbons (PAHs). _Journal of Molecular Structure: THEOCHEM_ 1999, _464_ , 257–266
  * Schleder et al. 2019 Schleder, G. R.; Padilha, A. C. M.; Acosta, C. M.; Costa, M.; Fazzio, A. From DFT to machine learning: recent approaches to materials science–a review. _Journal of Physics: Materials_ 2019, _2_ , 032001
  * Butler et al. 2018 Butler, K. T.; Davies, D. W.; Cartwright, H.; Isayev, O.; Walsh, A. Machine learning for molecular and materials science. _Nature_ 2018, _559_ , 547–555
  * Sakata et al. 2017 Sakata, R. P.; Figueiro, M.; Kawano, D. F.; Almeida, W. P. Effect on acetylcholinesterase and anti-oxidant activity of synthetic chalcones having a good predicted pharmacokinetic profile. _Medicinal Chemistry_ 2017, _13_ , 654–663
  * Ellman et al. 1961 Ellman, G. L.; Courtney, K.; Andres, V.; Featherstone, R. M. A new and rapid colorimetric determination of acetylcholinesterase activity. _Biochemical Pharmacology_ 1961, _7_ , 88–95
  * Neese 2012 Neese, F. The ORCA program system. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2012, _2_ , 73–78
  * Fischer 1977 Fischer, C. F. Hartree–Fock method for atoms. A numerical approach. 1977,
  * Szabo and Ostlund 2012 Szabo, A.; Ostlund, N. S. _Modern quantum chemistry: introduction to advanced electronic structure theory_ ; Courier Corporation, 2012
  * Ramachandran et al. 2008 Ramachandran, K.; Deepa, G.; Namboori, K. _Computational chemistry and molecular modeling: principles and applications_ ; Springer Science & Business Media, 2008
  * Schäfer et al. 1992 Schäfer, A.; Horn, H.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets for atoms Li to Kr. _The Journal of Chemical Physics_ 1992, _97_ , 2571–2577
  * Ouyang et al. 2018 Ouyang, R.; Curtarolo, S.; Ahmetcik, E.; Scheffler, M.; Ghiringhelli, L. M. SISSO: A compressed-sensing method for identifying the best low-dimensional descriptor in an immensity of offered candidates. _Physical Review Materials_ 2018, _2_ , 083802
  * Barone et al. 2000 Barone, P.; Braga, R.; Camilo Jr, A.; Galvão, D. Electronic indices from semi-empirical calculations to identify carcinogenic activity of polycyclic aromatic hydrocarbons. _Journal of Molecular Structure: THEOCHEM_ 2000, _505_ , 55–66
  * Braga and Galvão 2004 Braga, S. F.; Galvão, D. S. Benzo [c] quinolizin-3-ones theoretical investigation: SAR analysis and application to nontested compounds. _Journal of chemical information and computer sciences_ 2004, _44_ , 1987–1997
  * Coluci et al. 2002 Coluci, V. R.; Vendrame, R.; Braga, R.; Galvao, D. Identifying relevant molecular descriptors related to carcinogenic activity of polycyclic aromatic hydrocarbons (PAHs) using pattern recognition methods. _Journal of chemical information and computer sciences_ 2002, _42_ , 1479–1489
  * Troche et al. 2005 Troche, K. S.; Braga, S. F.; Coluci, V. R.; Galvão, D. S. Carcinogenic classification of polycyclic aromatic hydrocarbons through theoretical descriptors. _International journal of quantum chemistry_ 2005, _103_ , 718–730