Fate of the first-order chiral phase transition in QCD: Implications for dark QCD studied via a Nambu-Jona-Lasinio model

The first-order nature of the chiral phase transition in QCD-like theories can play crucial roles to address a dark side of the Universe, where the created out-of equilibrium is essential to serve as cosmological and astrophysical probes such as gravitational wave productions, which have extensively been explored. This interdisciplinary physics is built based on a widely-accepted conjecture that the thermal chiral phase transition in QCD-like theories with massless (light) three flavors is of first order. We find that such a first order feature may not hold, when ordinary or dark quarks are externally coupled to a weak enough background field of photon or dark photon (which we collectively call a ``magnetic"field). We assume that a weak ``magnetic"background field could be originated from some ``magnetogenesis"in the early Universe. We work on a Nambu-Jona-Lasinio model which can describe the chiral phase transition in a wide class of QCD-like theories. We show that in the case with massless (light) three flavors, the first-order feature goes away when $2 f_\pi^2 \lesssim eB ( \ll (4 \pi f_\pi)^2)$, where $eB$ is the ``magnetic"field strength and $f_\pi$ the pion decay constant at the vacuum. This disappearance is the generic consequence of the presence of the ``magnetically"induced scale anomaly and the ``magnetic"catalysis for the chiral symmetry breaking, and would impact or constrain modeling dark QCD coupled to an external ``magnetic"field.


Introduction
Pursuing the nature of the chiral phase transition in thermal QCD is crucial to deeply comprehend the origin of matter mass in light of the thermal history of the Universe.This research direction has extensively been explored since the pioneer work by Pisarski and Wilczek [1], where the analysis was done in the chiral limit for the lightest three or two flavors based on the renormalization group analysis by employing a chiral effective model for low-energy QCD.The remarkable extension was made by the Columbia university group utilizing lattice QCD calculations [2], which has been refined along with development in modern computational technologies.In lattice QCD, the QCD phase transition including the chiral and confinement-deconfinement phase transitions have been investigated taking into the quark mass and the number of flavor dependence.It is summarized as the socalled Columbia plot, which is drawn on the quark mass plane with the classification by the universality classes in terms of the renormalization group.
Currently, the interest has been extended from the original motivation of deeply understanding the QCD phase transition to applications to Beyond the Standard Models, which has opened a vaster market involving astrophysical and cosmological fields.Of late fashion, a dark side of the Universe is discussed by assuming what is called dark or hidden QCD which can leave various cosmological and astrophysical footprints in the present-day Universe, such as gravitational wave spectra, dark MACHOs (Massive Astrophysical Compact Halo Objects) including dark compact stars, and dark matter abundances.Then, QCD with light or massless three flavors is often referred to as a typical and promising candidate which can realize the first-order chiral phase transition to supply such footprints [3? -12].
This trend has been inspired by the pioneer work [1] showing the possibility of the first-order chiral phase transition for the number of flavors greater than two, though it is based on a perturbative analysis within the framework of a chiral effective model (linear sigma model).Actually, no clear lattice results on the full QCD simulations for the firstorder phase-transition has been obtained, and instead, only lower limits on the quark mass in the crossover regime has been placed [13,14].Thus the first-order nature has not been conclusive yet.
In the scenario of dark QCD, allowing dark quarks to couple with a dark photon would make it possible to address the dark pion as a dark matter candidate and produce its cosmological abundance due to introduction of the kinetic term mixing between the photon and dark photon, as discussed in the context of SIMPs (Strongly Interacting Massive Particles) [38][39][40].Therefore, the background "magnetic" field would be favored and naturally present in dark QCD models in a sense of the phenomenological and cosmological interests, as well as ordinary QCD.
The possible presence of a weak background "magnetic" field has naively been disregarded in addressing the order of the chiral phase transition.However, it turns out to be too naive.
In this paper, we find that when ordinary or dark quarks are externally coupled to a weak enough background field of photon or dark photon (which we hereafter call a "magnetic" field collectively), thermal QCD with light three flavors may not undergo the firstorder chiral phase transition.
To make an explicit demonstration, we work on an Nambu-Jona-Lasinio (NJL) model which can describe the chiral phase transition in a wide class of QCD-like theories including dark QCD.We include the contribution from the scale anomaly induced by quark couplings to a "magnetic" field, as discussed in the literature [41], which is present only in a weak field regime.It then turns out that the disappearance of the first-order nature in light three-flavor QCD arises as the generic consequence of the crucial presence of the "magnetically" -induced scale anomaly and the "magnetic" catalysis [42][43][44] for the chiral symmetry breaking at around the criticality, which is irrespective to details of the parameter setting in the NJL model description monitoring ordinary QCD.
The size of the critical "magnetic" field strength (eB cr ) normalized to f π is estimated as eB cr /f 2 π ∼ 2, above which (but still at weak eB (4πf π ) 2 ) the first-order feature goes away, where f π is the pion decay constant at the vacuum with T = eB = 0. Further allowing quarks to have finite masses, we show that the first-order feature in the light quark mass regime is completely gone at the critical field strength, which is captured in a couple of extended Columbia plots including the eB-axis.
Our finding not merely helps deeply understanding the chiral phase structure in QCD, but would thus impact modeling dark QCD in the NJL description and testing the model via the astrophysical probes and constrain the "magnetogenesis".

NJL model with electromagnetic scale anomaly
We begin by briefly introducing an NJL model with the light three flavors as a low-energy effective model of QCD, and present the Lagrangian and the thermodynamic potential with a constant magnetic field (which we shall hereafter call the thermomagnetic potential) containing the contribution of the electromagnetic scale anomaly.
Throughout the present paper, we refer to the ordinary up, down, and strange quarks with the electromagnetic charges 2/3, −1/3, and −1/3 in unit of the electric charge e and adopt ordinary QCD as the monitor for a wide class of QCD-like theories.However, as will be clarified later, the qualitative features associated with the chiral phase transition are substantially intact even if the charge assignments and model-parameter sets are taken to be adjusted to some dark QCD setup with photon or dark photon, differently from the ordinary QCD case.

Model description and parameter setting monitoring ordinary QCD
The three-flavor NJL model externally coupled to a background photon field is described by the following Lagrangian 1 : where ψ = (u, d, s) T represents the Dirac fermion field forming the flavor SU(3) triplet of up, down, and strange quarks (with QCD color and Dirac spinor indices suppressed); m denotes the quark mass matrix m = diag{m 0 , m 0 , m s }; λ a (a = 0, • • • , 8) are the Gell-Mann matrices in the flavor space, and λ 0 = 2/3 • 1 3×3 .In Eq.(2.1) the covariant derivative the coupling term between the quark field and the external photon field A µ along with the electromagnetic charges in Q em = diag{q u , q d , q s }, where q u = 2/3, q d = −1/3, and q s = −1/3.For simplicity, we introduce a constant magnetic field along the z-direction in space to be embedded in the external photon field as A µ = (0, 0, Bx, 0).The four-fermion interaction with the coupling G is invariant under the full chiral U (3) L ×U (3) R symmetry with the chiral transformation: ψ → U •ψ with U = e −iγ 5 8   a=0 (λ a /2)θ a 1 For reviews on the three-flavor NJL model, see, e.g., Refs.[45,46].and the chiral phases θ a .This four-fermion operator is supposed to be induced from the one-gluon exchange process at the confinement scale.
The Kobayashi-Maskawa-'t Hooft (KMT) determinant term [47][48][49][50] with the coupling K is a six-point interaction induced from the QCD instanton configuration coupled to quarks.The KMT term preserves SU (3) L × SU (3) R invariance (associated with the chiral phases labeled as a = 1, • • • , 8), but breaks the U (1) A (corresponding to a = 0) symmetry.This interaction gives rise to the mixing between different flavors and also uplifts the η mass to make η no longer a Nambu-Goldstone boson.
The quark mass terms with m f (f = u, d, s) and the couplings to the background photon field with the charges eq f are flavor dependent.Therefore, those explicitly, but weakly enough, break the chiral U (3) L × U (3) R symmetry as long as the strength of eB is perturbatively small, say, eB O(f 2 π ), which is indeed the case in the present setup.Since the U (1) A symmetry is to be badly broken by the KMT term so as to reproduce the observed η mass, the model keeps the approximate chiral SU (3) L × SU (3) R symmetry only.
When the couplings G and/or K get strong enough, the approximate chiral SU (3) L × SU (3) R symmetry is spontaneously broken by nonperturbatively developed nonzero quark condensates ψf ψ f = 0, down to the vectorial symmetry SU (3) V (corresponding to invairance associated with the vectorial phase rotation ψ → e −i 8 a=0 (λ a /2)θ a • ψ), to be consistent with underlying QCD.The present NJL model monitors the spontaneous breakdown by the large N c expansion or equivalently in the mean field approximation, where N c stands for the number of QCD colors, to be set to 3.
The NJL model itself is nonrenormalizable because the G and K coupling terms are higher dimensional operators greater than four in mass dimension.Therefore, a momentum cutoff Λ needs to be introduced to make the NJL model regularized.
Relevant NJL formulas to compute QCD hadronic observables at the vacuum are presented in Appendix A. We quote the literature [46] where the two coupling parameters G and K normalized to a three-dimensional momentum cutoff Λ are fixed as with Λ to the pion decay constant f π : This parameter set has been determined by fixing four representative hadronic observables in the isospin-symmetric limit at T = eB = 0, such as the pion decay constant f π = 92.4 MeV, the pion mass m π 135.0 MeV, the kaon mass m K 497.7 MeV, and the η mass m η 957.8 MeV, with m 0 = 5.5 MeV input, which also yields m s = 140.7 MeV, Λ = 602.3MeV, and the η mass m η 514.8 MeV.Those isospin-symmetric hadronic quantities are in good agreement with the recent lattice results at the physical point of 2 + 1 flavor QCD with m π = 138 MeV [51]: m K = (494.2± 0.4) MeV [51], m η = (554.7 ± 9.2) MeV [52], and m η = (930 ± 21) MeV [52].
When working at finite temperature, we will not consider intrinsic-temperature dependent couplings, instead, all the T dependence should be induced only from the thermal quark loop corrections to the couplings defined and introduced at the vacuum with T = eB = 0. Actually, the present NJL at eB = 0 qualitatively fits well with lattice QCD results at the physical point on the temperature dependence for the chiral, axial, and topological susceptibilities, as shown in Ref. [53] 2 .In this sense, we do not need to introduce such an additional T dependence into the model parameters.
The QCD-monitored parameter setup in Eqs.(2.2) and (2.3) forms a wide class of QCDlike theories in the NJL description obtained by scaling up of ordinary QCD with respect to Λ or f π .A similar parameter setup derived based on a different regularization scheme has been applied in dark QCD scenarios with use of the scaling up [3-5, 7, 9? ].

Thermomagnetic potential with electromagnetic scale anomaly
We work in the mean field approximation, so that from Eq.(2.1) the constituent quark masses (i.e. the full quark masses) are read off as [46] M where σ f s correspond to the dynamical quark masses, and S f (p) denotes the f -quark full propagator in the mean field approximation, which, in the Landau level representation, takes the form [54] with p ⊥ ≡ (p 1 , p 2 ), l f ≡ 1/ |q f eB|, and Here L α n (z) are the generalized Laguerre polynomials, To evaluate thermal and magnetic contributions to the thermomagnetic potential in the mean field approximation, we apply the imaginary time formalism and the Landau level decomposition.The following replacements can then be activated 3 : Here ω k denotes the Matsubara frequency with the level of k and n represents Landau levels.

Magnetically induced scale anomaly terms
As discussed in Ref. [41], in the thermomagnetic medium with a weak enough magnetic field strength (eB (4πf π ) 2 ), the electromagnetic scale anomaly induces tadpole potential terms for the chiral singlet component of the scalar mesons.At the quark one-loop level, the induced potential takes the form ϕ denotes the chiral singlet part of the scalar mesons arising as the radial component of the nonet-quark bilinear qf q f , so that ϕ ≈ σ 2 u + σ 2 d + σ 2 s , where σ f s are diagonal elements of qf q f .f ϕ is the vacuum expectation value of ϕ at T = eB = 0. Then the prefactor ϕ/f ϕ in Eq.(2.8) is given as where σ f 0 denotes the vacuum expectation value of σ f at T = eB = 0 with σ u0 = σ d0 = σ 0 , while σ f s are the background fields depending on T and eB.
T µ µ (M f , T, eB) in Eq.(2.8) is the trace anomaly induced from the quark loop with the constituent mass M f in Eq.(2.4) in the thermomagnetic medium with a weak eB (4πf π ) 2 , which takes the form [41] ) 3 In the present study, we simply apply a conventional soft cutoff scheme [55] to regularize the integration over p3 in Eq.(2.7).In this regularization scheme, the parameter setting in Eqs.(2.2) and (2.3) together with m0 = 5.5 MeV, ms = 140.7 MeV, and fπ = 92.4MeV yields mπ 138 MeV at T = eB = 0, which differs by about 2 % from the one (mπ 135 MeV) based on a three-momentum cutoff scheme with the same inputs used [See the discussion below Eq. (2.3)].The same level of discrepancy can be observed in other meson masses.This small dependence on the regularization scheme does not make substantial difference on the qualitative trend for the chiral-phase transition nature which we will mainly discuss later on. where In Eq.(2.10) β e in the first term of Eq.(2.10) stands for a beta-function coefficient for the electromagnetic coupling e divided by e 3 .When evaluated at the one-loop level, β e is independent of the electromagnetic coupling e: The first term in Eq.(2.10) is nothing but the well-known electromagnetic scale anomaly at the magnetized vacuum with a constant magnetic field eB, while the second term corresponds to the thermomagnetically induced part [41].The function F (T, M f ) is positive definite and monotonically increases with T , hence the tadpole term in Eq.(2.8) keeps developing with respect to both T and eB.
Both of two electromagnetic-scale anomaly terms originally take the form of the transverse polarization like F 2 µν in the operator form, reflecting the electromagnetic gauge invariance.Therefore, they will tend to vanish as the magnetic field gets strong enough, because the dynamics of quarks are governed by the lowest Landau level states longitudinally polarized along the direction parallel to the magnetic field.Thus, the electromagnetic scale anomaly is most prominent in a weak magnetic regime, as noted in the literature [41].
Before proceeding, we shall give several comments on the validity of the present formulation as the thermomagnetic chiral model in light of properly investigating the chiral phase transition around the three-flavor chiral limit.
• In Eq.(2.8) the lightest-chiral singlet-ϕ meson has been assumed to be mostly composed of the quarkonium state (q f q f ).A glueball field could mix with ϕ, hence might contribute to the scale anomaly relevant to the chiral phase transition.However, lattice simulations for 2 + 1 flavors have clarified that around the two-flavor chiral limit with m 0 less than the physical-point value, the fluctuation of Polyakov loop, which can be regarded as (an electric part) of glueball, does not have significant correlation with the chiral order parameter [56].This indicates negligible mixing between the glueball and ϕ in the two-flavor chiral limit that in the present work we focus mainly on.Hence such a gluonic term will not be taken into account later on.
• Though the mixing with a glueball is irrelevant, one might still suspect a direct contribution from the gluonic scale-anomaly (of the form ∼ ϕG 2 µν ) to the chiral phase transition, which is thought to be present all the way in thermal QCD including the vacuum.It is known that the gluonic scale-anomaly can be introduced as the log potential form in a scale-chiral Lagrangian [57].The log potential form like ϕ 4 log ϕ does not give any interference for the electromagnetic scale-anomaly induced tadpoleterm which dominates around ϕ ∼ 0.
• In deriving the thermally induced tadpole part in Eq.(2.10), we have ignored the next-to-leading order terms of O((eB) 2 ) [41], which might come in the formfactor F as the nonlinear eB dependence.Such higher order corrections would naively be suppressed by a loop factor, compared to the leading order term of O(eB) in Eq.(2.10).Even if it would not sufficiently be suppressed, and would constructively enhance or destructively weaken the tadpole contribution, the presently addressed chiral-phase transition nature would not substantially be altered, as long as the magnetic field is weak enough.To be more rigorous, it is subject to nonperturbative analysis on the photon polarization function coupled to the sigma field.In Ref. [58] the authors have made an attempt to straightforwardly compute the photon polarization in the weak field approximation at finite temperature.It would be interesting to apply such a technique described in the literature to our study, which is, however, beyond the scope of the present paper, instead deserve to another work.

