Anomalous low-temperature thermodynamics of QCD in strong magnetic fields

The thermodynamics of quantum chromodynamics at low temperatures and in sufficiently strong magnetic fields is governed by neutral pions. We analyze the interacting system of neutral pions and photons at zero baryon chemical potential using effective field theory. As a consequence of the axial anomaly and the external magnetic field, the pions and photons mix with one another. The resulting spectrum contains one usual, relativistic photon state, and two nonrelativistic modes, one of which is gapless and the other gapped. Furthermore, we calculate the leading, one-loop contribution to the pressure of the system. In the chiral limit, a closed analytic expression for the pressure exists, which features an unusual scaling with temperature and magnetic field, $T^3B/f_\pi$, at low temperatures, $T\ll B/f_\pi$. Finally, we determine the pion decay rate as a function of the magnetic field at the tree level. The result is affected by a competition of the anisotropic kinematics and the enlarged phase space due to the anomalous mass of the neutral pion. In the chiral limit, the decay rate scales as $B^3/f_\pi^5$.


Introduction
Extreme magnetic fields with strengths estimated to reach up to 10 19 G can exist in the universe, most notably in the terrestrial relativistic heavy ion collision experiments, and in the interiors of certain neutron stars: the magnetars. Since the energy scale associated with such magnetic fields is comparable to the characteristic scale of strong nuclear interactions, this fact has prompted intensive work on the structure of the phase diagram of quantum chromodynamics (QCD) in presence of magnetic fields (see ref. [1] for a recent review).
While most of the efforts have focused on the effect of the magnetic field on the chiral order parameter of QCD and the chiral phase transition at high temperature and zero baryon chemical potential, it was noticed early on that a nonzero chemical potential combined with the magnetic field leads to a qualitative change of the QCD ground state [2]. Namely, by virtue of the chiral anomaly, sufficiently strong magnetic fields lead to a spatially nonuniform condensate of neutral pions, which takes the form of a soliton lattice and can be energetically favored over normal nuclear matter [3,4]. A number of recent publications is devoted to a detailed investigation of this new state of matter [5][6][7].
In this paper, we draw inspiration from the above studies and investigate the effects of the chiral anomaly on the low-energy dynamics and low-temperature thermodynamics of QCD in strong magnetic fields and zero baryon chemical potential. In this case, the structure of the uniform, chiral-symmetry-breaking QCD vacuum remains unaltered, at least for magnetic fields B m 2 ρ ≈ 0.6 GeV 2 [8]. At the same time, charged pions become heavy as a consequence of the Landau level quantization. Hence, the low-energy physics of QCD in presence of a magnetic field is dominated by neutral pions and photons. This is an ideal system to test the effects of the chiral anomaly: the neutral pions do not couple minimally to the electromagnetic field, and the anomaly thus provides the only interaction between the two subsystems.
It is well known that in the vacuum, a neutral pion decays anomalously into two photons. At the effective-field-theory (EFT) level, this process can be described by an interaction term proportional to φ µναβ F µν F αβ , where φ represents the neutral pion and F µν ≡ ∂ µ A ν − ∂ ν A µ is the usual electromagnetic field strength tensor. Our key observation is that in presence of a background magnetic field, B ex , the same interaction leads to a mixing between pions and photons. This is detailed in section 3 upon a brief overview of the EFT setup in section 2. It turns out that one of the two photon polarizations is insensitive to the presence of the background field and retains its relativistic dispersion relation. It is easy to see that this mode carries electric field perpendicular to B ex , since for such field configurations, µναβ F µν F αβ ∝ E · B vanishes. The other photon polarization mixes with the neutral pion, giving rise to two modes with nonrelativistic and anisotropic dispersion relations. One of the two modes is gapless and, interestingly, its dispersion relation becomes quadratic for directions of propagation perpendicular to B ex .
The following sections then present an analysis of some more directly observable consequences of this anomaly-induced mixing between pions and photons. First, in section 4, we evaluate the pressure of the system at nonzero temperature in the leading, one-loop approximation. Intriguingly, there is a very simple closed analytic expression for the pressure in the chiral limit (zero pion mass) despite the complicated dispersion relation of the pion-photon modes. In section 5, we then analyze the effect of the magnetic field on the above-mentioned anomalous electromagnetic decay of the neutral pion. Finally, in the concluding section 6, we get back to some of the assumptions implicit to our analysis. First we discuss the separation of scales corresponding to the charged and neutral pion sectors, which defines the range of magnetic fields in which our EFT is applicable. Second, we compare the anomalous contributions to the pion spectrum at the tree level to the normal, one-loop corrections, neglected here.
Our main tool is the chiral perturbation theory [9][10][11], which governs the dynamics of QCD at low energies and temperatures. Its predictions for observables are organized in a derivative expansion, controlled by the parameter p/(4πf π ), where f π ≈ 92 MeV is the pion decay constant and p the characteristic momentum scale of the system. This limits the validity of our results to temperature and magnetic field scales well below the scale of chiral symmetry breaking. Moreover, we resort to the lowest-order approximations appropriate for the given physical observable. Thus, the excitation spectrum in section 3 and the pion decay in section 5 are analyzed at tree level, and thus include no feedback from the thermal medium. The result for the pressure presented in section 4 then corresponds to a gas of free quasiparticles with dispersion relations fixed to their zero-temperature values.
Throughout the paper, we use the natural units in which the Planck constant , speed of light c as well as the elementary electric charge e are all set to one. The magnetic field strength is given in the high-energy-physics units of 1 GeV 2 ≈ 1.7 × 10 20 G.

