Melting holographic mesons by cooling a magnetized quark gluon plasma

We extend our analysis of holographic meson dissociation in the presence of an intense magnetic field. In addition to the previously known critical temperature above which the mesons melt, we found that for certain magnetic field intensities there exists a second lower critical temperature, below which stable mesons cease to exist. While we showed before that there is a range of high temperatures for which mesons can be melted by changing the magnetic field intensity, here we show that, as a consequence of the second critical point, there is also a range of low temperatures for which this phenomenom, which we term Magnetic Meson Melting (MMM), can be triggered. Additionaly, we also show that the magnetic field decreases the mass gap of the meson spectrum along with their masses. We are able to observe this by constructing a configuration that makes it possible to apply gauge/gravity methods to study fundamental degrees of freedom in a quark-gluon plasma subject to a magnetic field as intense as that expected in high energy collisions. This is achieved by the confection of a ten-dimensional background that is dual to the magnetized plasma and nonetheless permits the embedding of D7-branes in it. The main difference with previous approaches, which in consequence gives the novel results, is that the magnetic field retroacts in the geometry itself, as opposed to be confined to the world volume of the probe D7-branes.

sufficiently high temperatures the opposite behavior is observed, and the quark condensate decreases as the magnetic field intensity increases. This phenomenom is now known as 'inverse magnetic catalysis' (IMC).
The dependence of the meson spectrum on the magnetic field has been studied in different frameworks. Chiral perturbation theory [12] as well as linear sigma models [13] show that charged and neutral pions respond oppositely to the magnetic field. While the former increase their masses in the presence of the magnetic field, the latter find their masses reduced. However, this doesn't means that all neutral mesons behave the same. Using lattice QCD it has been shown [14] that the neutral pion meson mass decreases monotonically as the magnetic field grows, while for the ρ-meson the oposite behavior is observed. Different results are also obtained for neutral heavy mesons using a potential model with constituent quarks [15].
The non-perturbative regime of the physics just described has been studied using the gauge/gravity correspondence [16]. It should be noted however, that the theory that can be studied through holographic methods is not properly quantum chromodynamics (QCD), because the exact gravity dual to this theory has not been found yet. Instead, what has been done is to consider theories similar to QCD, such as finite temperature N = 4 super Yang-Mills (SYM) with gauge group SU (N c ), and modify it to bring it as close to QCD as possible.
A relevant step in this direction was the inclusion of a small number of flavor degrees of freedom N f , or in other words, fields in the fundamental representation of the gauge group SU (N c ). In an abuse of language, we are going to refer to these fields generically as 'quarks' even if bosonic and fermionic degrees of freedom are included. Holographically, this is done by adding N f probe D-branes to the desired gravitational background [17], keeping N c N f to avoid backreaction. In particular, in order to add fundamental matter to SYM N = 4 it is necessary to study the dynamics of probe D7-branes on AdS 5 ×S 5 . The configuration is such that the D7-branes extend along the AdS 5 directions while wrapping a S 3 ⊂ S 5 .
In [18][19][20] it was found that if the background was deformed to allow a finite temperature, the system showed an interesting thermodynamic behavior. It was discovered that, for a critical value of the temperature T f un , the D7-branes undergoes a first order phase transition, taking place between a phase where they lie completely outside the horizon, in what it is called a Minkowski embedding, and a phase in which they fall into it, the black hole embedding. From the dual gauge theory perspective, this corresponds to a transition between a discrete meson spectrum (quark-antiquark bound states) and a gapless distribution of excitations. It is important to remark that what we have just described is a 'melting' or dissociation phase transition, not a confinement-deconfinement one. The holographic description of a confining phase involves a horizon-free geometry. At T deconf the gluons and adjoint matter become deconfined, at which point the dual geometry develops a black hole horizon. However, if the quark mass is large enough it is possible for the branes to lie outside the horizon and thus mesons are stable (to leading order within the approximations of large N c and strong coupling) for T deconf < T < T c , behavior consistent with lattice simulations [21][22][23].
The incorporation of a magnetic field in this setup was done in [24,25]. In these works the authors considered a pure gauge B-field along the gauge theory directions, which then enters as an excitation in the DBI action of the probe branes and thus it is equivalent to a magnetic field on the world-volume. The results obtained in [25] showed that the effect of the magnetic field was to increase the melting temperature of the mesons, to the point of not having dissociation at all for field intensities above a certain critical value. Not only that, but it was also found that the presence of the magnetic field induces chiral symmetry breaking, exhibited by a non-zero quark condesate even at zero bare quark mass. Regarding the meson spectrum, it was found that the magnetic field breaks the degeneracy of the levels given in [18] displaying a Zeeman like effect. It is important to remark that the excited modes on the S 3 were not considered in [24,25].
Another holographic setup for the magnetized quark gluon plasma was given in [26], where all the matter fields were in the adjoint representation. The five-dimensional theory presented there, that considers the full backreaction of the magnetic field, is dual to the desired gauge theory because it is a solution to a consistent truncation of supergravity (SUGRA) IIB [27]. While this five-dimensional truncation is enough in many scenarios, in order to add flavor degrees of freedom to the theory we require the ten-dimensional uplift. Embedding D7-branes into this background has proved difficult, because the compact part of the resulting geometry warps in a way that prevents an easy identification of the right 3-cycle that the D7-brane must wrap (see App. A).
We followed a different approach in [28], where we implemented a new consistent fivedimensional truncation. The family of solutions that we found features the full backreaction of a constant magnetic field and a scalar field dual to an operator of scaling dimension 2. From the five-dimensional perspective, we found that the presence of this scalar field induces the existence of a critical intensity for the magnetic field b c , above which the system becomes unstable. For intensities below b c , two branches of solutions exists, with one of them being thermodynamically preferred over the other. From the ten-dimensional perspective, the inclusion of this scalar field allows the compact space of the uplifted geometry to factorize as a warped 3-sphere, a 1-cycle and an angular coordinate θ that distributes the volume between these two spaces. Given this factorization, the probe D7-brane naturally wraps the 3-sphere while its embedding is described by the radial dependence of θ.
In [29] we presented our first results using the construction just described. We found that, in contrast to the results obtained with the setup from [25], the effect of the magnetic field is to decrease the critical temperature at which the mesons melt, that is, we observed IMC for meson dissociation. We also showed that the melting can be triggered at fixed temperature by adjusting the magnetic field intensity, but only for a certain range of temperatures. We named this phenomenom magnetic meson melting (MMM) in [29].
The main objective of this manuscript is to provide additional details and extend the analysis that we presented in [29]. Concretely, here we present a more in depth thermodynamic analysis and widen our study of the spectrum of mesons by considering excitations over the 3-sphere. We also give more details on the numerical construction of the D7-brane embeddings and the holographic renormalization of its action. Our main new result is that for a certain range of magnetic field intensities, in addition to the previously known critical temperature above which the mesons melt, there is a second critical temperature below which stable mesons cease to exist. While for the first hot temperature the effect of the magnetic field was the one of IMC, we observe MC for the new cold critical temperature 1 . As we will see, the transition is of first order in both cases. We will also show that, in addition to decrease the mass of the mesons, the magnetic field impose restrictions over the quantum numbers of the mesons over the 3-sphere.
The manuscript is organized as follows. In Sec. 2 we review the gravitational background, both from the 5D and 10D perspectives, while in Sec. 3 we introduce the probe D7-branes. In Sec. 4 we show the phase diagram of the flavor degrees of freedom, and in Sec. 5 we present the thermodynamic analysis. In Sec. 6 we compute the meson spectrum and we close in Sec. 7 with a discussion of our results. Some of the more technical details of our computations are contained in a series of appendices.

