The mass and spectral function of scalar and pseudoscalar mesons in a hot and chirally imbalanced medium using the two-flavor NJL model

Quantum chromodynamics (QCD) is a non-Abelian gauge theory characterized by a gauge group $SU(3)$. It provides a framework for studying the strong interactions occurring between quarks and gluons, originating through the exchange of color charges. The study of the QCD vacuum structure under extreme conditions of temperature and/or baryon density stands as a central objective in the context of relativistic heavy ion collision (HIC) experiments at RHIC and LHC. It is well established that the numerous energy-degenerate vacuum configurations of QCD at zero and low temperatures can be identified through topologically non-trivial gauge configurations with a non-zero winding number Shifman (1989). These particular gluon configurations are referred to as instantons, and they have the capability to induce transitions between distinct vacua by crossing potential barriers with heights approximately on the order of the QCD scale $\Lambda_{\text{QCD}}$. This phenomenon is known as instanton tunneling Belavin et al. (1975); ’t Hooft (1976a, b). However, at high temperatures, particularly in the quark-gluon plasma (QGP) phase observed in HICs, an abundant production of another class of gluon configurations, known as sphalerons, is anticipated Manton (1983); Klinkhamer and Manton (1984). It is conjectured that the prevalence of sphalerons can enhance the transition rate across energy-degenerate vacuum states by surmounting the energy barriers Kuzmin et al. (1985); Arnold and McLerran (1987); Khlebnikov and Shaposhnikov (1988); Arnold and McLerran (1988). The gauge field configurations with topological non-triviality have the capability to alter the helicities of quarks during interactions, thus leading to the violation of parity ($P$) and charge-parity ($CP$) symmetries by generating an asymmetry between left and right-handed quarks, as a consequence of the axial anomaly of QCD Adler (1969); Bell and Jackiw (1969). The production of chirality imbalance can occur locally, as there is no direct evidence of the global violation of $P$ and $CP$ in QCD Adler (1969); Bell and Jackiw (1969); McLerran et al. (1991); Moore and Tassler (2011). This locally induced chirality imbalance is characterized using a chiral chemical potential (CCP), which essentially quantifies the difference between the numbers of right and left-handed quarks.

In non-central heavy-ion collisions (HICs), magnetic fields of very high intensity, of the order of a few $m_{\pi}^{2}$, can arise Kharzeev et al. (2008); Skokov et al. (2009). When such strong magnetic fields coexist with chirality imbalance, they can result in the separation of positive and negative charges with respect to the reaction plane, giving rise to an electric current along the magnetic field direction. This phenomenon is known as the chiral magnetic effect (CME) Fukushima et al. (2008); Kharzeev et al. (2008); Kharzeev and Warringa (2009); Bali et al. (2012). Extensive efforts have been dedicated to the detection of the CME in HIC experiments conducted at RHIC at Brookhaven. However, a recent analysis by the STAR Collaboration has not provided any evidence of the CME occurring in these collisions Abdallah et al. (2021). Consequently, new experimental techniques for detecting the CME have been proposed An et al. (2022); Milton et al. (2021); Kharzeev et al. (2022); Abdulhamid et al. (2023).

Given the anticipation of locally generated chirality imbalance in QGP, in addition to the CME, considerable research endeavours have been devoted to understand various facets of this phenomenon. These studies encompass investigations into microscopic transport phenomena  Vilenkin (1979, 1980); Fukushima et al. (2008); Son and Surowka (2009), collective oscillations  Akamatsu and Yamamoto (2013); Carignano and Manuel (2019, 2021), fermion damping rate Carignano and Buballa (2020) and collisional energy loss of fermions Carignano and Manuel (2021) within chirally imbalanced media. Furthermore, chirally asymmetric plasmas are anticipated to form in the gap regions of the magnetospheres of pulsars and black holes Gorbar and Shovkovy (2021), as well as in other stellar astrophysical scenarios Charbonneau and Zhitnitsky (2010); Akamatsu and Yamamoto (2013); Kaminski et al. (2016); Yamamoto (2016); Shovkovy (2021). It’s worth noting that the CME has been observed in condensed matter systems, particularly in 3D Dirac and Weyl semimetals Li et al. (2016); Li and Kharzeev (2016); Kharzeev (2014); Kharzeev et al. (2016); Huang (2016); Landsteiner (2016); Gorbar et al. (2018); Joyce and Shaposhnikov (1997); Tashiro et al. (2012). As a result, the investigation of the properties of chirally imbalanced matter remains a topic of significant and ongoing interest.