Low-energy effective theory
The leading order of the chiral perturbation theory Lagrangian for two quark flavors reads Here Σ is the unimodular and unitary 2 × 2 matrix field that includes the three physical pion degrees of freedom and m π is the pion mass. For the neutral pion, the physical value is m π ≈ 135 MeV. The minimal coupling of the charged pions to the electromagnetic field A µ is introduced through the covariant derivative, 2 , and τ 3 is the third Pauli matrix.
Throughout this paper, we will restrict ourselves to the neutral pion degree of freedom φ, that is, replace Σ → exp i fπ φτ 3 . The master Lagrangian that forms the basis of all the subsequent arguments, is then given by 2) The first two terms arise from the chiral perturbation theory Lagrangian (2.1). The third term stems from the anomalous Wess-Zumino-Witten coupling of Σ to a background electromagnetic field and takes the given simple form when restricted to the neutral pion degree of freedom [2]; we have defined for the sake of brevity. Finally, the last two terms in the effective Lagrangian (2.2) introduce the dynamical electromagnetic field including the standard gauge-fixing term.
In the conventional derivative expansion, the pion field φ itself counts as order zero, whereas any derivatives acting on it count as order one, and so does the pion mass m π . In contrast to the more common way of counting the powers associated with the electromagnetic field, we will treat A µ as an object of order zero just like φ. Note that this is consistent with gauge invariance thanks to the fact that we only consider neutral pions here, and also expresses the fact that we are looking at pion dynamics in strong background fields. The master Lagrangian (2.2) then represents the complete effective Lagrangian at the leading, second order of the derivative expansion. This way, including the Wess-Zumino-Witten term in the leading-order Lagrangian is made consistent with power counting.
For the record, we write down explicitly the equations of motion following from the Lagrangian (2.2). In the covariant relativistic notation, the equations for φ and A µ read In order to find the physical interpretation of the various modes in the spectrum, it will also be convenient to have at hand the nonrelativistic, gauge-invariant form of the classical equations of motion. Dropping the gauge-fixing term and trading the field strength tensor F µν for the electric and magnetic field E and B, the equations of motion take the form where dots stand for time derivative. These are the equations of axion electrodynamics [12].