The gravitational background
The gravitational background is a solution to ten-dimensional SUGRA IIB that assymptotes AdS 5 × S 5 and admits a deformation that encodes the dual of a magnetic field in the gauge theory. Additionally, in order to be able to add flavor degrees of freedom it is crucial that the geometry its such that it permits to factorize the compact part of the space in the way described in the introduction.
It turns out that the line element that achieves this is of the form where ∆ is a wrapping factor given by

2)
ds 2 5 is the line element of a non-compact five-dimensional space that can be written as is the line element of a 3-cycle given by which for vanishing A reduces to the line element of a unit 3-sphere given in Hopf (toroidal) coordinates [30,31]. In what follows we will take, without loss of generality, L = 1. This implies that t'Hooft coupling λ is related to the string length l s by The r coordinate appearing in (2.3) meassures the radial distance such that we expect that the geometry assymptotes AdS 5 × S 5 as r → ∞. The 1-form A parametrizes an infinitesimal rotation along a periodic direction of the compact manifold that, in turn, codifies the internal degrees of freedom of the dual gauge theory. By keeping A and its exterior derivative F = dA in the cotangent space to the directions dual to those of the gauge theory, it will represent a U(1) vector potential. In order to introduce the desired magnetic field in the dual gauge theory we will look for solutions that admit and thus F will respresent a constant magnetic field in the z direction with intensity b.
The two properties of (2.1) that permit an easy embedding of a D7-brane are that the metric does not depend on the angular coordinate φ and that, according to (2.1), the direction that it represents remains orthogonal to the rest of the spacetime. Notice that the inclusion of ϕ was crucial in order to achieve this orthogonality. It is bacause of this that the brane can be placed at a fixed location in φ. Concerning the 3-cycle in (2.1), we notice that its volume depends on the position θ and, regardless of the intensity of the magnetic field b, this cycle becomes maximal at θ = 0, while for b = 0 it reduces to S 3 . For non-vanishing b, the 3-cycle gets tilted towards the five dimensional non-compact part of the spacetime in a manner that is volume preserving within the eight dimensions of this two spaces together.
In App. A we show that the configuration just described is part of the general truncation anzats presented in [27] and thus it constitutes a solution to type IIB supergravity, as long as the self-dual five-form is given by (A.5) while the five-dimensional metric and fields are a solution to five-dimensional gauged supergravity. In [28] we found a family of such solutions to the effective five-dimensional theory governed by the action where g is the determinant of the five-dimensional metric and R is the five-dimensional Ricci scalar. Notice that G 5 = π/2N 2 c because we set L = 1. The equations of motion that can be derived from this effective action are Additionally, in order for the solutions to be part of the truncation ansatz it is necessary to impose the constriction F ∧ F = 0, which our family of solutions satisfies inmediately. It is also important to note that a non-vanishing magnetic field implies a non-constant scalar field ϕ(r) in order for it to be a solution to (2.8). This means that the gravitational model found in [26] cannot be recovered from ours for b other than zero.
Given that the equations of motion (2.8) are highly non-linear, we resorted to numerical methods to solve them for non-vanishing intensity of the magnetic field. The general procedure is described in detail in [28] and here we will just review some steps that will be relevant to the following calculations. The near-boundary behavior of the solutions is presented in App. B.
To solve (2.8) we perform a numerical integration from the location of the horizon r h , at which the metric function U (r) vanishes, up to the boundary at r → ∞ where the geometry asymptotes AdS 5 . The symmetries of the equations of motion (2.8) allow us to completely specify any member of the family of solutions by the position of the horizon r h , the value that the scalar field takes at r h , and a parameter B. This three quantities are parameters over which we have computational control. However, what we want is to specify any member of the family of solutions by the values that the different physical quantities take in the dual gauge theory.
The parameter B is related to the magnetic field intensity by while r h is related to the temperature by (2.10) Note that althought T only depends on r h , the magnetic field b is a function of the position of the horizon r h and the parameter B. Thus, in order to vary T while keeping b constant, both r h and B need to be fine-tuned.
Regarding ϕ, the general behavior of the scalar field build by the procedure just described near the boundary is where ϕ 0 and ψ 0 are coefficients that can be read from the asymptotics of a specific solution. This particular behavior implies that ϕ is dual to an operator O ϕ of dimension 2, and thus it saturates the BF bound [32,33]. In consequence, ψ 0 is dual to the source of the operator and ϕ 0 to its vacuum expectation value O ϕ [33]. From the gauge theory perspective, it makes sense to specify the source of the operator and then compute the vacuum expectation value that it generates in response to said source. Thus, we use the freedom to choose ϕ(r h ) in order to fix ψ 0 to a given value.
It was found in [28] that for any given ψ 0 there exists a critical value for the intensity of the magnetic field b c /T 2 above which a singularity develops outside the horizon, indicating that the dual state in the gauge theory is unstable. For a given magnetic field intensity b/T 2 below this critical value, two different solutions exist, each one corresponding to a different value for the dimensionless ratio O ϕ /T 2 . The one with the higher value for O ϕ /T 2 corresponds to a state with negative specific heat, higher free energy and lower entropy than the other, showing that the solutions with smaller O ϕ /T 2 are thermodynamically preferred.
Finally, in [28] it was shown that because of the presence of the magnetic field the dual theory features a conformal anomaly, given by the fact that the trace of the stressenergy tensor does not vanish for a non-zero magnetic field. This causes some dimensionless observables not to depend only on dimensionless quantities such as b/T 2 and introduces a scheme dependent energy scale. This will have consequences on some of our numerical results, as some of them will be scheme dependent. The physical conclusion however, will be true regarding of the chosen scheme.