However, owing to the non-perturbative nature of QCD it is difficult to address these problems from first principles calculations. Most of our current understanding comes from Lattice QCD (LQCD) simulations de Forcrand and Philipsen (2007); Aoki et al. (2006, 2009); Bazavov et al. (2009); Cheng et al. (2008); Muroya et al. (2003). But, direct simulation of QCD thermodynamics at finite density is not possible, due to the sign problem. This arises when simulating QCD at non-zero baryon density because the fermion determinant, which is necessary for calculating observables, becomes complex and prevents the application of traditional Monte Carlo methods. An alternative for probing low-energy QCD involves the utilization of effective models. These models aim to isolate the relevant physics of a phenomenon by constructing mathematically manageable frameworks which can be used to probe the regions that are beyond the scope of perturbative QCD or LQCD. Often, these effective models are parametrized through fitting with more reliable theoretical frameworks or experimental data. In this work we will employ the Nambu$-$Jona-Lasinio (NJL) model Nambu and Jona-Lasinio (1961a, b) which respects the global symmetries of QCD, most importantly the chiral symmetry. Being an effective model, within its domain of applicability it provides a useful scheme to study some of the important non-perturbative properties of the QCD vacuum  Klevansky (1992); Vogl and Weise (1991); Buballa (2005); Volkov and Radzhabov (2006). As the gluonic degrees of freedom have been integrated out in this effective description, point-like interaction among the quarks appears Klevansky (1992) which makes the NJL model non-renormalizable. Therefore, one has to choose a proper regularization scheme to deal with the divergent integrals and subsequently, fix the model parameters by reproducing a set of phenomenological quantities, for example the pion-decay constant, quark condensate etc. This model has been extensively used to study the phase structure Ruggieri and Peng (2016); Ruggieri et al. (2016, 2020); Chaudhuri et al. (2021) as well as properties of electromagnetic spectral function Ghosh et al. (2022), dilepton production rate Chaudhuri et al. (2022) in a chirally asymmetric medium.

It is well known that, mesons are the bound/resonant states of quarks and antiquarks, so their propagation can be studied from the scattering of quarks/antiquarks in different channels employing the Bethe-Salpeter formalism. The physical properties of light scalar and pseudoscalar mesons have a direct relevance with the dynamics of chiral phase transition. In literature a lot of works can be found where NJL model is used to study different mesonic modes at finite temperature and/or density Hatsuda and Kunihiro (1985); Ayala et al. (2002); Hansen et al. (2007); Inagaki et al. (2008). However, we have not encountered any references that investigates spectral properties of mesons in the presence of chiral imbalance. Polarization functions of mesons play a crucial role in the computation of meson propagators within the framework of the Random Phase Approximation (RPA) Klevansky (1992). These propagators are instrumental in calculating quark and anti-quark in medium scattering cross-sections, which are essential for determining various transport coefficients using NJL-like models Zhuang et al. (1995) which involve meson exchange. Consequently, these polarization functions facilitate the study of dissipative phenomena in strongly interacting matter in presence of chiral imbalance.

The paper is organized as follows. In Sec. II we have discussed the gap equation in the presence of CCP. The meson propagators and spectral functions are derived in Sec. III. We present numerical results and related discussions in Sec. IV followed by summary and conclusions in Sec. V.

Let us start with the standard expression of the Lagrangian (density) for a 2-flavor NJL model at finite CCP $\mu_{5}$ and quark chemical potential $\mu$ as Nambu and Jona-Lasinio (1961a, b); Ghosh et al. (2022); Chaudhuri et al. (2022)

$\displaystyle\mathscr{L}_{\text{NJL}}=\overline{\psi}(i\cancel{\partial}-m+\mu\gamma^{0}+\mu_{5}\gamma^{0}\gamma^{5})\psi+G_{s}\left[(\overline{\psi}\psi)^{2}+(\overline{\psi}i\gamma^{5}\bm{\tau}\psi)^{2}\right]$

(1)

where $\psi=\left(u\leavevmode\nobreak\ d\right)^{T}$ is the quark iso-doublet, $m$ is the current quark mass, $\bm{\tau}$ is the Pauli isospin matrices and $G_{s}$ represents the scalar coupling. Throughout the article, we will be using metric tensor $g^{\mu\nu}=\text{diag}(1,-1,-1,-1)$.

In the Hartree mean field approximation (MFA), the constituent quark mass $M$ is obtained by solving the gap equation

$\displaystyle M=m-2G_{s}\left\langle\overline{\psi}\psi\right\rangle$

(2)

where, the chiral condensate is given by

$\displaystyle\left\langle\overline{\psi}\psi\right\rangle=\text{Re}\leavevmode\nobreak\ i\int\frac{d^{4}p}{(2\pi)^{4}}\text{Tr}_{\text{d,c,f}}\left[S_{11}(p)\right]$

(3)

in which the $\text{Tr}_{\text{d,c,f}}[...]$ implies taking trace in Dirac, color and flavor spaces and $S_{11}(p)$ is the $11$ component of the real time thermal quark propagator defined in the context of real time formulation (RTF) of finite temperature field theory. The explicit expression of $S_{11}(p)$ is given by Ghosh et al. (2022); Chaudhuri et al. (2022)

$\displaystyle S_{11}(p)=\mathcal{D}(p^{0}+\mu,\bm{p})\sum_{r\in\{\pm\}}\frac{1}{4|\bm{p}|r\mu_{5}}\left[\frac{-1}{(p^{0}+\mu)^{2}-(\omega^{r}_{\bm{p}})^{2}+i\epsilon}-2\pi i\eta(p^{0}+\mu)\delta\left\{(p^{0}+\mu)^{2}-(\omega^{r}_{\bm{p}})^{2}\right\}\right]\otimes\mathds{1}_{\text{color}}\otimes\mathds{1}_{\text{flavor}}$