Excitation spectrum
The excitation spectrum of a given system can be found using different approaches. In this section, we use for that purpose the field equations of motion, which give the best insight into the physical nature of the various excitation modes. Later on, we will rederive our result for the dispersion relations from the propagator of the pion and photon fields.
As the first step, we shift the magnetic field by the uniform background, B → B ex +B, and linearize the equations of motion (2.5) in the field fluctuations φ, E and B, Next, we carry out the Fourier transform to frequency-momentum space by introducing the conjugate variables ω and p, and the corresponding Fourier components of the fields, ϕ ω,p , e ω,p and b ω,p . 1 Upon using the Bianchi identity to eliminate b in favor of e, b = 1 ω p × e, and a brief further manipulation, the set of linearized equations of motion can be cast as This set of equations admits two classes of solutions: • One solution with ω γ (p) = |p|. This solution corresponds to the usual electromagnetic wave, characterized by ϕ = 0 and e · p = e · B ex = 0. Its electric component is thus linearly polarized in the plane perpendicular to the background field B ex .
• Two solutions with where p ⊥ is the component of momentum transverse to B ex . These solutions have nonzero ϕ as well as the electric field, related by The structure of the spectrum is best understood by focusing first on modes propagating along the magnetic field B ex . In this case, the residual SO(2) rotational symmetry ensures that helicity is a well-defined quantum number. Consequently, the pion and photon modes decouple. The dispersion relation ω − becomes degenerate with ω γ . Together, they define the two usual transverse electromagnetic waves, polarized in the plane perpendicular to B ex . Their propagation along the magnetic field B ex is unaffected by the anomaly as its contribution to the Lagrangian (2.2), proportional to E · B, now vanishes. The presence of two gapless degrees of freedom in the spectrum can thus be considered a consequence of gauge invariance just like in the vacuum of quantum electrodynamics (without background magnetic field).
The remaining excitation propagating along B ex corresponds to a ϕ-wave, accompanied by longitudinal fluctuations of the electric field, as dictated by eq. (3.5). The existence of this electric component follows from the fact that a gradient of the pion field induces nonzero electric charge density in presence of the magnetic background, see the equation of motion (2.5). This pion-like mode has a gap, which remains nonzero even in the chiral limit. This can be understood as Schwinger mass generation due to the 1+1-dimensional chiral anomaly, obtained from the 3+1-dimensional anomaly by dimensional reduction in presence of the background magnetic field. 2 Since the dispersion relations should be continuous functions of momentum, we can conclude from the above discussion that there will be two gapless and one gapped mode for arbitrary direction of propagation. The dispersion relation of the gapless ω − mode can be expanded in powers of momentum as where p is naturally the component of momentum in the direction of B ex . The anisotropy of the dispersion relation becomes maximal in the chiral limit where the p 2 ⊥ term vanishes and the dispersion relation in the transverse directions becomes quadratic.
At this point, we would like to remark that the observed anomaly-induced mixing of photons with a pseudoscalar (here the neutral pion) is of course not a new concept. It has in particular been known for a long time in the context of axion physics [13][14][15]. However, the physical scale hierarchy is quite different in that context, and the interplay of the extremely weak photon-axion interaction with the coherence of the photon beam results in a photon-axion conversion rather than full mixing; see also ref. [16] for a related discussion within neutrino physics. To the best of our knowledge, the actual mixing problem was first analyzed in ref. [17], and solved therein for the special case of propagation in the plane transverse to B ex . Our dispersion relations (3.4), valid for an arbitrary direction of momentum, are a generalization thereof.

Alternative derivation
It is instructive to rederive the dispersion relations (3.4) directly from the covariant equations of motion (2.4); this allows us to introduce already now some notation that we will make use of later in the calculation of the pion decay with. To make contact with that follow-up quantum-field-theoretic calculation, we keep the gauge-fixing term in the equation of motion for A µ , but select the Feynman gauge, ξ = 1, in which the equation takes a particularly simple form. 3 Denoting the Fourier component of A µ as a µ p , the linearized version of the equations of motion (2.4) can then be written compactly as The most straightforward way to solve these equations is to contract the second of them with n µ and thereby eliminate ϕ in terms of n·a. This leads to two classes of solutions, which are in a one-to-one correspondence with the solutions found above using the nonrelativistic equations of motion (2.5): • One solution with n · a = 0, for which ϕ = 0 and the momentum satisfies p 2 = 0. This is the relativistic linearly polarized photon.
• Two solutions with n · a = 0, for which momentum satisfies the covariant condition Given that n 2 = B 2 ex C 2 (p 2 − ω 2 ), it is easy to check that this reproduces the dispersion relations ω ± (p) in eq. (3.4).

