Quantum chromodynamics axion in a hot and magnetized medium

It is analyzed the effects of a hot and magnetized medium on the axion mass, self-coupling and topological susceptibility when in presence of an anisotropic external magnetic field along the $z$-direction, within the Nambu--Jona-Lasinio effective model for quantum chromodynamics. The effects of both Magnetic Catalysis and Inverse Magnetic Catalysis are explicitly taken into account through appropriate matching of parameters with those from lattice Monte-Carlo numerical simulations. It is also analyzed the dependence of the results with respect to different model parametrizations in the context of the Nambu--Jona-Lasinio model.

It is analyzed the effects of a hot and magnetized medium on the axion mass, self-coupling and topological susceptibility when in presence of an anisotropic external magnetic field along the zdirection, within the Nambu-Jona-Lasinio effective model for quantum chromodynamics. The effects of both Magnetic Catalysis and Inverse Magnetic Catalysis are explicitly taken into account through appropriate matching of parameters with those from lattice Monte-Carlo numerical simulations. It is also analyzed the dependence of the results with respect to different model parametrizations in the context of the Nambu-Jona-Lasinio model.

I. INTRODUCTION
The axion is a pseudo Nambu-Goldstone boson of a spontaneously broken global Abelian symmetry [1]. The axion is considered to be the most elegant and robust solution to the absence of the charge and parity (CP) violating effects (also dubbed the strong CP problem) in quantum chromodynamics (QCD) [2,3]. It has also been considered as a prime candidate for cold dark matter given that they are very weakly coupled to baryonic matter in general, besides of possibly being extremely light (see, e.g., Ref. [4] for a recent review and references therein).
The importance of the QCD axion as a solution to the strong CP problem and its potential in explaining the dark matter abundance in the universe makes it one of the most sought out prospects beyond the particle physics Standard Model. Recent studies also show that axions can thermalize and form a Bose-Einstein condensate [5,6] which in turn indicates the relevance of finite temperature extension of the axion properties. On the other hand, the coupling of QCD axions with an external electromagnetic field was first exploited long ago in Ref. [7] to make the axion experimentally detectable. Since then various other experimental techniques involving the axion-photon coupling have been proposed [8][9][10] and results of some such experiments are also available [11,12].
Axions have also been associated with the anomalous stellar cooling problem [13,14], including neutron stars [15][16][17]. This is because axions could be produced in hot astrophysical plasmas and, thus, could take part of the energy transport in stellar objects. This is in particular a motivation for studying the properties of the QCD axion when in presence of an environment that accounts not only thermal effects (which can also be of relevance to the physics of these particles in the early universe), but also density (chemical potential) and external magnetic fields (relevant also for the physics of compact stellar objects, like neutron stars).
As a pseudo Goldstone boson, the QCD axion acquires a mass, m a , from the QCD chiral symmetry breaking and this mass is typically of O(m π f π /f a ) [2], where m π is the pion mass, f π is the pion decay constant and f a is the axion decay constant, which is proportional to the Peccei-Quinn symmetry breaking scale. We can also define an effective self-coupling, λ a , to the axion field. Furthermore, there are derivative couplings that can also be present in a system involving axions. The axion mass and the couplings (including the self-coupling) are all controlled by the scale f a . Though exact determination of the value of f a is not yet achievable, present astrophysical constraints suggest the value for f a to be in between 10 8 GeV and 10 18 GeV [13,18,19]. This makes the couplings associated with axion field extremely small, making it experimentally challenging to be probed (see, e.g., Ref. [20] and references therein for some of the experimental proposals looking for axions).
The extension of axion studies in relatively higher temperatures can be done perturbatively, using, e.g., the dilute instanton gas approximation [21]. But around and below the pseudo critical temperature, T c ∼ 170 MeV, which are particularly important in relation with the chiral symmetry breaking, there are no reliable perturbative techniques. Non-perturbative methods making use of effective models come into play in this regime of study along with the first principle calculations of Lattice QCD. In this particular work, we have chosen to deal with one of such effective model, the Nambu-Jona-Lasinio (NJL) model [22] to study the QCD axions in a hot and magnetized medium (For a recent study of axion within a chiral effective Lagrangian model, see Ref [23]). Beside its relevance to spontaneous CP violation [24][25][26][27][28], the advantage of using a two-flavor NJL model to study a system of axions is that being a quark model, it incorporates the effects of axions via the effective ('t Hooft determinant) interaction [29,30]. As for the electromagnetic coupling of the axions, for this work we will only consider an anisotropic external magnetic field along the z-direction. In the future we plan to extend the present study for other general systems, where both electric and magnetic fields can be present in the system. As an additional remark, since the energy scales we are working with are much smaller than the axion symmetry breaking scale f a , we can consider the axion field to be in its (constant) vacuum expectation value, so it behaves like a CP violating term added to the QCD action. Thus, in this sense, our derivations can be performed in a way similar like in previous studies [24][25][26][27][28].
This paper is organized as follows. In Sec. II we discuss the formalism of the NJL model when incorporating the axion field. We also explain the derivation of the thermodynamic potential for the model in the cases without and with an external anisotropic magnetic field. The way the parameters are fixed is also discussed. In Sec. III we present our results concerning the analysis of the effects of temperature and magnetic field on the axion mass, coupling constant and also on the axion topological susceptibility, which is an important observable derivable for example from lattice QCD results. Finally, in Sec. IV we have our conclusions.

II. FORMALISM
We start with a brief review of the QCD axion and how its associated field is included in the NJL model quasi-particle description of quark matter. Then, we will discuss the incorporation of both temperature and an external anisotropic magnetic field in the derivation of the thermodynamic potential for the model. While parametrizing the model, focus will also be given on the recently discovered phenomena of Magnetic Catalysis (MC) [31] and Inverse Magnetic Catalysis (IMC) [32,33] on the thermodynamic potential when in presence of a magnetized medium.

A. The axion contribution
The study of CP violation within strongly interacting matter has been a subject of extensive scrutiny over the years [24][25][26][27][28][34][35][36][37]. It is a well known fact that within the regime of strong interaction instanton contributions can lead to CP violation [29,30]. In this kind of scenarios where gauge field configurations have nontrivial topologies, the QCD Lagrangian density generally contains an extra θ-term, where g is the coupling for the strong interaction and F andF are the gluonic field strength tensor and its dual respectively. The real parameter θ defines the choice of vacuum from infinite possibilities and its value subsequently dictates whether the corresponding theory is CP symmetric or not. As can be seen straightforwardly from Eq. (2.1), only for θ = 0 (mod π), the QCD Lagrangian is CP conserving. Again various experimental studies on pseudoscalar mass ratios [38], electrical dipole moments [39][40][41][42] as well as Lattice QCD calculations [43,44] conclude that the value of this real angular parameter θ is very close to 0 in nature. A simple and elegant way to explain why θ should be so small or null is giving θ a dynamical character, elevating it to a field, the axion, such as to have a vanishing vacuum expectation value [1]. The axion field a is the canonically normalized dynamical θ, θ(x) = a(x)/f a , where f a is the axion decay constant. Its only non-derivative coupling is to the QCD topological charge and it is suppressed by the scale f a . The interaction Lagrangian density in Eq. (2.1) can now be written as 2) can be effectively represented as an interaction of the QCD axion field a with the quarks by performing a chiral rotation [45] of the quark fields by the angle a/f a , which yields [29,30] where ψ L and ψ R are the left-and right-handed components of the quark wave function ψ and G 2 is a coupling constant.
B. Thermodynamic potential for the axion background field within the NJL model The effective Lagrangian density for the isospin symmetric two-flavor NJL model with the CP violating term is given by where ψ depicts the fermionic (quark) fields and m 0 is the current quark mass. The fermionic interaction part of the Lagrangian density is given by where τ k , k = 0, 1, 2, 3, represents the unit matrix (for k = 0) and the Pauli matrices (for k = 1, 2, 3) and G 1 is the coupling constant associated with the fermionic interaction. Lastly, the symmetry breaking 't Hooft determinant interaction term L a is given by Eq. (2.3). The axion contribution Eq. (2.3) effectively breaks down the global U(2) V symmetry into SU(2) V × U(1) B . Often, in recent studies involving flavor mixing within the NJL model and its extensions, the couplings G 1 and G 2 are also taken to be equal [24,45]. In the present work we choose G 1 = (1 − c)G s and G 2 = c G s , following Refs. [24,25,46], which assumes the connection of both the coupling constants with G s , the standard scalar channel coupling constant for the NJL model. The parameter c connecting the two couplings determines the strength of the axion interaction and we will discuss more about its value in Section II D.
Next, for the derivation of the thermodynamic potential, we use the usual mean-field approximation [22], where the scalar and pseudoscalar fields are replaced by their corresponding mean-field values, or condensates, with σ and η representing the chiral and pseudoscalar condensates, respectively. Note that in the following we will be considering only the isospin symmetric case, hence, only ψ ψ and ψ iγ 5 ψ survive from Eq. (2.5). Usually, for θ = 0, or in the present case, for a = 0, ψ iγ 5 ψ also vanishes. Considering a nonvanishing ψ iγ 5 ψ emphasizes explicitly the fact that the axion field couples to the axial current. Hence, in terms of σ and η, the thermodynamic potential for the QCD axion within the NJL model, following for instance Ref. [46], is given by where the quark contribution Ω q is given by [22,45] where N c = 3 is the number of colors, E p = p 2 + M 2 and and α 0 and β 0 in the effective mass M , Eq. (2.10), can be written in terms of the axion background field a and the condensates σ and η as [46] α 0 = −2 G 1 + G 2 cos a f a σ + 2G 2 η sin a f a , (2.11) The momentum integral in Eq. (2.9) consists of a vacuum part (the first term inside the argument of the integral) and a medium part (the second term inside the argument of the integral). The vacuum contribution is ultraviolet (UV) divergent and usually taken care of by a sharp cut-off regularization procedure, i.e., with a finite three-momentum upper cut-off Λ. It can be checked in a straightforward way that the vanishing axion field limit, i.e., setting a → 0 in Eq. (2.8), reproduces the usual NJL thermodynamic potential [22]. From the thermodynamic potential Ω, Eq. (2.8), we can now find the physical values for the condensates σ and η by solving the appropriate gap equations, (2.13) which also depend on the value of the axion background field through the ratio a/f a . The effective thermodynamic potential for the QCD axion in a hot medium and within NJL model is then given by (2.14) Since in the present study the axion is treated as a background field, the axion mass m a is simply defined by the second derivative of the effective potential at vanishing axion field, i.e., where χ top is the topological susceptibility, which is independent of the scale f a , as it is evident from Eq. (2.15). Similarly, the axion self-coupling λ a is defined as the fourth derivative of the effective potential at vanishing axion field limit, . (2.16)