Total form of thermomagnetic potential
The thermomagnetic potential containing the tadpole term in Eq.(2.8) in the mean field approximation thus takes the form where with E 2 q f ,n = p 2 z + 2n |q f eB| + M 2 f being the energy dispersion relation for the f -quark, and α n = 2 − δ n,0 , which counts the spin degeneracy with respect to Landau levels.The thermomagnetic potential is thus evaluated as a function of quark condensates φ u , φ d , and φ s through Eq.(2.4) with Eq.(2.5).Given m s , T , and eB, we search for the minimum of the thermomagnetic potential in Eq.(2.12) by focusing on the chiral limit for up and down quarks, m 0 → 0, through the stationary condition, Figure 1.An illustration of schematic view of the potential deformation around the chiral phase transition at T f π in the massless three-flavor QCD, with eB = 0 (left panel) and eB = O(f 2 π ).The thermodynamic (magnetic) potential Ω (the vertical axis) is described as a function of the chiral order parameters φ u and φ d (the horizontal axis).The arrows in the panels depict the evolution of the vacuum when T gets higher and higher.In the left panel with eB = 0, the first-order transition is realized due to the sizable cubic potential arising from the KMT term, while in the right panel with eB = O(f 2 π ) the transition becomes smoother, like the second order or crossover, because the cubic term fails to keep significant against the magnetically enhanced quadratic term which destructively contributes over the thermal correction (due to the magnetic catalysis), and only the single vacuum is favored due to the presence of the 'magnetic-scale anomaly-induced tadpole term.For more details, see the text.