Pressure at one loop
Pressure is one of the most important observables relevant for both model calculations and lattice simulations of the QCD phase diagram. In this section, we evaluate the thermal part of the pressure of our system; in other words, we do not investigate the effect of the magnetic field on the zero-temperature pressure of QCD. At the leading, one-loop order of the loop expansion, the pressure due to thermal excitations is given simply by where the sum runs over all physical excitations in the system and ω i (p) are their respective dispersion relations, here given by eqs. (3.3) and (3.4). This expression describes a free gas of noninteracting quasiparticles with dispersion relations fixed to their zero-temperature values. The feedback of the thermal medium into the spectrum only enters the pressure at the two-loop order through the quasiparticle interactions. Following our notation for the individual dispersion relations, we split the pressure into the contributions of the individual modes, P = P γ + P + + P − . The contribution of the relativistic photon mode with the linear dispersion relation is trivial to evaluate and amounts to the usual expression, As to the other two modes, we will in the following focus on the gapless mode ω − and discard the contribution P + . This is well justified at temperatures T m eff where this contribution is exponentially suppressed due to the gap of the ω + mode. By using spherical coordinates and integrating by parts with respect to the radial momentum, the momentum integral in eq. (4.1) can be rewritten in the dimensionless form where the dimensionless dispersion relationω − is defined bỹ θ is the angle between the momentum vector and the direction of the background field B ex , and the parameters α and τ are defined by While we cannot evaluate the thermal integral involving the dimensionless dispersion relationω − in a closed form, eq. (4.3) is suitable for an expansion in powers of τ , or equivalently in powers of temperature. (Recall that we assume the temperature to be small compared to m eff , and thus τ to be much smaller than one.) Both the radial and the angular integral is straightforward to carry out upon such expansion, thus leading to the series expansion of the pressure P − , with the expansion parameter B 2 ex C 2 T 2 /m 4 π . Interestingly, the expansion is simultaneously an expansion in powers of B ex : the pion mass in the denominators is the "bare" mass m π rather than the effective mass in the magnetic field, m eff . It is therefore mandatory to inspect separately what happens in the chiral limit where m π goes to zero.
The crucial observation is that in this limit, the dispersion relation in the transverse directions becomes quadratic in momentum, see eq. (3.7). Obviously, we have to treat the transverse and longitudinal directions separately, and it is therefore most natural to use cylindrical coordinates to carry out the momentum integration. In contrast to eq. (4.3), the pressure of the gapless mode can now be rewritten in the dimensionless form where x is the dimensionless longitudinal momentum, y is likewise the dimensionless transverse radial momentum squared, and Remarkably, the resulting two-dimensional integral can be evaluated in a closed form, for instance by introducing a new variable z via the substitution and subsequently using polar coordinates in the xz plane. The final result for the pressure due to the gapless mode ω − then reads As long as the contribution of the gapped mode ω + can be neglected, the full one-loop pressure of the system is given without further approximations by the closed expression (4.11) One can even say that this compact expression represents the complete asymptotic series expansion of the pressure at low temperatures. Namely, the contribution of the ω + mode is suppressed by the non-analytic Boltzmann factor e −m eff /T , as a consequence of which all of its Taylor coefficients at T = 0 vanish. We can conclude that the mixing of neutral pions with dynamical photons dramatically changes the thermodynamics of QCD at low temperatures and in strong magnetic fields: instead of the usual black-body scaling with T 4 , the pressure at temperatures T B ex C is dominated by a term that scales as T 3 B ex /f π . It is interesting to contrast this to the T 5/2 (B ex /f π ) 3/2 scaling at nonzero baryon chemical potential, found in ref. [7].