(4)

where, $r$ represents the helicity of the propagating quark, $\omega^{r}_{\bm{p}}=\sqrt{(|\bm{p}|+r\mu_{5})^{2}+M^{2}}$ is the single particle dispersion in presence of CCP, $\eta(x)=\theta(x)f_{+}(x)+\theta(-x)f_{-}(-x)$, $\theta(x)$ is the unit step function and $f_{\pm}(x)=\dfrac{1}{e^{(x\mp\mu)/T}+1}$ is the Fermi-Dirac thermal distribution function and $\mathcal{D}$ contains complicated Dirac structure as follows

$\displaystyle\mathcal{D}(p^{0},\bm{p})=\sum_{j\in\{\pm\}}\mathscr{P}_{j}\left[p_{-j}^{2}\cancel{p}_{j}-M^{2}\cancel{p}_{-j}+M(p_{j}\cdot p_{-j}-M^{2})+iM\sigma_{\mu\nu}p_{j}^{\mu}p_{-j}^{\nu}\right]$

(5)

in which $\mathscr{P}_{j}=\frac{1}{2}(\mathds{1}+j\gamma^{5})$ and $p^{\mu}_{j}\equiv(p^{0}+j\mu_{5},\bm{p})$. Substituting the propagator from Eq. (4) into Eq. (3) followed by performing the $dp^{0}$ integration, we get after some simplifications the chiral condensate as

$\displaystyle\left\langle\overline{\psi}\psi\right\rangle=-N_{c}N_{f}M\sum_{r\in\{\pm\}}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{\omega^{r}_{\bm{p}}}\left(1-f_{+}^{\bm{p}r}-f_{-}^{\bm{p}r}\right)$

(6)

where, $N_{c}=3$ is the number of colors, $N_{f}=2$ is the number of flavors and $f_{\pm}^{\bm{p}r}=f_{\pm}(\omega^{r}_{\bm{p}})$. The first term within the parenthesis in the above integral is the temperature independent vacuum contribution to the chiral condensate which is UV divergent; and as the NJL model is known to be non-renormalizable, a proper regularization prescription has to be employed to tackle the divergences. Various regularization schemes have been explored in the literature. One of the widely adopted approaches involves the use of a physically intuitive three-momentum cutoff, which is noncovariant in nature. This cutoff can be applied in two ways: by employing the regulator solely in the temperature-independent integrals or by applying it in both the vacuum and thermal parts, even if the thermal part is already UV finite. Additionally, several covariant schemes have been introduced, such as the four-momentum cutoff in Euclidean space, regularization in proper time, and the Pauli-Villars method. A comprehensive discussion on different regularization schemes can be found in references such as Klevansky (1992); Kohyama et al. (2015); Xue et al. (2022). In this work we will be using the smooth three-momentum cutoff regulator only in the divergent vacuum term following Ref. Fukushima et al. (2010) as the sharp momentum cutoff suffers from a regularization artifact. This is implemented by replacing the vacuum term as

$\displaystyle\int\dfrac{d^{3}p}{(2\pi)^{3}}1\to\int\dfrac{d^{3}p}{(2\pi)^{3}}f_{\Lambda}^{\bm{p}}$

(7)

where $f_{\Lambda}^{\bm{p}}=\sqrt{\dfrac{\Lambda^{20}}{\Lambda^{20}+|\bm{p}|^{20}}}$. Substituting Eq. (6) in Eq. (2) and making use of Eq. (7), we finally obtain the gap equation as

$\displaystyle M=m+2G_{s}N_{c}N_{f}M\sum_{r\in\{\pm\}}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{\omega^{r}_{\bm{p}}}\left(f_{\Lambda}^{\bm{p}}-f_{+}^{\bm{p}r}-f_{-}^{\bm{p}r}\right).$

(8)

Mesons are the bound/resonant states of quarks and antiquarks, so their propagation can be studied from the scattering of quarks/antiquarks in different channels using the Bethe-Salpeter formalism Hatsuda and Kunihiro (1985); Klevansky (1992); Ayala et al. (2002); Hansen et al. (2007); Inagaki et al. (2008); Ghosh et al. (2020). Using the random phase approximation (RPA), the meson propagators $D_{h}$ in the scalar ($\pi$) and pseudo-scalar ($\sigma$) channels can be written as

$\displaystyle D_{h}(q;T,\mu_{5})=\frac{2G}{1-2G\Pi_{h}(q;T,\mu_{5})}\leavevmode\nobreak\ \leavevmode\nobreak\ ;\leavevmode\nobreak\ \leavevmode\nobreak\ h\in\{\pi,\sigma\}$

(9)

where, $\Pi_{h}$ is the one-loop polarization function of meson $h$.

The one-loop polarization function $\Pi_{h}(q)$ at finite CCP and temperature can be calculated using the RTF of finite temperature field theory. Being a two-point correlation function, the real time polarization function acquires a $2\times 2$ matrix structure (in thermal space) whose $11$-component is given by

$\displaystyle\Pi_{h}^{11}(q)=i\int\!\!\frac{d^{4}k}{(2\pi)^{4}}\text{Tr}_{\text{d,c,f}}\left[\Gamma_{h}S_{11}(p=q+k)\Gamma_{h}S_{11}(k)\right]$