Disappearance of the first-order chiral phase transition
In this section, we show that the first-order feature of the chiral phase transition in the light or massless three-flavor case is wiped out in the presence of a weak magnetic field background.We start with a heuristic approach to qualitatively understand the essential difference on the type of the phase transition orders between the cases with and without a weak background magnetic field in the massless three-flavor limit [See Fig. 1].Then we confirm our observation explicitly by numerical analysis [See Fig. 2], and extend the parameter space by including finite strange quark mass m s , which is projected onto an extended Columbia plot in the (eB/f 2 π , m s /f π ) plane [See Fig. 3].

Heuristic observation: Overview of the disappearance
First of all, consider massless three-flavor QCD with eB = 0, and note the presence of the cubic potential (σ 3 ) induced from KMT determinant term, where σ stands for a symbolic order parameter regarding the chiral SU (2) L × SU (2) R symmetry, such as up and down quark condensates φ u and φ d in Eq.(2.4) as in the Ginzburg-Landau description.This cubic term will not get any correction from the thermal-quark loops, hence keeps contributing to the potential at any T with negative sign like V σ 3 = −|c| • σ 3 (with some coefficient |c|) which is definitely fixed by the positive sign of the η mass squared.This is because the thermal quark loop correction does not break the chiral SU (2) L × SU (2) R symmetry or the U (1) A axial symmetry in the mean field approximation.In contrast, the quadratic (σ 2 ) and quartic (σ 4 ) terms get sizable thermal corrections, which make the potential lifted up and the potential energy at the nontrivial vacuum smaller in magnitude as T gets higher.Thus the thermal corrections help keeping the potential barrier created by the cubic term sizable enough to stay even at high T , such that the size of the σ 3 term becomes comparable to those of the σ 2 and σ 4 terms in the potential, and then the phase transition happens at the first order with a jump from σ = 0 to σ = 0.This feature at around the critical temperature is generically seen irrespective to the sign of the original quadratic term at T = 0 and the size of the cubic term.See the left panel of Fig. 1, which would help qualitatively grasping the potential deformation around the phase transition era.
Next, turn on a weak eB, say, as weak as the strength of O(f 2 π ).Two remarkable natures, which act destructively against the thermal corrections, will then be seen: i) One is the magnetic catalysis, which makes the chiral-broken nontrivial vacuum and the vacuum-potential well deeper as eB increases against the thermal corrections (i.e.makes the chiral breaking promoted) 4 ; ii) The other is emergence of the thermomagnetic tadpole in Eq.(2.8), which makes the chiral-symmetric vacuum σ = 0 no longer the vacuum by the negative linear slope at σ = 0 [See Eq.(2.10)].
Due to the phenomenon i), the σ 2 term in the potential gets negatively and significantly enhanced.Convoluted with the second feature ii), as eB increases, the σ 3 term will then lose the chance to become comparable with the σ 2 term and fail to jump into another vacuum even at high T , thus there is only the single vacuum, and the phase transition will be smoother like the type of crossover or the second order.See the right panel of Fig. 1 5 .This is the overview of what will happen in the presence of a weak magnetic field along with the induced tadpole potential, in conjunction with the magnetic catalysis, which is robust and generic, and thus the first-order nature of the chiral phase transition in the massless three-flavor QCD disappears.