Anomalous pion decay
The presence of the strong background magnetic field affects also other physical observables than those relevant for equilibrium thermodynamics. In this section, we will focus on the decay properties of the neutral pion, bearing in mind that in the vacuum, the dominant, two-photon decay of the pion is one of the hallmarks of the chiral anomaly (see ref. [18] for a recent precision calculation of the vacuum anomalous pion decay rate). The magnetic field affects the neutral pion decay in several ways. First, it contributes to the pion mass through loop corrections [19][20][21]. Second, it affects, likewise through loop corrections, the pion decay constant f π [22]. Finally, it may open phase space for new decay processes, or increase the branching ratio of processes that in the vacuum are negligible compared to the two-photon decay [23,24]. Our discussion below focuses on the consequences of the chiral anomaly in a background magnetic field for the two-photon pion decay rate at tree level, in particular on the effects of the kinematic mixing of pions and photons. Non-anomalous loop corrections due to the magnetic field are not included here; their significance is discussed briefly in the concluding section 6.
For the sake of simplicity, we shall in this section refer to the ω + mode as the "pion", and denote the corresponding one-particle state with momentum p as |π, p . This is natural for this state smoothly interpolates to the vacuum neutral pion in the limit of vanishing magnetic field, B ex → 0. Likewise, we shall refer to the other two states as "photons", using the following notation: • |Γ, p for the "nonrelativistic" photon with the dispersion relation ω − (p).
Without further mentioning it explicitly, we shall use the Feynman gauge in which ξ = 1.
The two-photon decay of the neutral pion in principle includes three different processes: π → γγ, π → ΓΓ and π → γΓ; this corresponds to the sum over polarizations of the photons in the final state. However, only the last process, π → γΓ, is actually allowed. The reason for this is that QCD in a background magnetic field possesses a discrete symmetry, which can be thought of as modified parity, under which γ is even whereas π and Γ are odd [15]. Below, we therefore calculate the decay rate for the π → γΓ process alone and compare it to the vacuum decay rate of the neutral pion. (We have checked by an explicit calculation that the decay rates for the π → γγ and π → ΓΓ processes vanish.)

Relativistic photon polarization vector
In order to determine the probability rate for any decay or scattering process, it is mandatory to understand the properties of the asymptotic one-particles states. We already know the dispersion relations of all the one-particle states in our system. However, the calculation of decay rates or scattering cross-sections also requires the knowledge of the corresponding wave functions, or in other words of how the one-particle states couple to the elementary fields in our EFT. We start here with the discussion of the relativistic photon case, which is the most subtle of the three one-particle states as it is affected by gauge ambiguities.
Recall the field equations of motion in the relativistic notation and the Feynman gauge, eq. (3.8). In the second quantization, the solutions to these equations of motion enter the plane-wave expansion of the pion and electromagnetic fields. In particular, the relativistic photon mode has ϕ = 0, and thus only appears in the expansion of the field A µ via where a † p is the creation operator of the corresponding one-particle state, |γ, p = a † p |0 , a p is the associated annihilation operator, and µ (p) is the polarization vector of the oneparticle state. The four-momentum p in the integral is on-shell, that is, p 0 = ω γ (p) = |p|. The polarization vector µ must satisfy the following conditions: • The transversality constraint n · = 0, following from the property n · a = 0 of the relativistic photon solution to the classical equations of motion (3.8).
• The Feynman gauge constraint p · = 0. Note that upon using the on-shell condition p 2 = 0, this constraint is invariant under the gauge transformation µ → µ +λp µ with arbitrary complex λ; the same gauge invariance property must apply to all physical observables.
The freedom to shift µ without affecting physical observables can be used to set 0 = 0 without loss of generality. We thus have altogether three linear constraints on µ , which in four spacetime dimensions determine µ uniquely up to an overall factor. We note that given the explicit expression for the vector n µ , eq. (3.9), it is easy to find an explicit solution to these constraints in a covariant form, The overall normalization will be fixed using a different argument in the next subsection.