(10)

where, $h\in\{\pi,\sigma\}$, $\Gamma_{\pi}=i\gamma^{5}$, $\Gamma_{\sigma}=\mathds{1}$ and the quark propagator $S_{11}(p)$ is defined in Eq. (4). Now substituting Eq. (4) into Eq. (10) and by performing the $dk^{0}$ integration, we obtain the real and imaginary parts of the 11-component of mesonic polarization function as

	$\displaystyle\text{Re}\Pi_{h}^{11}(q_{0},\bm{q})$	$\displaystyle=$	$\displaystyle N_{c}N_{f}\int\frac{d^{3}k}{(2\pi)^{3}}\sum_{r\in\{\pm\}}\sum_{s\in\{\pm\}}\frac{1}{64rs\mu_{5}^{2}|\bm{p}||\bm{k}|}\mathcal{P}\Big{[}\frac{\mathcal{N}_{h}(k^{0}=\omega^{r}_{\bm{k}})(f_{\Lambda}^{\bm{k}}-2f_{+}^{\bm{k}r})}{\omega^{r}_{\bm{k}}\left\{(q^{0}+\omega^{r}_{\bm{k}})^{2}-(\omega^{s}_{\bm{p}})^{2}\right\}}+\frac{\mathcal{N}_{h}(k^{0}=-\omega^{r}_{\bm{k}})(f_{\Lambda}^{\bm{k}}-2f_{-}^{\bm{k}r})}{\omega^{r}_{\bm{k}}\left\{(q^{0}-\omega^{r}_{\bm{k}})^{2}-(\omega^{s}_{\bm{p}})^{2}\right\}}$		(11)
			$\displaystyle+\frac{\mathcal{N}_{h}(k^{0}=-q^{0}+\omega^{s}_{\bm{p}})(f_{\Lambda}^{\bm{k}}-2f_{+}^{\bm{p}s})}{\omega^{s}_{\bm{p}}\left\{(q^{0}-\omega^{s}_{\bm{p}})^{2}-(\omega^{r}_{\bm{k}})^{2}\right\}}+\frac{\mathcal{N}_{h}(k^{0}=-q^{0}-\omega^{s}_{\bm{p}})(f_{\Lambda}^{\bm{k}}-2f_{-}^{\bm{p}s})}{\omega^{s}_{\bm{p}}\left\{(q^{0}+\omega^{s}_{\bm{p}})^{2}-(\omega^{r}_{\bm{k}})^{2}\right\}}\Big{]},$
	$\displaystyle\text{Im}\Pi_{h}^{11}(q_{0},\bm{q})$	$\displaystyle=$	$\displaystyle-N_{c}N_{f}\pi\int\frac{d^{3}k}{(2\pi)^{3}}\sum_{r\in\{\pm\}}\sum_{s\in\{\pm\}}\frac{1}{16rs\mu_{5}^{2}|\bm{p}||\bm{k}|}\frac{1}{4\omega^{r}_{\bm{k}}\omega^{s}_{\bm{p}}}$		(12)
			$\displaystyle\times\Big{[}\mathcal{N}_{h}(k^{0}=-\omega^{r}_{\bm{k}})(1-f_{-}^{\bm{k}r}-f_{+}^{\bm{p}s}+2f_{-}^{\bm{k}r}f_{+}^{\bm{p}s})\delta(q_{0}-\omega^{r}_{\bm{k}}-\omega^{s}_{\bm{p}})+\mathcal{N}_{h}(k^{0}=\omega^{r}_{\bm{k}})(1-f_{+}^{\bm{k}r}-f_{-}^{\bm{p}s}+2f_{+}^{\bm{k}r}f_{-}^{\bm{p}s})\delta(q_{0}+\omega^{r}_{\bm{k}}+\omega^{s}_{\bm{p}})$
			$\displaystyle+\mathcal{N}_{h}(k^{0}=\omega^{r}_{\bm{k}})(-f_{+}^{\bm{k}r}-f_{+}^{\bm{p}s}+2f_{+}^{\bm{k}r}f_{+}^{\bm{p}s})\delta(q_{0}+\omega^{r}_{\bm{k}}-\omega^{s}_{\bm{p}})+\mathcal{N}_{h}(k^{0}=-\omega^{r}_{\bm{k}})(-f_{-}^{\bm{k}r}-f_{-}^{\bm{p}s}+2f_{-}^{\bm{k}r}f_{-}^{\bm{p}s})\delta(q_{0}-\omega^{r}_{\bm{k}}+\omega^{s}_{\bm{p}})\Big{]}$