Analysis in the massless three-flavor limit
Using the QCD-monitor parameter set in Eqs.(2.2) and (2.3), we evaluate the T dependence of the constituent mass of up quark M u in Eq.(2.4) with eB varied.Note that the chiral symmetry is explicitly broken even in the massless limit, because of the magnetic field, 4 Lattice QCD has observed the inverse magnetic catalysis which acts constructively with increase of T at high enough T = O(100 MeV) [59][60][61][62][63][64].The magnetic field strength is bounded from below as eB T 2 due to the currently applied method to create the magnetic field [65], which is only for the physical pion mass simulations at finite temperatures.Thus applying a weak eB at high T with much smaller pion masses, particularly close to the chiral limit, is challenging due to the higher numerical cost, and it is still uncertain whether the normal or inverse magnetic catalysis persists. 5The magnetic catalysis in i) would also play an crucial role even when the gluonic scale anomaly contribution is incorporated in the system following the procedure as in a chiral-scale effective theory [57]: the thermal mass term ∼ T 2 ϕ 2 gets destructively interfered with a (relatively) strong eB, enough not to build a potential barrier which could be created when the logarithmic potential term ∼ ϕ 4 log ϕ associated with the gluonic scale anomaly is present.hence the chiral restoration cannot exactly be achieved even at high enough T and the criticality becomes pseudo.We therefore determine the order of the phase transition by reading off the sharpness of dM u /dT at the pseudocritical temperature T pc defined at the inflection point as d 2 M u /dT 2 | T =Tpc = 0.When a discontinuous jump in dM u /dT | T =Tpc is observed, we interpret the phase transition to be of first order, otherwise a smooth transition including the second order and crossover.
The analysis has been done fully in terms of dimensionless quantities normalized to the pion decay constant at the vacuum (T = eB = 0), f π : all relevant things are scaled by the powers of f π and can be evaluated as a function of T fπ , eB Note that the electromagnetic coupling e itself is a redundant parameter, because it always arises with B as the combination of eB as long as the gauge interaction is perturbatively small enough to be well approximated at the one-loop level as in Eq. (2.11).
Thus our results will be fairly generic and can be applied to a wide class of QCD-like theories including dark QCD coupled to dark photon with charges (q f ) of O(1) like the electromagnetic charges of the ordinary quarks.
In Fig. 2 we show the T dependence of M u , normalized to f π , in the massless threeflavor limit (m 0 = m s = 0), for a couple of reference cases with eB/f 2 π = 1.8 (left panel) and 2.3 (right panel).We have confirmed that numerical results are unchanged even raising the Landau level up to n = 500 for eB/f 2 π < 1.8 including eB = 0, which means that our calculations are accurate enough.The first-order phase transition can still be seen in the left panel, however, in the right panel with the larger eB the transition gets changed to be smoother, which is of the type of the second order or crossover, in accordance with the heuristic observation in Sec.3.1 [See Fig. 1].
The strength of the first-order phase transition can be evaluated by measuring the size of the ratio ∆M u /T pc with ∆M u = M before u (T pc ) − M after u (T pc ) being the difference of the degenerated values of M u before and after the phase transition at T = T pc .From the left panel of Fig. 2, in the case with eB/f 2 π = 1.8 (blue curve) compared to the case with eB = 0 (black curve), it reads This implies that the strength of the first-order phase transition gets reduced to be about 37% of the one at eB = 0, where the potential is normalized as Ω[M u = M after u (T pc )] = 0.This happens due to two phenomena: i) the chiral symmetric vacuum (M u = 0) is no longer realized in a magnetic field even at high T ; ii) the temperature required to make M u jumped down gets higher to further decrease M u until the transition takes place, which is due to the magnetic catalysis, so that ∆M u becomes smaller.Thus the strong first-order feature in massless-three flavor QCD (with ∆M u /T pc > 1) significantly gets altered by increasing a weak eB to be finally gone as in the right panel of Fig. 2.
The reduction of the strength of the first-order phase transition further leads to weakening of the strength of the produced latent heat L(T ), which is evaluated with fixed eB/f 2 π as Then we find Thus the gravitational wave spectrum created via the bubble nucleation associated with this dark chiral phase transition gets weaker (by about 15%).The associated generation mechanism for the dark baryon -anti-dark baryon asymmetry and the stable formation of dark baryonic compact stars and MACHOs could also be affected.