Vacuum transition amplitudes
In presence of field mixing, the one-to-one correspondence between elementary fields and asymptotic one-particle states is lost. Physical scattering amplitudes then have to be extracted from the Green's functions of the fields using the Lehmann-Symanzik-Zimmermann reduction formula. To that end, we need the vacuum transition amplitudes 0|χ i (0)|n, p , where χ i (0) runs over all elementary fields of the theory and n over all one-particle states. Consider the matrix propagator of the elementary fields, that is, the time-ordered two-point Green's function, defined by D ij (x−y) ≡ −i 0|T {χ i (x)χ j (y)}|0 . By inserting the partition of unity in terms of the eigenstates of the Hamiltonian, one arrives at the Källén-Lehmann spectral representation of the propagator in its nonrelativistic form [25], where ω n (p) denotes the energy of the Hamiltonian eigenstate |n, p . These eigenstates are assumed to be normalized according to n, p|m, q = δ mn δ 3 (p − q). For one-particle states, n takes discrete values and the propagator correspondingly has simple poles at energies given by the dispersion relations of the one-particle states. For multi-particle states, on the other hand, the label n is continuous, resulting in a branch cut in the propagator. Equation (5.4) tells us that the vacuum transition amplitudes connecting elementary fields to a given one-particle state can be extracted from the residuum of the associated pole in the propagator. The inverse propagator can be directly read off the Lagrangian (2.2). It is a 5 × 5 matrix and upon Fourier transform to momentum space and setting ξ = 1 takes the form Note that we used a somewhat unusual notation wherein the first row and column corresponds to the φ field, and the remaining four rows and columns correspond to the A µ field and carry the index µ and ν, respectively. Inverting this 5 × 5 matrix and using the fact that n · p = 0 leads to an expression for the propagator, The pole structure of the propagator naturally reproduces our previous results for the quasiparticle dispersion relations, see e.g. eq. (3.10). Let us first have a look at the poles corresponding to the two nonrelativistic modes, π and Γ. Comparing eq. (5.6) to the general spectral representation of the propagator (5.4), we find that up to an overall phase, where we defined 8) and temporarily denoted the states |π, p and |Γ, p as |+, p and |−, p , respectively, in order to keep the expressions compact. The coupling of the fields to the relativistic photon state γ can likewise be extracted by looking at the pole at p 2 = 0 in the propagator (5.6). Here the situation is, however, complicated by the fact that the propagator in covariant gauges includes contributions from states in the unphysical sector of the Hilbert space. This is also easily seen from the fact ¡ k q p π Γ γ Figure 1. Kinematics of the π → γΓ process, indicating our notation for the momenta of the three states; the arrows indicate the flow of momentum. We only consider the decay of pion at rest here, in which case k µ = (m eff , 0). that the lower-right 4 × 4 corner of the propagator is proportional to −g µν + n µ n ν /n 2 at p 2 = 0, which up to a sign is a projector to the three-dimensional space of vectors orthogonal to n µ , although there is only one physical state with the dispersion relation given by p 2 = 0. The way out is to notice that by the argument of the preceding subsection, 0|A µ (0)|γ, p must be proportional to µ (p). This fixes the one-dimensional subspace, corresponding to the physical photon polarization. The overall normalization of the vacuum transition amplitude can then be fixed using eqs. (5.4) and (5.6), leading to 0|φ(0)|γ, p = 0, Equations (5.7) and (5.9) form the basis for our calculation of the pion decay rate below.

Kinematics
Before we move on to the calculation of the amplitude for the pion decay, we first discuss the kinematics of the π → γΓ process. Since Lorentz invariance is explicitly broken by the presence of the background magnetic field, the decay rate in principle has to be evaluated as a function of velocity or momentum. We limit our attention for the sake of simplicity to decay of a pion at rest, since this assumption, as we will see, leads to a simple, semi-analytic expression for the decay rate. We expect the result to give a reasonable approximation also for a nonzero velocity of the pion provided that it is much smaller than the speed of light. The kinematics of the decay process is displayed in figure 1. Momentum conservation in the rest frame of the pion leads trivially to p + q = 0. Imposing the energy conservation condition, m eff = |p| + ω − (q), then gives the magnitude of momentum of the particles in the final state as a function of the angle θ with respect to the magnetic field B ex , 10) using the notation introduced in eq. (4.5).
In the next subsection, we will calculate the amplitude M for the π → γΓ decay at tree level, including the corresponding vacuum transition amplitudes that couple fields to ¡ q p ν µ = −iC µναβ p α q β . Figure 2. Feynman rule for the interaction between pion and electromagnetic fields, following from the Lagrangian (2.2). The dashed line stands for the φ field, whereas the wavy lines for the A field. The momenta p and q flow out of the vertex.
one-particle states. The differential decay rate for the process then reads (5.11) The δ-function, imposing energy and momentum conservation, reduces the phase space integration to an angular integration over directions of momentum p, where |p| is determined by eq. (5.10).