where, the smooth-cutoff $f_{\Lambda}^{\bm{k}}$ has been used to regulate the UV-divergent part of $\text{Re}\Pi_{h}^{11}(q_{0},\bm{q})$, $\mathcal{P}$ denotes Cauchy principal value integration, and

	$\displaystyle\mathcal{N}_{h}(k)$	$\displaystyle=$	$\displaystyle 4a_{h}M^{6}+2M^{4}\big{\{}-2a_{h}(k_{+}\cdot k_{-})+(k_{+}\cdot p_{-})+(k_{-}\cdot p_{+})-2a_{h}(p_{+}\cdot p_{-})\big{\}}-2M^{2}\big{\{}-2a_{h}(k_{+}\cdot p_{-})(k_{+}\cdot p_{+})$		(13)
			$\displaystyle-2a_{h}(k_{+}\cdot k_{-})(p_{+}\cdot p_{-})+2a_{h}(k_{+}\cdot p_{+})(k_{-}\cdot p_{-})+(k_{-}\cdot p_{-})(k_{+}^{2}+p_{+}^{2})+(k_{+}\cdot p_{+})(k_{-}^{2}+p_{-}^{2})\big{\}}$
			$\displaystyle+2\big{\{}(k_{+}\cdot p_{-})k_{-}^{2}p_{+}^{2}+(k_{-}\cdot p_{+})k_{+}^{2}p_{-}^{2}\big{\}}$

in which $a_{\pi}=-1$ and $a_{\sigma}=1$. Having calculated the 11-components, it is now easy to diagonalize the real-time thermal polarization matrix Mallik and Sarkar (2016), and extract the diagonal components which are analytic functions that appear in Eq. (9). The real and imaginary parts of the analytic function $\Pi_{h}(q)$ are related to the 11-components via relations

$\displaystyle\text{Re}\Pi_{h}(q)=\text{Re}\Pi_{h}^{11}(q)\leavevmode\nobreak\ \leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \text{Im}\Pi_{h}(q)=\text{sign}(q^{0})\tanh\left(\frac{q^{0}}{2T}\right)\text{Im}\Pi_{h}^{11}(q)$

(14)

respectively. For the calculation of meson masses, we only require the special cases: $\Pi_{h}(q^{0},\bm{q}=\bm{0})$. For $q^{0}\neq 0,\bm{q}=\bm{0}$, substitution of Eqs. (11) and (12) into Eq. (14) yields after simplification:

	$\displaystyle\text{Re}\Pi_{h}(q_{0},\bm{q}=\bm{0})$	$\displaystyle=$	$\displaystyle\frac{N_{c}N_{f}}{2\pi^{2}}\int_{0}^{\infty}d|\bm{k}|\sum_{r\in\{\pm\}}\sum_{s\in\{\pm\}}\frac{1}{64rs\mu_{5}^{2}}\mathcal{P}\Big{[}\frac{\mathcal{N}_{h}(k^{0}=\omega^{r}_{\bm{k}})(f_{\Lambda}^{\bm{k}}-2f_{+}^{\bm{k}r})}{\omega^{r}_{\bm{k}}\left\{(q^{0}+\omega^{r}_{\bm{k}})^{2}-(\omega^{s}_{\bm{k}})^{2}\right\}}+\frac{\mathcal{N}_{h}(k^{0}=-\omega^{r}_{\bm{k}})(f_{\Lambda}^{\bm{k}}-2f_{-}^{\bm{k}r})}{\omega^{r}_{\bm{k}}\left\{(q^{0}-\omega^{r}_{\bm{k}})^{2}-(\omega^{s}_{\bm{k}})^{2}\right\}}$		(15)
			$\displaystyle+\frac{\mathcal{N}_{h}(k^{0}=-q^{0}+\omega^{s}_{\bm{k}})(f_{\Lambda}^{\bm{k}}-2f_{+}^{\bm{k}s})}{\omega^{s}_{\bm{k}}\left\{(q^{0}-\omega^{s}_{\bm{k}})^{2}-(\omega^{r}_{\bm{k}})^{2}\right\}}+\frac{\mathcal{N}_{h}(k^{0}=-q^{0}-\omega^{s}_{\bm{k}})(f_{\Lambda}^{\bm{k}}-2f_{-}^{\bm{k}s})}{\omega^{s}_{\bm{k}}\left\{(q^{0}+\omega^{s}_{\bm{k}})^{2}-(\omega^{r}_{\bm{k}})^{2}\right\}}\Big{]}_{\bm{q}=\bm{0}},$
	$\displaystyle\text{Im}\Pi_{h}(q_{0},\bm{q}=\bm{0})$	$\displaystyle=$	$\displaystyle-\frac{N_{c}N_{f}}{2\pi}\int_{0}^{\infty}d|\bm{k}|\sum_{r\in\{\pm\}}\sum_{s\in\{\pm\}}\frac{1}{16rs\mu_{5}^{2}}\frac{1}{4\omega^{r}_{\bm{k}}\omega^{s}_{\bm{k}}}$		(16)
			$\displaystyle\times\Big{[}\mathcal{N}_{h}(k^{0}=-\omega^{r}_{\bm{k}})(1-f_{-}^{\bm{k}r}-f_{+}^{\bm{k}s}+2f_{-}^{\bm{k}r}f_{+}^{\bm{k}s})\delta(q_{0}-\omega^{r}_{\bm{k}}-\omega^{s}_{\bm{k}})+\mathcal{N}_{h}(k^{0}=\omega^{r}_{\bm{k}})(1-f_{+}^{\bm{k}r}-f_{-}^{\bm{k}s}+2f_{+}^{\bm{k}r}f_{-}^{\bm{k}s})\delta(q_{0}+\omega^{r}_{\bm{k}}+\omega^{s}_{\bm{k}})$
			$\displaystyle+\mathcal{N}_{h}(k^{0}=\omega^{r}_{\bm{k}})(-f_{+}^{\bm{k}r}-f_{+}^{\bm{k}s}+2f_{+}^{\bm{k}r}f_{+}^{\bm{k}s})\delta(q_{0}+\omega^{r}_{\bm{k}}-\omega^{s}_{\bm{k}})+\mathcal{N}_{h}(k^{0}=-\omega^{r}_{\bm{k}})(-f_{-}^{\bm{k}r}-f_{-}^{\bm{k}s}+2f_{-}^{\bm{k}r}f_{-}^{\bm{k}s})\delta(q_{0}-\omega^{r}_{\bm{k}}+\omega^{s}_{\bm{k}})\Big{]}_{\bm{q}=\bm{0}}.$