Fate of the first-order regime in a view of Columbia plot
Further allowing m s to take finite values, in Fig. 3 we make an extended Columbia plot on the (m s /f π , eB/f 2 π ) plane in the two-flavor chiral limit with m 0 = 0. We have taken the representative four points (m s /f π , eB/f 2 π ) (0.162, 0.586), (0.0758, 1.17), (0.0325, 1.76), and (0, 2.08), and interpolated those reference points to draw the phase boundary separating the first order and second order -crossover regimes.
The range of m s /f π to realize the first-order phase transition gets narrower and narrower, as eB increases.This trend confirms the conjecture made in Ref. [41]  The phase diagram in a view of an extended Columbia plot projected onto the (m s /f π , eB/f 2 π ) parameter space on the m 0 = 0 plane.The four reference points (denoted as red blobs) have been chosen at (m s /f π , eB/f 2 π ) (0.162, 0.586), (0.0758, 1.17), (0.0325, 1.76), and (0, 2.08), to make the interpolation curve for the phase boundary (corresponding to the black dots) which separates the first order and second order -crossover phases.
to be present at a finite m s on the m s axis in the extended Columbia plot, which corresponds to the "2nd order -crossover" domain in Fig. 3 including the massless two-flavor limit with extrapolation of m s → ∞ 6 .
The first-order regime completely disappears when eB reaches the critical strength, eB cr /f 2 π 2.1.This probes the disappearance of the first-order regime due to the convolution of the magnetic catalysis with the electromagnetic scale anomaly, as has heuristically been viewed in Sec.3.1.
We have also examined the critical quark mass in the three-flavor symmetric limit (m 0 = m s ), at which the phase transition feature turns from the first order to the secondorder, to be crossover.In the case with eB = 0 we have m cr 0 | eB=0 /f π 0.048.Applying a weak eB, we find that m cr 0 monotonically decreases as eB increases as plotted in Fig. 4.This is consistent with Fig. 3, which shows that the first-order domain (the range of m s /f π ) gets shrunk as the strength of the magnetic field gets stronger.Indeed we have observed m cr 0 at the critical magnetic field eB cr /f 2 π 2.1.Increasing the size of the quark mass allows even weaker eB to cease the phase transition at the first order, as seen in the phase diagram, Fig. 4.
Figures 3 and 4 thus imply that the critical magnetic field eB cr /f 2  plete disappearance of the first-order nature in the small quark mass regime present in the conventional Columbia plot without a magnetic field [2].