C. Adding an external magnetic field
We now turn to the modification of the thermodynamic potential when in the presence of an external magnetic field. The effective Lagrangian density of the QCD axion in the presence of an external electromagnetic (EM) field, within the two-flavor NJL model and at leading order in 1/f a , can be written as where D µ = ∂ µ − iqA µ , q being the electric charge and A µ the EM gauge field. F µν is the field strength tensor given by F µν = ∂ µ A ν − ∂ ν A µ andF µν is its dual. The axion-photon coupling constant g γa appearing in the interaction term is given by [47] g γa = α em 2π where α em is the EM fine-structure constant and E/N is the EM to color anomaly ratio (for example, with the value of 8/3 in light of the Grand Unification Theory (GUT) SU (5) model [47]). As mentioned in Sec. I, for our present study we consider an external anisotropic magnetic field along the zdirection. In this case, the axion-photon coupling term vanishes (F µνF µν ∝ E.B → 0) and the effective Lagrangian density simplifies to In the presence of the anisotropic external magnetic field along the z-direction, the transverse plane in the momentum space gets quantized and the dispersion relation for the quarks modifies to [48] where n and s represent the Landau levels and the spin states, respectively, and q f is the absolute charge of the fermion with flavor 1 f ≡ u d. Using 2l = (2n + 1 − s), the above relation can also be written as l being the redefined index for the Landau levels. As we shall see, this reorganization of variables produces the degeneracy factor (2−δ l,0 ) in the expression, counting the spin states for all except the lowest Landau level. Thus, by incorporating the quantized transverse momenta, the three momentum integral becomes The quark part of the thermodynamic potential Ω q now gets modified in the presence of the magnetic field and becomes and the total thermodynamic potential for the present system can be written as, It is important to note here that in the presence of the external magnetic field, the condensates corresponding to the u and d quarks are not the same anymore due to the presence of the factor q f in Ω q (B, T ). Hence, in Eq. (2.24), the new σ and η condensates are basically the average with respect to the quark flavors, i.e., σ = (σ u + σ d )/2 and η = (η u + η d )/2. For convenience, the thermodynamic potential can also be rearranged in three separate parts using the Magnetic Field Independent Regularization (MFIR) procedure. MFIR was proposed in [48][49][50] and has been recently applied in several works [51][52][53][54][55][56][57][58][59][60][61][62][63]. Using this procedure the thermodynamic potential becomes where Ω V corresponds to the vacuum part, Ω B (B) a term that carries the dependence on B only and Ω M (T, B) is the medium part which is dependent on both T and B. Each one of these terms are given, respectively, as and Λ is the finite threemomentum cut-off introduced in Section II B. The first derivative of the Hurwitz zeta function, ζ (−1, x f ), can be written in a simplified form that helps to evaluate the derivatives numerically without hassle. It is done by differentiating and integrating the function with respect to x f , leading to Hence, we get for Ω B the result where ψ (m) (x f ) represents the m-th polygamma function and the x f independent term ζ (−1, 0) has been neglected. Subsequently, using Eq. (2.25) as the thermodynamic potential for the case of axions in a hot and magnetized medium within the NJL model, we can solve the gap equations (2.13) to get the T and B dependent condensates vis-a-vis effective potential and, hence, obtain the T and B dependent axion mass, the axion self-coupling and the topological susceptibility from Eqs. (2.15) and (2.16).