We note that, the remaining $d|\bm{k}|$ integral of Eq. (16) can be performed using the Dirac delta functions present in the integrand.

$\text{Im}\Pi_{h}(q_{0},\bm{q}=\bm{0})$ in Eq. (16) contains sixteen Dirac delta functions and they give rise to branch cuts of the polarization function in the complex energy plane. The terms with $\delta(q_{0}-\omega^{r}_{\bm{k}}-\omega^{s}_{\bm{k}})$ and $\delta(q_{0}+\omega^{r}_{\bm{k}}+\omega^{s}_{\bm{k}})$ are the Unitary-I and Unitary-II cuts respectively, whereas the terms with $\delta(q_{0}+\omega^{r}_{\bm{k}}-\omega^{s}_{\bm{k}})$ and $\delta(q_{0}-\omega^{r}_{\bm{k}}+\omega^{s}_{\bm{k}})$ are respectively termed as Landau-I and Landau-II cuts. Each of these Unitary and Landau cuts in turn consists of four sub-cuts corresponding to different helicities $(r,s)$. Each of the sixteen Dirac delta functions in Eq. (16) are non-vanishing at different respective kinematic domains which are tabulated in Eq. (27)

$\displaystyle\begin{tabular}[]{|c|c|}\hline\cr Dirac Delta Function&\noindent\hbox{}\hfill{{\hbox{\begin{tabular}[c]{@{}c@{}}\nobreak\leavevmode\hfil\\ Kinematic Regions\\ \nobreak\leavevmode\hfil\end{tabular}}\hbox{}\hfill}}\\ \hline\cr\hline\cr$\delta(q_{0}\mp\omega^{r}_{\bm{k}}\mp\omega^{s}_{\bm{k}})$&\noindent\hbox{}\hfill{{\hbox{\begin{tabular}[c]{@{}c@{}}$2M\leq\pm q_{0}<\infty$ \leavevmode\nobreak\ \leavevmode\nobreak\ if\leavevmode\nobreak\ \leavevmode\nobreak\ $(r,s)=(-,-)$,\\ $2\sqrt{M^{2}+\mu_{5}^{2}}\leq\pm q_{0}<\infty$ \leavevmode\nobreak\ \leavevmode\nobreak\ otherwise\nobreak\leavevmode\nobreak\leavevmode\end{tabular}}\hbox{}\hfill}}\\ \hline\cr$\delta(q_{0}\pm\omega^{r}_{\bm{k}}\mp\omega^{s}_{\bm{k}})$&\noindent\hbox{}\hfill{{\hbox{\begin{tabular}[c]{@{}c@{}}$0\leq\pm q_{0}\leq 2\mu_{5}$ \leavevmode\nobreak\ \leavevmode\nobreak\ if\leavevmode\nobreak\ \leavevmode\nobreak\ $(r,s)=(-,+)$,\\ $-2\mu_{5}\leq\pm q_{0}\leq 0$ \leavevmode\nobreak\ \leavevmode\nobreak\ if\leavevmode\nobreak\ \leavevmode\nobreak\ $(r,s)=(+,-)$\end{tabular}}\hbox{}\hfill}}\\ \hline\cr\end{tabular}\leavevmode\nobreak\ .$

(27)

In Eq. (16), when helicity indices $(r,s)$ are summed over, we find from Eq. (27) that, the kinematic domain for Unitary-I cut is $2M\leq^{0}<\infty$ whereas the same for Unitary-II cut is $-\infty<q^{0}\leq 2M$; on the other hand, both the Landau cuts lie in the region $|q^{0}|\leq 2\mu_{5}$. The complete analytic structure of the polarization function $\Pi_{h}(q^{0},\bm{q}=\bm{0})$ has been shown in Fig. 1. It is evident from Fig. 1 that, for non-zero $\mu_{5}$ the Landau cuts contribute to the physical time-like region defined in terms of $q^{0}>0$ and $q^{2}>0$ which is a purely finite CCP effect. Moreover, for small values of $\mu_{5}<2M$, there exist a forbidden gap between the Unitary and Landau cuts. However for sufficiently high values of $\mu_{5}>2M$, this forbidden gap becomes zero which is also a purely finite CCP effect.