Conclusion and outlook
In conclusion, thermal QCD with light three flavors may not undergo the first-order chiral phase transition, when ordinary or dark quarks are externally coupled to a weak enough background field strength of photon or dark photon, which we collectively call a background "magnetic" field.The critical "magnetic" field strength (eB cr ) normalized to the pion decay constant at the vacuum with T = eB = 0, f π , is found: eB cr /f 2 π ∼ 2, above which (but at still weak eB (4πf π ) 2 ) the first-order feature goes away [See Figs. 3, 4, and also Eqs.The disappearance of the first-order regime is the generic consequence of the crucial presence of the "magnetically" -induced scale anomaly and the "magnetic" catalysis for the chiral symmetry breaking at around the criticality, which is irrespective to details of the parameter setting such as those in Eqs.(2.2) and (2.3) [See also Fig. 1].Indeed, the trend of the disappearance of the first-order feature as shown in Fig. 3 (or 2) is fairly generic and can be applied to a wide class of QCD-like theories, which might still hold even beyond the mean-field approximation in the NJL description.
The dark photon is typically the portal between the dark and the Standard Model sectors.In particular, the dark photon coupling e is crucial to drain the entropy of the Standard Model particles in the Universe.This generically constrains the value of e for dark QCD models such as SIMP models to be phenomenologically and cosmologically viable.As noted in Sec.3.2, the present analysis is irrespective to the size of e itself, because the external-magnetic field strength contributes to dark QCD models with keeping the gaugeinvariant form F 12 = eB without separation between e and B. Therefore, the presence of the critical strength eB cr /f 2 π and the extended Columbia plots in Figs. 3 and 4 do not give any constraint on a singled e, or a dark photon portal, hence are generic consequences applicable for any size of dark photon coupling e in any dark QCD scenario.
When the "magnetogenesis" is realized by some inflationary scenario along with reheating Universe, the initially produced strength of √ eB 0 could be on the order of the reheating temperature.Passing through the redshift down to the dark QCD epoch, one then considers the critical condition (eB 0 )eB cr ∼ 2f 2 π .This implies that our finding can constrain the magnetogenesis and the reheating scenario, instead of the dark photon physics in dark QCD.
The present conclusion would thus give significant impacts on modeling and testing dark QCD via the astrophysical probes and/or constraining the magnetogenesis, as well as would help deeply understanding the chiral phase structure in QCD.
In closing, we give several comments on the future prospects in order.
• The chiral extrapolation on lattice in absence of a magnetic field has systematically been studied, which has also been improved so much in the case with 2+1 flavors [66].
In addition, a remarkable progress on creating small magnetic fields at the physical point has been reported, and the current lower bound is ∼ 100 MeV [65].Calculations on derivatives with respect to magnetic fields have also been worked out directly on the lattice, so in principle small magnetic fields become accessible [67].Therefore, we expect that the disappearance of the first-order nature could be explored on lattice QCD in near future.
• The disappearance of the first-order feature also indicates null tricritical point at which all the first, second, and crossover regimes merge [2], so that the rich phase structure in the conventional Columbia plot is completely gone, and only the crossover nature survives.This phenomenon has been conjectured in the literature [41], in which it is dubbed like scale-anomaly pollution.Currently this pollution gets more quantified in a sense that it has been found to take place at eB above the critical strength, eB 2f 2 π .It would be intriguing to examine how this pollution will go away as eB increases up to the strength as much as the strong field regime (∼ (4πf π ) 2 ), and the crossover feature can keep persisting.To this end, it would be necessary to formulate the electromagnetic-scale anomaly term as a formfactor depending on the photontransfer momentum and polarization tensor structure, so that the decoupling of the Landau levels in the photon polarization function is made automatically embedded to interpolate the weak and strong magnetic field regimes continuously.Presumably, we then need to nonperturbatively compute the ϕ -photon -photon triangle diagram involving the nonperturbative quark propagators, which has never been carried out.Our present result would motivate people to work on this issue to clarify the fate of the scale-anomaly pollution as increasing eB and also existence of a new tricritical point on the eB -axis, as can be expected from Figs. 3 and 4.
• The current analysis has been done based on the NJL model in the mean field approximation.This approach is on the same level of reliability as in the literature [3-5, 7, 9? ] in a sense of employing low-energy effective models to address the chiral phase transition instead of the underlying thermal dark QCD.Since the linear sigma model belongs to the same universality class as the NJL model in terms of the renormalization group, the current conclusion on the disappearance of the first-order nature could directly be applied also to the other dark QCD scenario [6-8, 10-12, 68? , 69] argued based on the linear-sigma model description including the pioneer work by Pisarski and Wilczek [1].Though being effective enough to take the first step, the current model analysis should anyhow be improved to be more rigorous by going beyond the mean field approximation, say, the nonperturbative renormalization group method.This way we would get a more conclusive answer to whether the first-order nature of the chiral phase transition surely goes away in thermomagnetic QCD with massless three flavors.
• When contributions beyond the mean field approximation, or equivalently the subleading corrections in the 1/N c expansion are taken into account, the meson loops would come into play in the chiral phase transition.If one then considers dark QCD with many flavors, the elimination of the first-order feature by a background "magnetic" field could be avoided because a sizable thermal-cubic term amplified by the large number of flavors would be generated to build a strong enough potential barrier at around the critical temperature.This implies that the lower bound on the number of flavors might be placed with eB/f 2 π fixed in such a way that a strong enough first order (with ∆M u /T pc 1) is realized.This would be interesting to study and give further impact on astrophysical and cosmological probes from dark QCD with many flavors [6,8,10,11].The issue can be simplified by working in the linear sigma model description, as an application of studies on many-flavor QCD in the literature [70,71] with a magnetic field added, which would be worth pursuing elsewhere.
• It is of interest to clarify the full phase structure of the extended Columbia plot in three dimensions spanned by the (m 0 /f π , m s /f π , eB/f 2 π ) parameter space (such as an extension of Figs. 3 and 4), whichever the approach goes beyond the mean field approximation, or not.
the Seeds Funding of Jilin University (S.M.).The work by M.K. was supported by the Fundamental Research Funds for the Central Universities.

Figure 2 .
Figure 2. The T dependence of M u , normalized by the pion decay constant at the vacuum with T = eB = 0, f π , in the case with massless three flavors (with m 0 = m s = 0) for eB/f 2 π = 1.8 (blue curve), in comparison with the case with eB = 0 (black curve) (left panel: the first order transition), and eB/f 2 π = 2.3 (right panel: the second order transition or crossover).The close-up window in the left panel zooms-in the difference between the cases with eB = 0 and eB/f 2 π = 1.8 around the critical temperature.
Figure 3.The phase diagram in a view of an extended Columbia plot projected onto the (m s /f π , eB/f 2 π ) parameter space on the m 0 = 0 plane.The four reference points (denoted as red blobs) have been chosen at (m s /f π , eB/f 2 π ) (0.162, 0.586), (0.0758, 1.17), (0.0325, 1.76), and (0, 2.08), to make the interpolation curve for the phase boundary (corresponding to the black dots) which separates the first order and second order -crossover phases.

2 m 0
=m s plane

Figure 4 .
Figure 4. Another look at the phase diagram in terms of an extended Columbia plot, on the m 0 = m s plane (in the three-flavor symmetric limit) in the (m 0 /f π , eB/f 2 π ) space.The black curve corresponds to the critical surface separating the two phases, the first order and second ordercrossover regimes.