Decay rate
The decay of the pion into a photon pair is driven by a single interaction vertex in the Lagrangian (2.2), containing both fields. The corresponding Feynman rule is shown in figure 2. The calculation of the decay amplitude is, however, complicated by the kinematic mixing between the fields. As a consequence, the φ field in the interaction operator can couple to both the π in the initial state and the Γ in the final state, and so can one of the A fields. The other A must necessarily couple to the γ in the final state, since 0|φ|γ, p = 0. The decay amplitude therefore consists of two contributions, Next one inserts the vacuum transition amplitudes from eqs. (5.7) and (5.9) and takes the squared absolute value of the amplitude. What follows is a rather lengthy calculation, including manipulation of products of Levi-Civita tensors and kinematic properties of the momenta k, p, q, the vectors n µ k and n µ q , and the polarization vector µ p . At the end of the calculation, a very compact result surfaces, (5.14) Note that in the limit B ex → 0, this expression further simplifies to C 2 m π /[8(2π) 9 ], and upon trivial angular integration following eq. (5.12) gives which agrees with the well-known expression for the neutral pion decay rate in the vacuum.
For nonzero values of the background field, one simply has to insert the squared amplitude in eq. (5.12) and simplify the result. It is natural to normalize the decay rate by its value Γ vac in zero magnetic field. Upon some manipulation, we thus obtain the final result for the magnetic field dependence of the decay rate of the neutral pion, where the angular brackets indicate angular averaging over the full solid angle corresponding to the variable θ, and the angle α is defined in eq. (4.5). 4 The dependence of the pion decay rate on the magnetic field is shown numerically in figure 3. The increase of the decay rate with increasing B ex is, in fact, a result of a competition of two effects. First, the phase space for the decay products is increased due to the increase of the pion mass from m π to m eff ; this corresponds to the leading 1/ cos 3 α factor in eq. (5.16). Second, the anisotropic kinematics reduces somewhat the result by the angular average factor in eq. (5.16). Obviously, the anomalous contribution to the pion mass plays a dominant role for the pion decay.
Apart from the full numerical result, analytical approximations for the decay rate may also be of some interest. Given eq. (5.16), it is straightforward to obtain power expansions for the decay rate in both weak and strong magnetic fields, The latter expansion is particularly relevant for the chiral limit where only the leading term survives and we get a closed expression for the decay rate, this time in absolute units, 5