D. Parametrization
Before getting the results from the thermodynamic potential, we need to discuss the choice of parametrization  [32,33]. This counterintuitive behaviour was dubbed as IMC. The dominance of sea contributions over the valence contributions of the condensate around T c is one of the probable reasons [68] behind the IMC. After the recognition of the presence of IMC through lattice QCD for both chiral and deconfinement transitions [68,69], several attempts have been made to understand it through different effective QCD models [53,59,63,[70][71][72][73][74][75][76][77]. At this point we want to emphasize that in the present work we will be discussing about both the effects of MC and IMC and the way of incorporating them in the axion thermodynamic potential, which is one of the novel features related to the present work.
In the present study, to observe the parametrization dependence of the various quantities evaluated in the vanishing magnetic field limit, we have used three different parameter sets, given in Table I. By using different parameter sets, it will enable us to quantify the dependence of the results on them.
For MC we have considered Set I from Refs. [53,59] with a fixed value of the coupling constant G s . To in- corporate IMC in the model we have also chosen the procedure taken in Ref. [59], i.e., by keeping the other parameters same as in Set I and using a four-fermion scalar coupling constant that is dependent on both the temperature and magnetic field. This was proposed for the case of the SU (2) NJL model in Ref. [59], such that where the aforementioned parameters b, β, T a and s (see Table II) are obtained by fitting the lattice data for the average of the quark condensates [33]. At this point we would also like to mention that scarcity of lattice data at lower values of temperature eventually requires this procedure to be performed with the fitting at T > 110 MeV and to extrapolate through the region of lower temperatures, which can give rise to ambiguities over the value of the coupling constant G s (eB, T = 0). This dilemma can be avoided all together following, e.g., the procedure explained in Ref. [58], where the authors have generated a coupling constant G 0 (eB) for T=0, as a good fit to the lattice simulations using selected values of eB from 0 to 1 GeV 2 , given as with values of α, β and γ given, respectively, as 1.44373 GeV −2 , 3.06 GeV −2 and 1.31 GeV −4 . Using this magnetic field dependent coupling G 0 (eB) in our present study 2 , we will separately show the behavior of the magnetic field dependence of the topological susceptibility χ top at T = 0. Finally, for the choice of the parameter c introduced in subsection II B that controls the strength of axion interaction, there have been many discussions in the past. Comparing with the three-flavor NJL model parameters [24] one gets the value of c to be around 0.2. On the other hand, successful models like the instanton liquid model, suggests c ∼ 1, thereby indicating L a as the dominant part of the Lagrangian density. In the present work most of our results consider c = 0.2 following the Refs. [25,46], though we have also shown a case comparison for the topological susceptibility with two different values of c, i.e., c = 0.2 and c = 0.8, which can be considered as representative cases when comparing the results. Furthermore, choosing the vacuum expectation value (VEV) of the scaled axion field a/f a = π as an illustration, we have done a systematic analysis of the minimum value of c (which we call as c min ), i.e., the value of c after which the condensates σ and η becomes nonvanishing, hence signalling the chiral and spontaneous CP symmetry violating phases respectively. In Fig. 1, we have shown the T − c min phase diagram for different values of external magnetic field and with different parametrizations. From this phase diagram we can clearly identify the three different phases appearing in the T − c min plane due to second order CP phase transition and the chiral crossover. So the bottom right part of Fig. 1, i.e., σ = 0, η = 0 represents both the chiral and the spontaneous CP symmetry violating phase, the bottom left part i.e., σ = 0, η = 0 stands for the phase where chiral symmetry is still violated but spontaneous CP symmetry is restored and finally the top part i.e., σ ≈ 0 shows the almost chiral symmetric phase (chiral symmetry is not fully restored due to the nonvanishing quark masses). Figure 1 also shows that for each case corresponding to different values of eB and parametrization, after certain values of higher temperatures (dependent on eB and parametrization) the spontaneous CP violating phase vanishes. These behaviors will be more transparent through the discussions of the next section. The results shown in Fig. 1 are also consistent with the ones shown for example, in Ref [25], for the values of current quark masses considered in our present work (i.e., given by Sets I, II and III). Finally all of the cases shown in Fig. 1 suggest that the chosen value of c = 0.2 is well within the viable regime.

III. RESULTS AND DISCUSSIONS
In this section we will present our results for the different quantities associated with the thermodynamics of the QCD axion in a hot and magnetized medium within the NJL model. Firstly, we will revisit the case with vanishing magnetic field [46] and discuss some relevant points related to these results. Next, we will move into the scenario with an external anisotropic magnetic field. Then, we will deal with the two different procedures, both related to the NJL model. First we will be working with a fixed coupling constant G s , which accounts for only MC. Then, finally, we will also incorporate the IMC effect in our results using the effective B and T dependent coupling constant G s (eB, T ) defined in the previous section, and discuss about the effects of both MC and IMC on the axionic QCD system.

A. The eB = 0 case
The case of a vanishing magnetic field has recently been studied with a different parametrization in Ref. [46]. In the current study we will present some different aspects along with some of the quantities already evaluated in Ref. [46], for the purposes of comparing between different parametrizations. This will be particularly useful as a way of quantifying the differences that different parametrizations can make on the results.
In Fig. 2 we have shown the variation of the physical condensates σ and η with the axion field a scaled with f a for both T = 0 and T = 180 MeV. The sinusoidal behavior can be attributed to the sin(a/f a ) and cos(a/f a ) functions present in the effective mass M , appearing explicitly in the defined quantities in Eqs. (2.10)-(2.12). It is evident from Fig. 2 that the behaviors of σ and η differ whether the results are obtained for T = 0 or for T = 0. As one can note from Fig. 2(b), at T = 0 η displays discontinuities for a/f a = (2j + 1)π, for j = 0, 1, 2, . . ., whereas for a/f a = 2jπ it vanishes. This result agrees with the Dashen's phenomena [78], which predicts the existence of two degenerate vacua at T = 0 and a/f a = (2j + 1)π due to spontaneous CP symmetry violation 3 . Different signs of the two vacua indicates that they differ by a CP transformation between them. But Dashen's phenomena starts to break down with the increase in temperature. After a certain critical temperature T c , the spontaneous CP symmetry is restored and from the CP violating phase (η = 0, σ = 0), it returns to the phase of the ordinary chiral condensate (η = 0, σ = 0). This value of T c depends on the external magnetic field and parametrization. In fact, it can be realized from the phase diagram showed in Fig. 1 also, e.g., the solid curve in Fig. 1(a) confirms the fact that for eB = 0 and using the parameters from the Set I case, the spontaneous CP violating phase has already been disappeared before T = 180 MeV. This is evident from Fig. 2(b), looking at the dashed curve at T = 180 MeV, where now we have only one vacuum. On the other hand σ attains its minimum value for a/f a = 2jπ and reaches the maximum value for a/f a = (2j + 1)π for both T = 0 and T = 0. In Fig. 3 we have shown the comparison between three different sets of parametrizations listed in Table I for the variations of the temperature dependent axion mass m a (T ) and the temperature dependent axion selfcoupling λ a (T ), scaled with their respective zero temperature values. One notices that after the value of T 0.1 GeV, the different parametrizations start to deviate from each other. For both the case of the axion mass and the axion self-coupling, a relatively rapid decrease is noticed around the chiral (pseudo critical) transition temperature T c , while for the self-coupling a kink-like feature appears in this region.
In Fig. 4 we show the results for the topological susceptibility χ top (T ), for which we also have available the lattice QCD results [80] at zero magnetic field. We have displayed the variation of χ top with respect to the temperature for the three different parametrizations sets previously mentioned and compared them both with the available lattice QCD results. As one can see from Fig. 4(a), the parameter sets I and III given in the Table I look more compatible with the lattice data. In this case, we can also see that χ top (T ) has decreased substantially around T c . This is a behavior matched by the corresponding lattice QCD data. In Fig. 4(b), we have plotted χ top (T ) in the case of choosing two different values of the constant c, corresponding to two different strengths of the axion interaction. When the strength of the axion interaction is weaker, like in the case where c = 0.2, the topological susceptibility is dragged down a little bit more at higher T than in the higher strength case of c = 0.8.
B. The case of finite magnetic field, eB = 0, but fixed coupling constant Gs In this subsection we work with the same parametrization used for the eB = 0 case, i.e., with a fixed coupling constant G s = 5.022 GeV −2 . As discussed earlier, with the fixed coupling constant we only consider the effects of magnetic catalysis in this case. Below we will also con-  Table I. sider the more physical scenario, which accounts for the inverse magnetic catalysis effect that appears around the pseudo critical temperature.
In Fig. 5, the variation of the scaled axionic part of the effective thermodynamic potential is shown with the scaled axion background field a/f a for three different values of the external magnetic field, namely, eB = 0, 0.2 GeV 2 and 0.4 GeV 2 , respectively. The aforementioned variation is shown for two different temperatures in the two panels shown in Fig. 5. It is apparent from these results that for both T = 0 and T = 180 MeV, with increasing magnetic field the value of the effective thermodynamic potential gets increased, which effectively also follows a valley-hill like structure. This behavior can be traced again due to the presence of the terms sin(a/f a ) and cos(a/f a ) in the effective mass M . It is also noticeable that with increasing temperature, the thermal effects start to melt down the amplitude of the axionic  Table I). (b) Same comparison done for the axion self-coupling ratio λa(T )/λa. part of the thermodynamic potential, which is visible most prominently when we compare the case of eB = 0 in Fig. 5 (a) and Fig. 5 (b) .
In Fig. 6 we show the effects of the temperature and magnetic field on the effective dynamical mass M (T, B), which in turn depends on the physical condensates σ(T, B) and η(T, B). The effect of MC is clearly evident from the results displayed in the figure: By increasing the magnetic field magnitude eB, the effective mass increases throughout the temperature range T = 0 − 250 MeV. The effect of the axion field is also noticeable comparing the two plots in Fig. 6 as we see that the effective dynamical mass M (T, B) gets slightly decreased with a finite value of the scaled axion field.
In Fig. 7 it is shown the temperature behavior for the axion mass and for the self-coupling, both scaled with  Table I) along with the lattice QCD results from [80]. (b) Comparison of χtop with respect to T with two different values for the constant c, i.e., two different strengths of the axion interactions along with the lattice QCD results. In this case it is used the parameters from Set I shown in Table I. their respective zero temperature values, for three different values for the external magnetic field. For both the cases shown in the figure, the magnetic catalysis effects are visible around the inflection point close to T c . We can see the inflection points shifting towards higher values of T for higher values of eB, a clear sign of MC. Similar behavior of MC is obtained in Fig. 7(c), where it is shown the topological susceptibility as a function of the temperature. The MC effect is visible throughout the temperature range. We can notice that the value of χ  Table I for the magnetized medium and for two different values of temperature, i.e., for T = 0 (a) and for T = 180 MeV (b).
C. The case of considering a temperature and magnetic field dependent coupling Gs(eB, T ) constrained by lattice QCD Finally, let us here explore the more physical scenario involving the effect of inverse magnetic catalysis around the transition temperature, originally predicted by lattice QCD for a hot and magnetized medium. As mentioned previously, for the purposes of incorporating IMC in the ambit of the NJL model, in this subsection we have used a temperature and magnetic field dependent coupling constant G s (eB, T ).
In Fig. 8 we have again plotted the scaled axion dependent part of the thermodynamic potential varying with  Table I and for two different values of the scaled axion field, i.e., for a/fa = 0 and for a/fa = 2π/3. respect to the scaled axion field. This time, using the fully temperature and magnetic field dependent coupling constant G s (eB, T ), given by Eq. (2.30). The main qualitative and quantitative difference in this case with respect to that shown in the previous Fig. 5 can be noticed in Fig. 8, where the effect of IMC is evident throughout the whole range of the axion background field. For the case of T = 180 MeV, i.e., T ∼ T c , we can see that with increasing magnetic field the amplitude of the thermodynamic potential gets decreased. For this case, the overall value of the scaled thermodynamic potential also becomes very low because of the low value of G s (eB, T ) at higher temperatures. The strong dependence on the way the coupling constant is parametrized is evident comparing Fig. 5 with Fig. 8.
In Fig. 9 we show the variation of the effective constituent quark mass M with respect to T and using lower temperatures it shows the catalysis effect as higher eB values correspond to higher effective mass. But at the higher temperatures, it goes through a crossover and around T c the behavior is completely opposite, emphasizing the occurrence of the inverse magnetic catalysis in that regime.
In Fig. 10 we have shown the variation of scaled axion mass, self-coupling and topological susceptibility with T for the fully magnetic field and temperature dependent coupling constant G s (eB, T ). The MC effect at lower temperature is not visible in Figs. 10(a) and 10(b) because of the scaling used, but for both these cases the IMC is clearly noticeable around T c : The inflection points, indicating the pseudo critical temperature location, is shifted towards lower temperature with increasing magnetic fields. The kink-like feature in the position of the pseudo critical temperature location appearing in the axion self-coupling plot shows a similar  behavior. In Fig. 10(c) the topological susceptibility is also shown in view of the fully B and T dependent coupling constant G s (eB, T ), which again shows both the MC and the IMC effects in different temperature regimes with a crossover happening in between. The quantitative difference between the zero temperature values of χ top shown in Fig. 7(c) and Fig. 10(c) is solely due to the difference in the coupling constant G s parametrization, i.e., the difference between the extrapolation of lattice fitted G s (eB, T ) T,eB=0 = 4.6311 GeV −2 with the fixed G s = 5.022 GeV −2 case shown in Fig. 7(c). Finally, as mentioned in Sec. II D, at T = 0 we can now predict the behavior of the magnetic field dependence of the topological susceptibility χ top using the magnetic field dependent coupling G 0 (eB). This is explicitly shown in Fig. 11, where the variation for both G 0 (eB) and when taking a fixed coupling constant, putting B = 0 in Eq.(2.31), i.e., G 0 (eB = 0) = 4.50373 GeV −2 , are considered.
As a final remark here, we should point out that as a central theme of our present work, let us discuss about the effect of the spontaneous CP violation observed in our results. In Fig. 12, we have shown the variations for both the σ and η condensates with respect to the temperature and for three different values of the exter- nal magnetic field. As a consequence of our discussions and results up to this point, we have only considered the running coupling G s (eB, T ) when getting these results. Also, for illustration purposes, we have chosen two specific values for a/f a , i.e., 0 and π, which we believe should suffice for our discussion. The result shown in Fig. 12(a) illustrates that there is no CP violation happening at a/f a = 0, irrespective of the increasing temperature and external magnetic field. This result is compatible with the well known Vafa-Witten theorem [79]. Now, the results displayed in Fig. 12(b) clearly show the CP violating η = 0 phase for a/f a = π in the lower temperature region. The external magnetic field dependence of the critical temperature T c (eB) is also quite visible from Fig. 12(b), which gradually decreases with increasing eB (also in agreement with the IMC phenomena). This result helps to show both the spontaneous CP symmetry breaking at T = 0, i.e., Dashen's phenomena, and the gradual restoration of the same at T = T c (eB). Then, in Figs. 12(c) and 12(d), they show similar variations for the chiral condensate. They mimic the nature of the constituent quark mass shown in Fig. 9(a). As a summary of the possible phases involved in our work, we can say that for a/f a = 2jπ we have only the chiral symmetry breaking phase, i.e., σ = 0 and η = 0, which becomes restored at a temperature T > 0.2 GeV, i.e., leading to σ = η = 0. But for a/f a = (2j + 1)π, apart from the chiral phase transition (a crossover), we have another second order phase transition at a magnetic field dependent critical temperature T c (eB), when the spontaneous CP violating phase gets restored. For a/f a = π, these three phases have explicitly been shown in the T − c min plane in Fig. 1.