Having calculated the real and imaginary parts of the polarization functions $\Pi_{h}(q;T,\mu_{5})$ from Eqs. (15) and (16), they can be substituted into Eq. (9) to obtain the meson propagators $D_{h}$. Using $D_{h}$, we will in turn calculate the in-medium mesonic spectral function $S_{h}$ defined as

$\displaystyle S_{h}(q;T,\mu_{5})=\frac{1}{2G\pi}\text{Im}D_{h}(q;T,\mu_{5})=\frac{1}{2G\pi}\frac{2G\text{Im}\Pi_{h}}{(1-2G\text{Re}\Pi_{h})^{2}+(2G\text{Im}\Pi_{h})^{2}}.$

(28)

The in-medium mesonic spectral function defined in Eq. (28) contains both the real and imaginary parts of the polarization function and thus it carries all the information about the propagation of the mesonic excitation in the medium. In this work, we also aim to calculate the masses of the mesons in scalar and pseudoscalar channels. However mass can not be uniquely defined corresponding to such collective excitation. For example, the mesonic mass can be defined as the position of the pole of the propagator in the complex energy plane at vanishing three-momentum in which case, the real(imaginary) part of the pole corresponds to the mass(width) of the mesons Zhuang et al. (1994); although the problem of finding a complex pole is a bit involved, it can be simplified by considering the small width approximation in which case the equations for the mass and width are decoupled Zhuang et al. (1994). Alternatively, one can also define the mesonic mass as the position of the global maximum of the spectral function at vanishing three-momentum. In this work, we will be using a hybrid method to calculate the mass of the mesons as we describe below.

At vanishing three-momenta ($\bm{q}=\bm{0}$), the position of the pole $q^{0}=M_{h}$ of the propagator [from Eq. (9)]

$\displaystyle D_{h}(q^{0},\bm{q}=\bm{0};T,\mu_{5})=\frac{2G}{1-2G\text{Re}\Pi_{h}(q^{0},\bm{q}=\bm{0};T,\mu_{5})-2iG\text{Im}\Pi_{h}(q^{0},\bm{q}=\bm{0};T,\mu_{5})}$

(29)

gives the mesonic mass $M_{h}$ which essentially requires solving the transcendental equation

$\displaystyle 1-2G\text{Re}\Pi_{h}(q^{0}=M_{h},\bm{q}=\bm{0};T,\mu_{5})=0$

(30)

for $M_{h}$. It is easy to realize that, for $\text{Im}\Pi_{h}(q^{0}=M_{h},\bm{q}=\bm{0};T,\mu_{5})\approx 0$ near the pole (small widths), the pole position coincides with the peak position of the spectral function as can be understood from Eq. (28). It has been observed that at finite CCP Eq. (30) may not admit a real solution for $M_{h}$. In such a scenario the mesonic mass will be obtained from the the position of global maximum $q^{0}=M_{h}$ of the spectral function $S_{h}(q^{0}=M_{h},\bm{q}=\bm{0};T,\mu_{5})$ at vanishing three-momenta ($\bm{q}=\bm{0}$).

Let us first specify the values of the parameters of the NJL model that are used in this work. We have used the NJL coupling $G=5.742$ GeV^-2, the current quark mass $m=5.6$ MeV, and smooth three-momentum cutoff $\Lambda=568.69$ MeV; these values of the parameters successfully reproduce the experimental/phenomenological vacuum values of pion decay constant $f_{\pi}\simeq 92$ MeV, quark-condensate $\left\langle\bar{u}u\right\rangle=\left\langle\bar{d}d\right\rangle\simeq-(245)^{3}$ MeV³, and the pion-mass $m_{\pi}\simeq 140$ MeV. Further, with these choice of the parameters, the dressed quark mass in vacuum comes out to be $M\simeq 343.6$ MeV. Although the analytical calculation in Secs. II-III have been done keeping the chemical potential $\mu$, we plan to give numerical results in this section by taking $\mu=0$ in order to see the interplay between $\mu_{5}$ and $T$. Numerical results taking all parameter $\{T,\mu,\mu_{5}\}$ will be reported elsewhere which definitely will be interesting to explore.

A general comment on the use of two flavor NJL model in the present context is in order here. It is well known that, the inclusion of Polyakov loop leads to a better description of the QCD phenomenology, in particular the thermodynamic quantities. However, as highlighted in Hansen et al. (2007); Chaudhuri et al. (2019, 2020), incorporating the Polyakov loop does not bring about any qualitative differences in the variation of the masses of neutral $\pi$ and $\sigma$ mesons as a function temperature. This is attributed to the fact that the constituent quark mass is the key quantity governing the temperature evolution of mesonic masses. Furthermore, the spectral properties of $\sigma$ and neutral $\pi$ mesons are expected to remain qualitatively unchanged even with the inclusion of strange quarks in a three-flavor model Rehberg et al. (1996); Costa et al. (2009). Instead, use of a two-flavor NJL model facilitates a less cumbersome calculation enabling a focused exploration of the interplay between $T$ and $\mu_{5}$.