Summary and discussion
In this paper, we have analyzed the low-temperature thermodynamics of QCD in strong magnetic fields. This is dominated by neutral pions and photons since the charged pions acquire a large gap due to Landau level quantization. We showed that the chiral anomaly leads in presence of the background magnetic field to mixing of neutral pions and photons, and worked out the consequences of this mixing for physical observables. 4 In fact, the angular averaging indicated in eq. (5.16) can be carried out analytically in a closed form.
However, the resulting expression is rather cumbersome, and we therefore prefer the simple form of eq. (5.16); the angular average can, if needed, be done numerically with no effort. 5 In the chiral limit, the phase space for the pion decay is closed in the vacuum, see also eq. (5.15). The chiral anomaly then provides two key ingredients that make the pion decay possible in background magnetic fields: both the interaction with photons and the phase space by giving the pion a nonzero mass.  Figure 3. Thick red line: dependence of the pion decay rate on the scaled magnetic field, given by eq. (5.16). The decay rate is normalized to unity in zero magnetic field; the absolute value of the decay rate in the vacuum is determined by eq. (5.15). Dashed black line: polynomial approximation for the decay rate given by the first two terms in eq. (5.17). Solid black line: the decay rate without the 1/ cos 3 α prefactor, that is, just the anisotropy factor in eq. (5.16). While all the numerical values only depend on the dimensionless ratio B ex C/m π , the range on the horizontal axis was chosen so that, for physical values of m π and f π , its upper limit corresponds to B ex ≈ 10 20 G.
Our first main result is an expression for pressure of the system in the leading, one-loop approximation, see eqs. (4.6) and (4.10). The softening of the dispersion relation of one of the photon polarizations due to the mixing leads to an enhancement of pressure at low temperatures, which is most dramatic in the chiral limit, where the leading contribution to pressure scales as B ex T 3 /f π .
Our second main result is a formula for the dependence of the neutral pion decay rate on the magnetic field, eq. (5.16). The effect of the magnetic field is again most dramatic in the chiral limit. The B 3 ex /f 5 π scaling, displayed in eq. (5.19), is actually easy to understand. Namely, at tree level, the decay rate must consist of a factor C 2 from the interaction vertex times a kinematic function of the product B ex C, entering the dispersion relations of the pion and photon. Dimensional analysis then fixes the powers of B ex and f π in the final result. The nontrivial part of our result therefore is the numerical factor in eq. (5.19).
We would now like to discuss some of the assumptions and approximations underlying our analysis in the form of a set of concluding remarks.
First of all, we worked for the sake of simplicity strictly at the tree level, that is, we neglected one-loop corrections to the vacuum pion mass in presence of a magnetic field, see e.g. ref. [21]. While those are necessarily proportional to the pion mass itself in accord with the chiral symmetry, the anomaly makes the pion massive even in the chiral limit. Hence its  effect will certainly be dominant in sufficiently strong magnetic fields, or sufficiently close to the chiral limit. At the physical point, the anomaly becomes dominant in moderate fields, B 0.1 GeV 2 , see figure 4. Second, we only kept neutral pions in our EFT, which assumes that there is sufficient scale separation between the neutral and charged pion sectors. The effective mass of the charged pion in the background magnetic field is determined by the standard Landau level problem, m π ± (B ex ) = m 2 π + B ex . The graphical illustration of the numerical values of the neutral and charged pion masses as a function of the magnetic field in figure 5, makes it clear that there is a large range of magnetic fields in which the requirement of scale separation is satisfied.
In ref. [24] it was argued that in sufficiently strong magnetic fields, the dilepton decay π 0 → e + e − will be the dominant decay channel for neutral pion. However, it seems that these authors only included the effect of the magnetic field on the amplitude for such decay, not its consequences for the phase space of the decay products. Namely, the energy levels of the electron-positron pair also undergo Landau level quantization. By spin conservation, the decay is only possible into one gapless and one gapped fermion. Hence the threshold energy for the π 0 → e + e − decay channel to be open altogether equals m e + m 2 e + 2B ex , where m e is the vacuum electron mass. The position of the threshold is indicated by the blue line in figure 5, which makes it clear that the dilepton channel is actually closed in  Figure 5. Comparison of various energy scales in our system as a function of the external magnetic field. The thick black line corresponds to the anomalous neutral pion mass, eq. (3.6). The red line stands for the charged pion mass, lifted by the magnetic field as a consequence of Landau level quantization. Finally, the blue line indicates the threshold for the dilepton (e + e − ) decay of the neutral pion in presence of the external magnetic field. There is clearly a range of magnetic fields in which the pion spectrum features scale separation, that is, the neutral pion is considerably lighter than the charged pion. For illustration, demanding that the ratio of the pion masses is at most 1/2 requires that 0.06 GeV 2 B ex 3.2 GeV 2 , or 10 19 G B ex 5 × 10 20 G. The numerical results were obtained using the vacuum mass m π ≈ 135 MeV.
most of the range displayed therein, except for fields below ca 10 18 G. Finally, we neglected nonlinear effects within electrodynamics, which are induced by loop corrections and in presence of the background magnetic field lead to vacuum birefringence and photon splitting [26][27][28]. While these effects in general modify the photon polarization tensor, and thus its propagation in the background field, they will not affect the main qualitative conclusions of our paper.