IV. CONCLUSION
In conclusion, we would like to convey that in the present work we have studied, for the first time as far as it is of our knowledge, the thermodynamics of QCD axions within the NJL model in a hot and magnetized medium by explicitly emphasizing the effects of two of the most appealing phenomena in this context, i.e., the magnetic catalysis and the inverse magnetic catalysis effects. In the process, we have studied the different behaviors shown by the condensates for two different values of the scaled axion field, i.e. a/f a = 0 and π, consistent with the well known Vafa-Witten theorem and Dashen's phenomena respectively, leading to a discussion about the occurrence of different possible phases. Furthermore we have investigated the effect of the temperature and the external magnetic field explicitly on the measure of the spontaneous CP violation, thereby showing the magnetic field dependence of the critical temperature for the spontaneous CP symmetry restoration. Choosing the VEV of the axion field as a CP violating term in the QCD action our results extend the previous results in the context of strong CP violation [24][25][26][27][28] by incorporating Inverse Magnetic Catalysis effect in the medium using thermomagnetic field dependent coupling constant G s (eB, T ). On top of that, the axion mass, self-coupling and topological susceptibility have their own significance regarding the study of cold dark matter, axion stars, the cooling anomaly problem in astrophysics, etc, and we have explored the effects of temperature and external magnetic field on these quantities. Finally throughout this present work, we have also explored the importance of evaluating the differences in the results when different parametrizations are considered in the context of an effective model for QCD. In this work, we have explicitly made this study in the context of the two-flavor NJL model applied to the QCD axion thermodynamics. These results thereby strengthen the claim of a strong dependence on parametrization within the NJL model.