The constituent quark mass $M$ has been obtained by solving the gap equation given by Eq. (8). As can be seen in Fig. 2(a), $M$ is large ($\sim 300$ MeV) in the low temperature region (irrespective of the value of $\mu_{5}$) due to the spontaneous breaking of the chiral symmetry characterized by the large values of the quark condensate. As the temperature increases, $M$ initially remains constant up to a certain temperature. Beyond this, $M$ gradually decreases, eventually reaching the values of the bare quark mass. This behaviour is a consequence of the pseudo-chiral phase transition. In the high temperature region, irrespective of the values of CCP, the constituent quark mass $M$ asymptotically approaches to the current quark mass $m$ value. If we compare the curves for different CCP of Fig. 2(a), we notice that, with the increase in CCP, $M$ also increases in the low temperature region. Moreover, the transition temperature of the pseudo-chiral phase transition decreases with the increase in CCP i.e. and the sudden drop of the constituent quark mass occurs at relatively smaller temperature for higher values of CCP. Hence, CCP has the tendency to make the chiral condensate stronger at low temperature region indicating a “chiral catalysis” (CC) where the chiral imbalance (characterized by $\mu_{5}$) catalyzes the dynamical symmetry breaking. On the other hand, an opposite effect is observed on the high temperature region in which the $\mu_{5}$ weakens the chiral condensate so that chiral symmetry is restored at a relatively lower temperature as compared to the $\mu_{5}=0$ case; thus indicating an “inverse chiral catalysis” (ICC) in which the dynamical symmetry breaking is opposed by the presence of a chiral imbalance Ghosh et al. (2022).

Having obtained the constituent quark mass, we have now calculated the chiral susceptibility

$\displaystyle\chi=\frac{1}{2G}\left(\frac{\partial M}{\partial m}-1\right)$

(31)

which is plotted as a function of temperature at different CCP in Fig. 2(b). The position of the maximum of the susceptibility corresponds to the pseuodo-chiral phase transition temperature $T_{c}$ Zhao et al. (2008); Chaudhuri et al. (2019). As can be observed in Fig. 2(b), with the increase in CCP, the position of the maximum moves towards lower values of temperature indicating an ICC as described in the previous paragraph. The explicit variation of $T_{c}$ as a function of CCP has been depicted in Fig. 2(c) from which we notice that $T_{c}$ decreases monotonically with the increase in $\mu_{5}$ owing to the ICC.

In this work, the neutral meson properties such as mass, polarization functions and spectral functions have been studied in presence of chiral imbalance using two-flavor Nambu–Jona-Lasinio model. To begin with, we have studied the the temperature dependence of the constituent quark mass for different values of CCP. It is found that CCP has the tendency to enhance the chiral condensate in the low temperature region indicating a chiral catalysis. Conversely, at higher temperatures, an opposing effect is observed in which $\mu_{5}$ weakens the chiral condensate. Consequently, chiral symmetry is restored at a relatively lower temperature compared to the case when $\mu_{5}=0$, indicating a manifestation of inverse chiral catalysis.

Both real and imaginary part of the polarization function as well as the spectral functions are studied in detail for both scalar and pseudo-scalar channels as a function of $q^{0}$ for different values of CCP for three different values of temperature capturing different stages of chiral phase transition. A comprehensive analysis of the real part of the scaled polarization function for the $\pi$ meson reveals multiple poles (three) in the propagator of the pseudo-scalar meson. This is observed for specific values of CCP, both in the chiral symmetry-broken phase and in the vicinity of the transition temperature. One of these multiple poles exhibits negative residue and thus corresponds to unphysical mesonic excitation. However, no such effect is observed for the scalar meson. Additionally, the imaginary part of the polarization function for both the $\sigma$ and $\pi$ mesons shows non-trivial Landau cut contributions in addition to the usual Unitary cut contributions. These non-trivial Landau cut contributions are purely a result of finite CCP effects. Moreover, the imaginary part of the polarization function exhibits a highly non-monotonic behavior with respect to $q^{0}$. At high temperatures, both the real and imaginary parts of the $\sigma$ and $\pi$ mesons are observed to coincide with each other. This merging indicates partial restoration of chiral symmetry. As a consequence of the behavior of the real part of pion polarization function mentioned earlier, at low and intermediate temperatures the spectral function of pions exhibits multiple poles along with the continuum stemming from the two-quark threshold for specific values of CCP. In contrast, the spectral function of sigma mesons shows a Breit-Wigner-like structure that begins from the two-quark threshold, indicating that sigmas are resonant states with finite decay widths. As chiral symmetry is partially restored at higher temperatures, both pions and sigmas appear as resonant states, and their spectral functions become identical. This phenomenon is known as the Mott transition, where the pion spectral function loses its Dirac delta structure and transforms into a Breit-Wigner form similar to that of the sigma meson. The Mott transition temperature shows slightly oscillatory and overall decreasing trend with the increase in CCP. Furthermore, it is observed that an increase in CCP leads to an increase in the width of the spectral function and a decrease in the peak. This implies that both pions and sigmas become more unstable as CCP increases, indicating an enhancement in their decay processes. The variations of pion and sigma masses as functions of temperature as well as CCP is also studied and at higher values of $T$ and $\mu_{5}$, an abrupt jump-like structure can be seen in the masses of both the mesonic modes.

S.G. is funded by the Department of Higher Education, Government of West Bengal, India. N.C., S.S. and P.R. are funded by the Department of Atomic Energy (DAE), Government of India.