Flavor D7-brane embeddings
In order to add flavor degrees of freedom to the theory, it is necessary to add probe D7branes to the ten-dimensional gravitational backgrounds just described. This system can be thought of as a D3/D7 system with N c D3-branes and N f flavor D7-branes, with N c N f to be able to describe the latter as probe objects moving in the fixed geometry generated by the former. These objects extend in the following directions: In general, the position of the D7-branes is described by φ and θ as functions of the eight world-volume coordinates. This functions are in turn determinated by the equations of motion derived from the DBI action 2 where T D7 is the D7-brane tension given by

3)
g D7 denotes the induced metric on the world-volume, and the integration is performed over the eight directions on which the D7-branes extends (3.1). However, it turns out that it is possible to set φ = 0, χ(r) = sin θ(r). (3.4) consistently with the equations of motion, meaning that this choice automatically satisfy the equation for φ and leaves a second order ordinary differential equation for the profile function χ which can be solved after suitable initial conditions are given. With this choice the induced metric over the D7-brane world-volume, which can be computed from (2.1), is given by (3.5) where the prime denotes the derivative with respect of r and the profile function χ also appears in the wrapping factor ∆ = X + χ 2 (X −2 − X). (3.6) The remaining equation of motion for χ can also be derived from the general DBI action (3.2) after substitution of (3.5) and integration along the compact directions where is the Lagrangian density. The equation of motion for χ is given by where (3.10) The near-boundary behavior of the embedding profile χ can be obtained by solving (3.9) for large r, which can be acomplished by using the asymptotic expansions for the metric functions and the scalar field given in App. B. The result is where U 1 is one of the coefficients appearing in the expansions (B.1) for the background metric, m is related to the bare quark mass by 3 and as we will show below, c is related to the quark condensate. In order to directly compare our results to the ones presented in [20] we will report them in terms of the rescaled mass which is related to the mass gap of the meson spectrum [18,[34][35][36].
By solving (3.9) across all the bulk, there are two different kinds of embeddings that can be obtained: one where the branes lie completely outside the horizon, called Minkowski embedding, and another one where they fall through it, called black hole embedding. The general solution to (3.9) for different values for the magnetic field needs to be computed numerically, given that the background itself is not known analytically. We now procced to explain how to perform said numerical integration to compute both types of embeddings.
Given that in the case of black hole embeddings the D7-branes fall into the horizon, the numerical integration is performed from there to the boundary. In practice, we solved the equation of motion (3.9) as a power series of r around r h (3.14) With this procedure it is possible to write all the undeterminated coefficients up to any desired order in terms of χ h and the values that the metric functions and the scalar field take at the horizon. Even if the equation of motion is of second order, our choice of coordinates is such that the value of χ at the horizon is fixed by requiring the solution to be regular at r = r h . Then (3.14) is used to provide initial data for the numerical integration, starting at r = r h + with r h all the way to the boundary at r → ∞.
For the case of Minkowski embeddings the D7-branes lie completely outside the horizon, localized at θ = π/2 for a given r i > r h . We then solve (3.9) as a power series of r around r i and write all the undeterminated coefficients in terms of the values that the metric functions take at r i . Again, the value for χ (r i ) is fixed by requiring the solution to be regular at r = r i . By using (3.15) as initial conditions it is possible to numerically integrate (3.9) from r = r i + to the boundary r → ∞.
Finally, given that we wish to study the effect of the magnetic field and the temperature on the quarks, we need to work at fixedM . In practice, this is done by choosing the value of χ h or r i that gives the desired value of m (and in consequence ofM by (3.13)) for any given values of T and b.

Phase diagram
Following the procedure described in the previous section, it is possible to compute the D7-brane embeddings for any b/T 2 below its critical value. In the following we fix the source term for the scalar operator ψ 0 equal to zero, which implies that the maximum intensity for the magnetic field that the background can hold is b c /T 2 ≈ 11.24.
We present the phase diagram of the system in Fig. (1), which visually summarizes many of our results. In practice we generated the phase diagram fixingM = 1, but we explicitly checked that it is independent of the chosen value forM . The red curve indicates the critical value b/T 2 = 11.24, thus the region to its left it's not accesible with our current gravitational backgrounds. Note that we cannot explore the region given by T /M = 0 and b/M 2 = 0, which is exactly the physical scenario considered in [24]. Nonetheless, we can study the effect of any value of b/M 2 and T /M to the right of the red curve, which is included in the region considered in [25]. The blue curve denotes the interphase between the two types of embeddings.
Many interesting phenomena can be read from the phase diagram. The first can be seen by keeping b/M 2 fixed. For b/M 2 = 0 we recover the familiar results from [20], where from T /M = 0 up to the melting temperature T hot /M = 0.765, the D7-branes are in the Minkowski phase, while for higher temperatures the D7-branes are in the black hole phase. As long as the intensity of the magnetic field is such that the value of b/M 2 remains smaller than approximately 2.82, the behavior does not change much with respect to the b/M 2 = 0 case, in the sense that from the minimum T /M for that given b/M 2 up to a melting temperature T hot /M the branes are in the Minkowski phase, while for higher temperatures the branes are in the black hole phase. From the gauge theory perspective this means that, for 0 < b/M 2 < 2.82, stable mesons states exist for low temperatures while for high temperatures the mesons melt.
However, for 2.82 < b/M 2 < 3.44 something interesting happens: black hole embeddings exist for T /M near the minimum temperature. This new cold black hole phase, which is a consequence of the presence of the background magnetic field, is dual to a low temperature phase in which stable mesons also cease to exist, and hence their melting can alternatively be achieved by lowering the temperature. Thus, for 2.82 < b/M 2 < 3.44, in addition to the previous high melting temperature T hot /M , there exist now a low melting temperature T cold /M . For b/M 2 = 3.44 both temperatures coincide, and for b/M 2 > 3.44 no Minkowski embedding exist, meaning that stable meson states don't exist for high intensities of the magnetic field.
We would like to remark that our findings are very different, and even completely opposite when comparable, to the ones obtained using the setup from [25]. In particular, the cold phase transition is a novelty of our construction. While for vanishing magnetic field [25] also recovered the results from [20], they found that the melting temperature T hot /M increases with the magnetic field. This effect is such that, above a certain magnetic field intensity, there was no meson melting at all and the D7-branes were in the Minkowski phase regardless of the value of T /M . Thus while [25] found that the magnetic field prevents the melting of the mesons, here we found that it facilitates it.
The last affirmation is more evident if we return to the phase diagram in Fig. (1), but now exploring what happens if we keep the temperature fixed while changing the magnetic field. For temperatures below the intersection of the red and blue curves, that is T /M < 0.501, no black hole embeddings are permited and the branes are in the Minkowski phase for all intensities of the magnetic field up to its maximum value. The situation changes for 0.501 < T /M < 0.765, since for this temperatures the D7-branes transition from Minkowski to black hole embeddings at a melting magnetic field intensity, which we denote b M M M /M 2 . We named this phenomenom magnetic meson melting (MMM) in [29], because what happens from the dual theory perspective is that the mesons are being melted at a fixed temperature by increasing the magnetic field only. For higher values of the temperature, namely T /M > 0.765, only black hole embeddings exist regardless of the magnetic field intensity. From the gauge theory perspective, what happens is that the mesons are already melted because of the high temperature.
To close this section, in Fig. (2) we present the profiles of some characteristic embeddings as parametric plots in the (r cos θ, r sin θ)-plane, since this presentation makes it easier to visualize how the transition between both types of embeddings takes place. For this plot we fixed the position of the horizon at r h = 1/2, meaning that the temperature of the system is fixed at T = 3/4π. This is represented in Fig. (2) with a black line. However, each embedding in Fig. (2) is at a different value of T /M because each has a different value ofM .
We can see from

Thermodynamic analysis
In this section we will perform a thermodynamic analysis of the phase transition between Minkowski and black hole embeddings described above. This will include the computation of the quark condensate, free energy, entropy, and energy densities of the branes. According to the holographic dictionary, all of this quantities can be computed from the on-shell D7brane Euclidean action where the integration is performed over one period β = 1/T of Euclidean time t E = it and the metric is replaced by its Euclidean version. However, the result of this direct evaluation diverges as the integration is taken all the way up to the boundary at r → ∞. To deal with this issue it is necessary to substract this divergent behavior by adding covariant boundary terms to the action, in a process known as holographic renormalization [33,37]. We include the actual computation of the counterterm action in App. C, while here we just present the result expressed in the Fefferman-Graham coordinate u defined in (C.1), in terms of which the boundary is located at u = 0: where γ is the determinant of the induced metric at the boundary located at a radial cutoff u = . The action (5.2) contains the minimum terms that need to be added to the action in order to render it finite. However, we still have the freedom to add terms whose contributions by themselves are finite. In the current case the only ones are Choosing a specific value for the free coefficients C 1 and C 2 is equivalent to fixing the renormalization scheeme. As explained in [20], we can fix C 1 = 1/4 by demanding that the on-shell action vanishes for the supersymmetric embedding. Because we don't have a physical reason to choose a specific value for C 2 , some of our next results will be scheme dependent. Concretely, we will find that the free energy and energy densities are scheme dependent, but neither the quark condensate or the entropy are. Having said that, the full renormalized action, defined as is then finite in the limit → 0. We would like to remark that in order to properly evaluate (5.4) both the DBI action and boundary terms need to be expressed in the same radial coordinate. Given that our numerical solutions are naturally computed using the r coordinate, we will work using it in what follows.

Quark condensate
With the renormalized action at hand, we can compute the quark condensate by taking the variation of (5.4) with respect to the quark mass. The computation is performed explicitly in App. D, while here we just present the end result Once a given numerical embedding is known, it is possible to extract the value of the coefficients m and c by analyzing its near-boundar behavior. Thus, in practice m and c are both functions of b, T and either χ h or r i depending on the type of embedding. By inverting this relations we can eliminate the latter in favor of m and then express the quark condensate as a function of b/M 2 and T /M . Notice that, given that U 1 is a function of the magnetic field and the temperature, said term will affect the value of the quark condensate in a non-trivial way. The plots in this subsection were all produced at fixed M = 1. However, we explicitly checked that the result is independent of this choice as long as everything is plotted in terms of dimensionless quantities, as the quark condensate is scheme independent 4 . For any magnetic field intensity other than zero, the quark condensate never vanishes for the temperatures that we are allowed to explore with our gravity construction, which signalizes that chiral symmetry is broken. Additionally, the condensate diverges as T /M aproaches its minimum value for the given magnetic field 5 . Note that this happens whether the low temperature phase is Minkowski or black hole. This behavior is consistent with the fact that for that temperature the background itself is expected to undergo a phase transition, as explained in detail in [28]. A consequence of this is that we cannot check if the quark condensate vanishes at T /M = 0 for b/M 2 = 0.
In Fig. (4) we zoom-in near the region where the phase transition takes place. From Fig. (4) (a) we see that we recover the results from [20] when the magnetic field vanishes. For non-vanishing magnetic field intensities the curves retain the same shape, which sig-  nalizes that the phase transition remains first order regardless of the magnetic field. This is also true for the new cold phase transition, as it can be seen from Fig. (4) (c) and (e).
To conclude this subsection we show the quark condensate qq /κM 3 as a function of the magnetic field at fixed temperature in Fig. (5). The values of T /M displayed are in the range in which the phase transition can be triggered by changing the magnetic field intensity. We can see that the quark condensate decreases as the magnetic field increases for all the temperatures that we explored, which means that we are observing IMC in regard to the quark condensate. This is consistent with lattice results [8][9][10], as it is known that for some quarks IMC is observed at high temperatures. Lastly, in Fig. (6) we zoom-in near the region where the transition takes place. The shape of the curves signalize that this is a first order phase transition.

Free energy density
According to the holographic dictionary, the free energy of the system is related to the on-shell Euclidean D7-brane action (5.4) through Note that because there is no dependence on the gauge theory directions, it is possible to factorize the infinite three-dimensional volume vol(x) and a factor of 1/T coming from the integration over the Euclidean time from (5.4). While the remaining radial dependence in principle can be computed just by direct evaluation, for numerical purposes it is more convenient to rewrite (5.4) as a finite integral, instead of two infinite integrals with a finite diference.
To this end we express the counterterm action (5.2) and the finite term (5.3) using the near-boundary expansions of the fields (3.11) and (B.1) evaluated at r = r max . Given that  r max is understood to be taken to the boundary, we can keep only the leading terms. The result can be organized as where In order to create a finite integral, we rewrite the r max terms as an integral from r min to r max minus the lower integration limit terms 6 . This gives where we have denoted In conclusion, the final expression for the free energy that we evaluate numerically is where we normalized our result with respect to The integrand L + L ct goes to zero in the limit r max → ∞ by construction, so the integral has been regularized. Notice that fixing a specific value for the free coefficient C 2 will change the dependence of the free energy on the magnetic field. If we work at a fixed field intensity hovewer, this term amounts only to an additive constant.
In Fig. (7) we show F/NM 4 as a function of T /M for different values for b/M 2 at C 2 = 0. Note that the free energy seemingly diverges as T /M aproaches its minimum value for the given magnetic field, displaying the same behavior as the quark condensate. As in that case, the explanation is that the background itself is expected to undergo a phase transition at that point. In Fig. (8) we zoom over the parameter region at which the transition takes place. For b/M 2 = 0 the free energy displays the "swallow tail" shape characteristic of a first order phase transition (Fig. (8) (a)), which agrees with the findings in [20]. The same behavior is observed for b/M 2 = 2 ( Fig. (8) (b)). For b/M 2 = 3.3 and b/M 2 = 3.4 the the swallow tail seems to disappear both for the cold (Fig. (8) (c) and (e)) and hot (Fig. (8) (d) and (f)) phase transitions. However, from our analysis of the quark condensate we know that this is still a first order phase transition, even if the free energy difference between both phases is very small.
It is also important to note that the free energy is a scheme dependent quantity. This means that, even when F/M 4 is a dimensionless quantity, it will depend on b, T , and M independently, and not only on the two independent dimensionless ratios that can be built with them. The plots presented here were all computed forM = 1. Additionaly, the scheme can be changed by fixing a different value for C 2 . Given that the curves shown in Fig. (7) and (8) correspond to fixed magnetic field intensities, changing the renormalization scheme will only modify the curves by an additive constant. This in turn doesn't change the dependence of the critical temperature on the magnetic field nor the order of the phase transition 7 . We also show the free energy as a function of b/M 2 for four different values of T /M in Fig.(9). The values of T /M displayed are in the range in which the phase transition can be triggered by changing the magnetic field intensity. Again, we note that the free energy seemengly diverges as b/M 2 aproaches its maximim value for the given temperature. In Fig. (10) we zoom over the parameter region at which the phase transition takes place. This, along with our findings for the quark condensate, confirms that the transition is of first order.

Entropy and energy densities
The entropy is given by the derivative of the free energy with respect to the temperature at fixed magnetic field and quark mass As explained in Section 2, we don't have direct numerical control over b andM , as both quantities need to be extracted once a numerical solution is known. In consequence, we computed the derivative in (5.14) numerically, as opossed to the analytical aproach followed in [20]. We present the result of said numerical derivative in Fig. (11), where it can be seen that the general effect of the magnetic field is to increase the entropy, becoming infinite at the minimum temperature. It is important to note that the entropy is a scheme independent quantity, as the term proportional to C 2 in the free energy vanishes when taking the derivative (5.14). This is also realized in the fact that, like in the case of the quark condensate, Fig. (11) doesn't depend on the specific value ofM .
With both the entropy and free energy at hand, we can compute the internal energy by using the simple thermodynamic identity As the energy inherits the scheme dependence from F, we need to specify the value of C 2 .
In particular, the plots in Fig. (12) were generated for C 2 = 0 andM = 1. We can also read the qualitative behavior of the specific heat from Fig. (12), as it is given by The main conclusion from this is that, as can be seen from the slope of the curves, the specific heat becomes negative at sufficiently low temperatures for any non-vanishing magnetic field. Note that this is independent of the renormalization scheme, as the finite term vanishes when taking the derivative at fixed b. The temperature at which the specific heat becomes negative is displayed in Fig. (12) as a vertical line. This happens in either phase of the system, as for b/M 2 = 3.3 and b/M 2 = 3.4 cases the specific heat becomes negative in the black hole phase, while for b/M 2 = 2 it occurs in the Minkowski phase. A negative specific heat signalizes a thermodynamic instability, and thus, perturbing the system while in any of this states will cause it to decay into a stable one. However, this stable states, that should have positive specific heat at the same T /M and b/M 2 of the metastable state, are unknown to us since they do not appear as dual to any configuration of our gravity construction (see Fig. (12)).

Meson spectrum
From the gauge theory perspective, the phase transition just described takes place between a discrete meson spectrum (Minkowski embeddings) and a gapless distribution of excitations (black hole embeddings). Holographically, the mesons are related to stable excitations over the D7-branes. Thus in order to study the meson spectrum, it is necessary to study the dynamics of perturbations over the Minkowski embeddings.
Scalar perturbations are the ones related to vibrational modes of the D7-brane itself. A general excitation of this kind can be implemented as where x denotes the gauge theory spatial directions (x, y, z), Σ 3 denotes the coordinates of the 3-cycle (ψ, ϑ 1 , ϑ 2 ), 1 is a dimensionless parameter and χ (0) and φ (0) are the exact Minkowski solutions discussed earlier. In order to study these perturbations we substitute (6.1) and (6.2) in the most general EOM coming from the DBI action and solve them at first order in . This reveals that the χ (1) and φ (1) perturbations decouple and can be solved separately. We look for separable solutions of the form where ω and k are dual to the energy and momentum of the mesons. Given that Lorentz invariance is broken by the finite temperature of the plasma, the notion of meson mass is frame dependent. We will work with the same convention as [20] and define the meson mass in its rest-frame given by ω at k = 0.
By substituting (6.3) and (6.4) in their respective equations of motion, we obtain two equations, one for χ r and one for φ r , that only depend on the radial coordinate r, and a third equation that both F (Σ 3 ) and G(Σ 3 ) need to satisfy where ∇ 2 S 3 denotes the Laplace operator on the 3-sphere. Note however that this equation not only depends on the 3-sphere coodinates, but also on the radial coordinate and one of the gauge theory coordinates by the presence of the second term. This was not a problem for vanishing magnetic field, as the second term on the right hand side of (6.5) is eliminated if b = 0, and the expression reduces to Laplace equation on the 3-sphere as in [18,20,38]. Thus the solutions of (6.5) for b = 0 are the spherical harmonics Y l m 1 ,m 2 given by [30,31] Y l m 1 ,m 2 = C l m 1 ,m 2 (e iϑ 1 cos ψ) m 2 +m 1 (e iϑ 2 sin ψ) m 2 −m 1 P where C l m 1 ,m 2 is a normalization constant and P is a Jacobi polynomial. This spherical harmonics transform in the ( l 2 , l 2 ) representation of SO(4) and satisfy The quantum numbers (l, m 1 , m 2 ) are such that l ∈ N and |m 1,2 | ≤ l/2 vary independently on integer steps, hence taking integer or half-integer values according to the parity of l.
For any non-vanishing magnetic field both terms on the left hand side of equation (6.5) are present, and hence the only way that F , as a function of the angular coordinates exclusively, can be a solution is to satisfy on top of Laplace equation. By plugging (6.6) in (6.8) we conclude that this can only happen if F is a spherical harmonic with m 2 = 0. Given the range of values that this quantum number can take, this means that the presence of the magnetic field limits the angular momentum l of the perturbations on the 3-sphere to take only even values.
The remaining equations 8 for χ r and φ r are two independent ordinary second order differential equations, which depend on the frequency ω and the angular momentum on the 3-sphere l. By solving said equations near the boundary we obtain that the behavior of a generic solution is of the form χ r (r) = Ar −1+l + Br −3−l , (6.9) φ r (r) = Cr l + Dr −2−l , (6.10) where A, B, C, and D are free coefficients that can be read once a particular numerical solution its known. We are looking for the solutions that remain normalizable when we take the limit r → ∞, which are the ones for which A = 0 and C = 0 [18,20,38].
The numerical procedure to solve the equations for χ r and φ r is analogous to the one used in the previous sections. Given that the Minkowski embeddings are the ones for which θ = π/2 for a r i > r h , the equation of motion is degenerated at r i . Thus, we first solve the equations by a power series method around r i which allows to solve for the a j and b j coefficients in terms of the value that the metric functions, the scalar field and the unperturbed Minkowski embedding functions take at r i , as well as the frequency ω and the angular momentum l. Using this as initial conditions we integrate the equations for χ r and φ r numerically from r = r i + to the boundary at r → ∞. For fixed T , b andM , we then look for the frequencies that corresponds to normalizable solutions. This discrete set ω n,l,m 1 ,m 2 corresponds to the meson masses, where n = 0, 1, 2, ... denotes the radial quantum number.
In Fig. (13) we show the meson spectrum ω 2 normalized with respect toM 2 for the χ (left column) and φ (right column) perturbations, as a function of T /M for fixed magnetic field and different values for the quantum numbers n and l. Fig. (13) (a) and (b) display the case of vanishing magnetic field, which coincides with the values reported in [20]. It can be seen how the mass gap goes to zero as we aproach the melting temperature, and also how the meson masses increase as the temperature is reduced, asymptoting a constant value as T /M → 0. Fig. (13) (c) and (d) show the results for b/M 2 = 2. As we aproach the melting temperature T hot /M the meson masses go to zero, just like for vanishing magnetic field. However, as we reduce the temperature and aproach its minimum value the meson masses start decreasing without asymptoting to a specific value. Fig. (13) (e) and (f) display the results for b/M 2 = 3.4, for which there are two melting temperatures, T hot /M    and T cold /M . We can see that as we aproach both temperatures the meson masses go to zero.
We close this section by showing the meson spectrum as a function of b/M 2 at fixed  temperature in Fig. (14). While we specifically present the results for T /M = 0.65, we checked that the qualitative behavior is the same for any temperature in the range in which the phase transition can take place. We can see how the mass gap aproaches zero as the magnetic field goes to its critical value, confirming that the mesons can be melted by means of the magnetic field alone. Additionally, we recover the results from [20] for vanishing magnetic field.

Discussion
In this paper we employed holographic methods to study the effects of an external magnetic field over flavor degrees of freedom. While previously this magnetic field had been introduced as an excitation over the world volume of the probe D7-branes [24,25], here we followed a different aproach and considered a family of solutions in which the magnetic field backreacts in the geometry itself. Notably, the results obtained with our setup are very different, and even completely opposite when comparable, to the ones obtained in [24,25].
Some of our main results are conveniently summarized in the phase diagram in Fig.(1), where we can see that the presence of the magnetic field induces a very rich thermodynamic behavior. While previously [29] we had only observed IMC for meson dissociation, we see now that the physics of the system is much more interesting, as the dependence of the phase transition on the magnetic field is highly non-trivial.
One of the main novelties is that, for some magnetic field intensities, mesons can be melted by decreasing the temperature of the plasma. This can be seen explicitly from our thermodynamic analysis, as there is now a first order phase transition at low temperatures in addition to the expected first order phase transition at high temperatures. While for the first hot temperature the effect of the magnetic field is the one of IMC, we observe MC for the cold critical temperature. This is reminiscent of the effect of the magnetic field on the quark condensate for light quarks [8][9][10], as there are some temperatures for which both behaviors can be observed. However, it is important to remark that this behavior is not only a consequence of the magnetic field, but of its interplay with a scalar field that saturates the BF bound. It is possible that for low temperatures the energy of the whole system, consisting of the scalar field and the quarks, is such that it is unfavorable to have quark-antiquark bound states. In order to properly check that this is the case one would need to isolate the energy associated with the scalar field from the stress-energy tensor of the background solution. Nonetheless, this is not a trivial matter given the intricate interplay between the metric, scalar and magnetic field.
We would also like to remark that what we study here is neither a confinement/deconfinement transition nor a chiral symmetry breaking transition. This is because the existence of a horizon implies that the adjoint degrees of freedom are deconfined and chiral symmetry is broken as in both phases of our system the expectation value of the chiral condensate is different from zero for finite values of b/M 2 and T /M . In fact, lattice calculations [21][22][23] show that the temperature at which the mesons dissociate is higher than that at which deconfinement happens, thus it is expected that the transition we study does not coincide with any of the former two.
Another important result is that, for a certain range of temperatures, the mesons can be melted by just changing the intensity of the magnetic field. We first observed this in [29], where we named this phenomenom magnetic meson melting. For sufficiently high temperatures the mesons are melted regardless of the magnetic field. The opposite behavior is observed from the setup in [25], where it was imposible to melt the mesons for magnetic field intensities above a certain critical value.
From our results for the quark condensate in Fig. (5), we can see that the quark condensate decreases with the magnetic field for the temperatures that we explored. Thus we conclude that we observe IMC in regard to chiral symmetry. This is consistent with lattice calculations [8][9][10], which predict that for sufficientely high temperatures and light quarks IMC is observed.
By performing a thermodynamic analysis we showed that the phase transition is of first order at both critical temperatures. Interestingly, we also found that for any nonvanishing magnetic field intensity there exists a temperature below which the specific heat of the flavor degrees of freedom becomes negative. This of course signals a thermodynamic instability, but the stable state to which the system decays is not part of the states we have access with our construction. However, determining said state is beyond the scope of this paper and thus we leave that as future work.
By studying perturbations of the D7-brane position we computed the spectrum of scalar mesons. We showed that in the case of the Minkowski embeddings the meson spectrum is discrete, as it was known for zero magnetic field [18,20,38]. However, we found that a non-vanishing magnetic field has many important effects over the meson spectrum. The first one is that it reduces the number of possible excited states on the 3-sphere by imposing that the quantum number m 2 has to vanish for any b = 0. The second one is that the magnetic field reduces the meson masses for all T /M and any quantum number, as shown in Fig. (14). For the values of T /M for which the transition can be triggered using the magnetic field we showed that the meson masses, including its mass gap, tend to zero as the magnetic field reaches its melting value.
Our results for the meson spectrum are in agreement with those found in [13], where a one-loop calculation is done in a linear sigma model to find that the mass of the neutral pion is reduced by the application of a magnetic field. The same behavior was observed using lattice calculations in [14] for the same pion. The spectrum of the η meson was already computed in the holographic context [39]. However, their calculation indicates an increase in the mass with the magnetic field and the authors concluded that it should be corrected by taking backreaction effects into account. Similar corrections should be applied to [40], where the melting transition is studied out of equilibrium. Given that the magnetic field backreacts on the geometry of our construction, we provide such correction, showing that the effect of the magnetic field is the inverse to the one reported in [39].
It is important to remark that our gravitational setup does not allow us to explore all values of the temperature and the magnetic field, as shown in Fig. (1). However, our thermodynamic analysis shows that the quark condensate, free energy, entropy, and internal energy, diverge as we aproach the critical value of b/T 2 , showing that an interesting thermodynamic behavior is expected.

Acknowledgments
We acknowledge partial financial support from PAPIIT IN113618, UNAM. We thank Omar del Pilar and Aslan García for their help with the execution of numerical computations. We also thank David Mateos, Alejandro Ayala and Maria Elena Tejeda-Yeomans for useful discussion and comments.

A General truncation ansatz
The family of solutions to 10-dimensional SUGRA IIB we found is part of the general truncation ansatz presented in [27]. In this apprendix we will show how to obtain it and the family of solutions from [26] from this general setting. We will use the notation employed in [27], which is sligthly different than the one used in the main text. We will remind the reader when necessary.
The five-dimensional family of solutions presented in [27] can be uplifted to a tendimensional family of solutions to the SUGRA IIB equations of motion by considering the metric and self-dual five-form as the only non-vanishing fields. The truncation ansatz for the ten-dimensional metric is given by where L is the AdS radius, ds 2 5 denotes the line element of a 5-dimensional non-compact spacetime, A i are three independent Maxwell fields living in the non-compact part of the geometry, µ i are parametrized by two angular coordinates {θ, ψ} as µ 1 = sin θ, µ 2 = cos θ sin ψ, µ 3 = cos θ cos ψ, while ∆ is given by ϕ 1 and ϕ 2 are two independent scalar fields living in the non-compact part of the spacetime. Meanwhile the 5-form is given by where denotes the Hodge dual associated with (A.1) and where 5 and¯ are the volume form and the Hodge dual associated with the five-dimensional metric, and F i = dA i . Note that G 5 reduces to 5 when the scalar and Maxwell fields are taken equal to zero.
The family of solution we found can be recovered from the general truncation ansatz by taking 2 and the vectors a i as This in turn implies that With this choice the line element (A.1) simplifies to which is precisely the line element (2.1) from the main text after renaming the angular coordinates {φ, ψ, ϑ 1 , ϑ 2 } → {φ 1 , ψ, φ 2 , φ 3 }. As it is explained in the main text, this metric permits an easy embedding of a D7-brane because it doesn't depend on φ 1 and the direction that it represents remains orthogonal to the rest of the spacetime. This allows to put the D7-brane at a fixed position in φ 1 while it extends along the AdS directions and wraps the 3-cycle, which gets maximal for θ = 0. The family of solutions presented in [26] doesn't possesses this characteristics.
The family of solutions from [26] is also part of the general truncation ansatz described above. It can be recovered by taking which implies that X i = 1 and ∆ = 1 and thus the 10D metric is given by (A.14) Note that, crucially, the φ 1 coordinate is no longer orthogonal to the rest of the spacetime and thus the D7-brane cannot lie at constant φ 1 unless the embedding profile depends on the gauge theory coordinates.

B Boundary expansions of the background fields
Given that the background solution is naturally computed numerically using the r coordinate defined in (2.3), for many purposes it is convenient to know how a generic solution behaves at the boundary at r → ∞. This can be obtained by solving the equations of motion (2.8) by a power series method around the boundary. The only restrictions we impose is that the metric asymptotes exactly the metric of AdS 5 and the non-normalizable mode of the scalar field is turned off. The result reads where U 1 , U 4 , W 4 and ϕ 0 are coefficients that are not fixed by the equations of motion, but can be read as functions of the magnetic field intensity b and the temperature T once a particular numerical solution is known. As explained in detail in [28], ϕ 0 is dual to the VEV of the scalar operator O ϕ while U 4 and W 4 are both related to the stress-energy tensor.

C Holographic renormalization
As it is commonly the case in the gauge/gravity correspondence, the Euclidean DBI action (5.1) suffers from long-distance divergences when evaluated on-shell, making the variational principle ill defined. There exists a systematic method to deal with this issue, known as holographic renormalization, that has been studied extensively [33,37]. The first step is to write the near-boundary behavior of all the fields that enter in the Euclidean action, which is more conveniently done in the Fefferman-Graham coordinate where the five-dimensional metric (2.3) takes the form Next we solve the equations of motion (2.8) by a power series method around the boundary, located in this particular coordinate at u = 0. The resulting expansions for g, ϕ, and F are given by g ij (u) = g ij (0) + (g ij (4) + h ij log u + H ij log 2 u)u 4 + O(u 6 ), ϕ(u) = u 2 (ϕ (0) + ψ (0) log u + (ϕ (2) + ψ (2) log u + Ψ (2) log 2 u)u 2 ) + O(u 6 ) F uν = 0, F ij = F ij (t, x, y, z), and any not listed coefficient up to the specified order is equal to zero. Doing the same for the D7-brane profile χ(u) we obtain χ(u) = uχ (0) + u 3 χ (2) + O(u 5 ), (C. 4) where χ (0) and χ (2) are both free coefficients not fixed by the equations of motion.
The next step is to evaluate the Euclidean DBI action (5.1) on the solutions (C.2) and (C.4). After that we integrate from a radial cut-off to an arbitrary u max , which even if bigger than , is still close to the boundary and remains fixed. This allows us to collect all the divergent terms in the limit → 0, which can be organized as follows: where a (1) = 1 4 4 , a (2) = − 1 2 are independent of the scalar field, while are divergent terms that only depend on ϕ, and are depend only on the source ψ (0) . In order to properly substract this divergent terms from the action, it is necessary to invert the series for all the fields, and express the a (i) , b (i) , and c (i) terms in a covariant manner. To this end we first note that up to leading order in 2 γ ij = g ij (0) + . . . , γ ij 2 = g ij (0) + . . . , χ = χ (0) + . . . , (C.9) where the dots denote terms that vanishes in the limit → 0. This allows us to rewrite the a (i) terms in a covariant manner as a (1) = 1 4 4 , a (2) = − χ 2 2 4 , a (3) = −F ij F ij log 8 4 , (C. 10) where the indexes are raised and lowered using the boundary metric γ ij .
Let us turn now to the b (i) terms. Given that both ϕ (0) and ψ (0) appear at the same order in the expansion of ϕ, the inversions need to be performed carefully. To leading order we have ϕ 2 = ϕ (0) + ψ (0) log + . . . , (C.11) from where ϕ 2 4 = ϕ 2 (0) + 2ϕ (0) ψ (0) log + ψ 2 (0) log 2 + . . . , (C. 12) and discarding the finite term in the limit → 0 we conclude that b (1) = ϕ 2 12 4 . (C. 13) Note however that if the source term ψ (0) vanishes, the finite term we just discarded is the only contribution. Thus, for vanishing source, b (1) gives only a finite term which is fixed for ϕ (0) = 0. For b (2) , from (C.3) in (C.2) log (2ϕ (0) ψ (0) +ψ 2 (0) log )−ϕ (2) +. . . , (C.14) and then using (C.12) Discarding the finite term ϕ (2) we conclude that (C. 16) It turns out that for ψ (0) = 0 the term that goes like 1/ log vanishes and once again the ϕ 2 term gives only a finite contribution. Even if our case of interest is the one with no source for the scalar field, we conducted the analysis up to this point not setting ψ (0) to zero because it induces a fixed a finite counterterm that should be taken into account also in the sourceless case. Nothing similar will happen for the remaining c (i) terms, which vanish identically for ψ (0) = 0, so from here onwards, as in the main text, we will limit ourselves to that case.
Finally, by using (C.2) and (C.3) the expansion of the determinant of the boundary metric is found to be g (0) = 4 √ γ 1 + 1 6 ϕ 2 . (C.17) With this ingredients we can write (C.5) in a covariant manner, and thus the counterterms are given by where any term that goes to zero in the limit → 0 was discarded.

D Computation of the quark condensate
The quark condensate is given by the variation of the action with respect to the quark mass [37,41] where we have used (3.12) to express everything in terms of the parameter m. This can be rewritten as a variation with respect to the embedding profile as qq = π 2 λ lim