{'fraction_non_alphanumeric': 0.04729785140232901, 'fraction_numerical': 0.018625963588650155, 'mean_word_length': 3.9301561631374136, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Unquenched massive flavors and flows in Chern-Simons matter theories

We construct a holographic dual to the three-dimensional ABJM Chern-Simons matter theory with unquenched massive flavors. The flavor degrees of freedom are introduced by means of D6-branes extended along the gauge theory directions and delocalized in the internal space. To find the solution we have to solve the supergravity equations of motion with the source terms introduced by the D6-branes. The background we get is a running solution representing the renormalization group flow between two fixed points, at the IR and the UV, in both of which the geometry is of the form AdS_4 x M_6, where M_6 is a six-dimensional compact manifold. Along the flow, N=1 supersymmetry is preserved and the flavor group is Abelian. The flow is generated by changing the quark mass m_q. When m_q is very large we recover the original unflavored ABJM solution, while for m_q very small our solution becomes asymptotically equivalent to the one found recently for massless smeared flavors. We study the effects of the dynamical quarks as their mass is varied on different observables, such as the holographic entanglement entropy, the quark-antiquark potential, the two-point functions of high dimension bulk operators, and the mass spectrum of mesons.


Introduction
Three-dimensional Chern-Simons matter theories have been studied extensively in the last few years due to their rich mathematical structure and their connection with different systems of condensed matter physics. In particular, the Aharony-Bergman-Jafferis-Maldacena (ABJM) theory [1] has provided a highly non-trivial example of the AdS/CFT correspondence [2,3]. The ABJM theory is an N = 6 supersymmetric U(N) × U(N) gauge theory with Chern-Simons levels k and −k, coupled to matter fields which transform in the bifundamental representations (N,N) and (N, N) of the gauge group. The ABJM construction was based on the analysis of [4,5], in which the supersymmetric Chern-Simons theories were proposed as the low energy theories of multiple M2-branes. When N and k are large the ABJM theory admits a gravity dual in type IIA supergravity in ten dimensions. The corresponding background is a geometry of the form AdS 4 × CP 3 with fluxes (see refs. [6][7][8][9] for reviews of different aspects of the ABJM theory).
One of the possible generalizations of the ABJM theory is the addition of flavor fields transforming in the fundamental representations (N, 1) and (1, N) of the gauge group. In the supergravity description these flavors can be added by considering D6-branes extended along the AdS 4 directions and wrapping a three-dimensional submanifold of CP 3 . By imposing the preservation of N = 3 supersymmetry one finds that the D6-brane must wrap a RP 3 submanifold of the internal space [10,11]. When the number of flavors is small one can treat the D6-branes as probes, which is equivalent to the quenched approximation on the field theory side. This is the approach followed in refs. [12][13][14][15] (see also [16]).
In order to go beyond the quenched approximation, one must be able to solve the supergravity equations of motion including the backreaction induced by the source terms generated by the flavor branes. The sources modify the Bianchi identities satisfied by the forms and the Einstein equations satisfied by the metric. These equations with sources are, in general, very difficult to solve, since they contain Dirac δ-functions whose support is the worldvolume of the branes. In order to bypass this difficulty we will follow here the approach proposed in [17] in the context of non-critical holography, which consists of considering a continuous distribution of flavor branes. When the branes are smeared in this way there are no δfunction sources in the equations of motion and they become more tractable. Substituting a discrete set of branes by a continuous distribution of them is only accurate if the number N f of flavors is very large. Therefore, this approach is valid in the so-called Veneziano limit [18], in which both N and N f are large and their ratio N f /N is fixed. The smearing procedure was successfully applied to obtain supergravity solutions that include flavor backreaction in several systems [19][20][21] (see [22] for a detailed review and further references).
A holographic dual to ABJM with unquenched massless flavors in the Veneziano limit was found in [23]. In this setup the flavor branes fill the AdS 4 and are smeared in the internal CP 3 space in such a way that N = 1 supersymmetry is preserved. Notice, that since the flavor branes are not coincident, the flavor symmetry is U(1) N f rather than U(N f ). A remarkable feature of the solution found in [23] is its simplicity and the fact that the tendimensional geometry is of the form AdS 4 × M 6 , where M 6 is a compact six-dimensional manifold whose metric is a squashed version of the unflavored Fubini-Study metric of CP 3 . The radii and squashing factors of this metric depend non-linearly on the flavor deformation parameter N f N λ, where λ = N/k is the 't Hooft coupling of the theory. Moreover, the dilaton is also constant and, since the metric contains an AdS 4 factor, the background is the gravity dual of a three-dimensional conformal field theory with flavor. Actually, it was checked in perturbation theory in [24] that the ABJM theory has conformal fixed points even after the addition of flavor. This solution captures rather well many of the effects due to loops of the fundamentals in several observables [23]. Its generalization at non-zero temperature in [16] leads to thermodynamics which pass several non-trivial tests required to a flavored black hole.
Contrary to other backgrounds with unquenched flavors, the supergravity solutions dual to ABJM with smeared sources are free of pathologies, both at the IR and the UV. This fact offers us a unique opportunity to study different flavor effects holographically in a wellcontrolled setup. In this paper, we will study such effects when massive flavors are considered. The addition of massive flavors breaks conformal invariance explicitly and, therefore, the corresponding dual geometry should not contain an Anti-de Sitter factor anymore. Actually, for massive flavors the quark mass is an additional parameter at our disposal which we can vary and see what is the effect on the geometry and observables. Indeed, let m q denote the quark mass. In the IR limit in which m q is very large we expect the quarks to be integrated out and their effects to disappear from the different observable quantities. Thus, in the IR limit we expect to find a geometry which reduces to the unflavored ABJM background. On the contrary, when m q → 0, we are in the UV regime and we should recover the deformed Anti-de Sitter background of [23]. The important point to stress here is that the quark mass triggers a non-trivial renormalization group flow between two fixed points and that we can vary m q to enhance or suppress the effects due to the loops of the fundamentals.
To find the supergravity solutions along the flow, we will adopt an ansatz with brane sources in which the metric and forms are squashed as in [23]. By imposing the preservation of N = 1 supersymmetry, the different functions of the ansatz must satisfy a system of first-order BPS equations, which reduce to a single second-order master equation. The full background can be reconstructed from the solution to the master equation.
The flavor branes corresponding to massive flavors do not extend over the full range of the holographic coordinate. Indeed, their tip should lie at a finite distance (related to the quark mass) from the IR end of the geometry. Moreover, in the asymptotic UV region, the geometry we are looking for should reduce to the one in [23], since the quarks should be effectively massless in that region. Therefore, we have to solve the BPS equations without sources at the IR and match this solution with another one in which the D6-brane charge is non-vanishing and such that it reduces to the massless flavored solution of [23] in the deep UV. Amazingly, we have been able to find an analytic solution in the region without sources which contains a free parameter which can be tuned in such a way that the background reduces to the massless flavored geometry in the asymptotic UV. This semi-analytic solution interpolates between two different conformal AdS geometries and contains the quark mass and the number of flavors as control parameters.
With the supergravity dual at our disposal, we can study the holographic flow for different observables. The general picture we get from this analysis is the following. Let l be a length scale characterizing the observable. Then, the relevant parameter to explore the flow is the dimensionless quantity m q l. When m q l is very large (small) the observable is dominated by the IR unflavored (UV massless flavored) conformal geometry, whereas for intermediate values of m q l we move away from the fixed points. We will put a special emphasis on the study of the holographic entanglement entropy, following the prescription of [25]. In particular, we study the refined entanglement entropy for a disk proposed in [26], which can be used as a central function for the F-theorem [27]. We check the monotonicity of the refined entropy along the flow (see [28] for a general proof of this monotonic character in three-dimensional theories). Other observables we analyze are the Wilson loop and quarkantiquark potential, the two-point functions of high-dimension bulk operators, and the mass spectrum of quark-antiquark bound states.
The rest of this paper is divided into two parts. The first part starts in Section 2 with a brief review of the ABJM solution. In Section 3 we introduce the squashed ansatz, write the master equation for the BPS geometries with sources, and classify its solutions according to their UV behavior. In Section 4 we write the analytic solution of the unflavored system that was mentioned above while, in Section 5 we construct solutions which interpolate between an unflavored IR region and a UV domain with D6-brane sources. The backgrounds corresponding to ABJM flavors with a given mass are studied in Section 6.
In the second part of the paper we study the different observables. In Section 7 we analyze the holographic entanglement entropy for a disk. Section 8 is devoted to the calculation of the quark-antiquark potential from the Wilson loop. In Section 9 we study the two-point functions of bulk operators with high mass, while the meson spectrum is obtained in Section 10. Section 11 contains a summary of our results and some conclusions. The paper is completed with several appendices with detailed calculations and extensions of the results of the main text.

Review of the ABJM solution
The ten-dimensional metric of the ABJM solution in string frame is given by: where ds 2 AdS 4 and ds 2 CP 3 are respectively the AdS 4 and CP 3 metrics. The former, in Poincaré coordinates, is given by: where dx 2 1,2 is the Minkowski metric in 2+1 dimensions. In (2.1) L ABJM is the radius of the AdS 4 part of the metric and is given, in string units, by: where N and k are two integers which correspond, in the gauge theory dual, to the rank of the gauge groups and the Chern-Simons level, respectively. The ABJM background contains a constant dilaton, which can be written in terms of N and k as: Apart from the metric and the dilaton written above, the ABJM solution of type IIA supergravity contains a RR two-form F 2 and a RR four-form F 4 , whose expressions can be written as: where J is the Kähler form of CP 3 and Ω AdS 4 is the volume element of the AdS 4 metric (2.2). It follows from (2.5) that F 2 and F 4 are closed forms (i.e., dF 2 = dF 4 = 0). The metric of the CP 3 manifold in (2.1) is the canonical Fubini-Study metric. Following the approach of [23], we will regard CP 3 as an S 2 -bundle over S 4 , where the fibration is constructed by using the self-dual SU(2) instanton on the four-sphere. This representation of CP 3 is the one obtained when it is constructed as the twistor space of the four-sphere. As in [23], this S 4 -S 2 representation will allow us to deform the ABJM background by squashing appropriately the metric and forms, while keeping some amount of supersymmetry. More explicitly, we will write ds 2 CP 3 as: where ds 2 S 4 is the standard metric for the unit four-sphere, x i (i = 1, 2, 3) are Cartesian coordinates that parameterize the unit two-sphere ( i (x i ) 2 = 1) and A i are the components of the non-Abelian one-form connection corresponding to the SU(2) instanton. Let us now introduce a specific system of coordinates to represent the metric (2.6) and the two-form F 2 . First of all, let ω i (i = 1, 2, 3) be a set of SU(2) left-invariant one-forms satisfying Together with a new coordinate ξ, the ω i 's can be used to parameterize the metric of the four-sphere S 4 as: where 0 ≤ ξ < +∞ is a non-compact coordinate. The SU(2) instanton one-forms A i can be written in these coordinates as: Let us next parameterize the x i coordinates of the unit S 2 by two angles θ and ϕ (0 ≤ θ < π, 0 ≤ ϕ < 2π), Then, it is straightforward to demonstrate that the S 2 part of the Fubini-Study metric can be written as: where E 1 and E 2 are the following one-forms: (2.11) Therefore, the CP 3 metric can be written in terms of the one-forms defined above as: We will now write the expression of F 2 in such a way that the S 4 -S 2 split structure is manifest. Accordingly, we define three new one-forms S i (i = 1, 2, 3) as: Notice that the S i are just the ω i rotated by the two angles θ and ϕ. In terms of the forms defined in (2.13) the line element of the four-sphere is obtained by substituting ω i → S i in (2.7). Let us next define the one-forms S ξ and S i as: 14) in terms of which the metric of the four-sphere is ds 2 Moreover, the RR two-form F 2 in (2.5) can be written in terms of the one-forms defined in (2.11) and (2.14) as: The solution of type IIA supergravity reviewed above is a good gravity dual of the U(N) k × U(N) −k ABJM field theory when the AdS radius L ABJM is large in string units and when the string coupling constant e φ is small. From (2.3) and (2.4) it is straightforward to prove that these conditions are satisfied if k and N are in the range N 1 5 ≪ k ≪ N.

Squashed solutions
Let us consider the deformations of the ABJM background which preserve the S 4 -S 2 splitting. These deformed backgrounds will solve the equations of motion of type IIA supergravity (with sources) and will preserve at least two supercharges. We will argue below that some of these backgrounds are dual to Chern-Simons matter theories with fundamental massive flavors. The general ansatz for the ten-dimensional metric of our solutions in string frame takes the form: where the warp factor h and the functions f and g depend on the holographic coordinate r. Notice that f and g determine the sizes of the S 4 and S 2 within the internal manifold.
Actually, their difference f −g determines the squashing of the CP 3 and will play an important role in characterizing our solutions. We will measure this squashing by means of the function q, defined as: Clearly, the ABJM solution has q = 1. A departure from this value would signal a non-trivial deformation of the metric. Similarly, the RR two-and four-forms will be given by: where k is a constant and K = K(r), η = η(r) are new functions. The background is also endowed with a dilaton φ = φ(r). As compared with the ABJM value (2.15), the expression of F 2 in our ansatz contains the function η(r) which generically introduces an asymmetry between the S 4 and S 2 terms. Moreover, when η = 1 the two-form F 2 is no longer closed and the corresponding Bianchi indentity is violated. Indeed, one can check that: The violation of the Bianchi identity of F 2 means that we have D6-brane sources in our model. Indeed, since F 2 = * F 8 , if dF 2 = 0 then the Maxwell equation of F 8 contains a source term, which is due to the presence of D6-branes since the latter are electrically charged with respect to F 8 . The charge distribution of the D6-brane sources is determined by the function η, which we will call the profile function. The function K of the RR four-form can be related to the other functions of the ansatz by using its equation of motion d * F 4 = 0 and the flux quantization condition for the integral of * F 4 over the internal manifold. The result is [23]: where the integer N is identified with the ranks of the gauge groups in the gauge theory dual (i.e., with the number of colors). It is convenient to introduce a new radial variable x, related to r through the differential equation: From now on, all functions of the holographic variable are considered as functions of x, unless otherwise specified. The ten-dimensional metric in this new variable takes the form: It was shown in [23] that the background given by the ansatz written above preserves N = 1 supersymmetry in three dimensions if the functions satisfy a system of first-order differential equations. It turns out that this BPS system can be reduced to a unique secondorder differential equation for a particular combination of the functions of the ansatz. The details of this reduction are given in Appendix A. Here we will just present the final result of this analysis. First of all, let us define the function W (x) as: Then, the BPS system can be reduced to the following second-order non-linear differential equation for W (x): We will refer to (3.10) as the master equation and to W (x) as the master function. Interestingly, the BPS equations do not constrain the profile function η. Therefore, we can choose η(x) (which will fix the type of supersymmetric sources of our system) and afterwards we can solve (3.10) for W (x). Given η(x) and W (x) one can obtain the other functions that appear in the metric. Indeed, as proved in Appendix A, g(x) and f (x) are given by: while the warp factor h can be written as: where β is a constant that determines the behavior of h as x → ∞ (β = 0 if we impose that h → 0 as x → ∞). Finally, the dilaton is given by: From the expression of f and g in (3.11) it follows that the squashing function q can be written in terms of the master function W and its derivative as: . (3.14)

Classification of solutions
Let us study the behavior of the solutions of the master equation in the UV region x → ∞. This analysis will allow us to have a classification of the different solutions. We will assume that the profile function η(x) reaches a constant value as x → ∞, and we will denote: Let us restrict ourselves to the case in which η 0 = 0. We will assume that W (x) behaves for large x as: where A 0 and α are constants. It is easy to check that this type of behavior is consistent only when the exponent α ≥ 1 or, in other words, when W (x) grows at least as a linear function of x when x → ∞. We will also characterize the different solutions by the asymptotic value of the squashing function q, which determines the deformation of the internal manifold in the UV. Let us denote It follows from (3.14) that the asymptotic value of the squashing function and that of the profile function are closely related. Actually, this relation depends on whether the exponent α in (3.16) is strictly greater or equal to one. Indeed, plugging the asymptotic behavior (3.16) in (3.14) one immediately proves that: This result indicates that we have to study separately the cases α > 1 and α = 1. As we show in the next two subsections these two different asymptotics correspond to two qualitatively different types of solutions.

The asymptotic G 2 cone
Let us assume that the master function behaves as in (3.16) for some α > 1. By plugging this asymptotic form in the master equation (3.10) and keeping the leading terms as x → ∞, one readily verifies that the coefficient A 0 is not constrained and that the exponent α takes the value: Therefore, it follows from (3.18) that the asymptotic squashing is: Let us evaluate the asymptotic form of all the functions of the metric. From (3.11), we get, at leading order: where C is a constant of integration. Moreover, since q 0 = 2, the asymptotic value of the function f is: Let us now evaluate the warp factor h from (3.12). Clearly, we have to compute the integral: which vanishes when x → ∞. Therefore, by choosing the constant β in (3.12) to be nonvanishing we can neglect the integral (3.23) and, since e g W ′ → constant, then the warp factor h becomes also a constant when x → ∞. To clarify the nature of the asymptotic metric, let us change variables, from x to a new radial variable ρ, defined as ρ = 2C √ x. Then, after some constant rescalings of the coordinates the metric becomes: where ds 2 7 is: The metric (3.25) is a Ricci flat cone with G 2 holonomy, whose principal orbits at fixed ρ are CP 3 manifolds with a squashed Einstein metric. In the asymptotic region of large ρ the line element (3.25) coincides with the metric of the resolved Ricci flat cone found in [29], which was constructed from the bundle of self-dual two-forms over S 4 and is topologically S 4 × R 3 (see [30] for applications of this manifold to the study of the dynamics of M-theory).

The asymptotic AdS metric
Let us now explore the second possibility for the exponent α in (3.16), namely α = 1. In this case the coefficient A 0 cannot be arbitrary. Indeed, by analyzing the master equation as x → ∞ we find that A 0 and η 0 must be related as: On the other hand, A 0 should be related to the asymptotic squashing q 0 as in (3.18), which we now write as: By plugging (3.27) into (3.26) we arrive at the following quadratic relation between q 0 and η 0 : Using this equation we can re-express A 0 as: Moreover, we can solve (3.28) for q 0 and obtain the following two possible asymptotic squashings in terms of η 0 : Thus, there are two possible branches in this case, corresponding to the two signs in (3.30).
In this paper we will only consider the q + 0 case, since this is the one which has the same asymptotics as the ABJM solution when there are no D6-brane sources. Indeed, (3.30) gives q + 0 = 1 when η 0 = 1, which means that the internal manifold in the deep UV is just the un-squashed CP 3 (when η 0 = 1 there are no D6-brane sources in the UV, see (3.5)).
Let us now study in detail the asymptotic metric in the UV corresponding to the x → ∞ squashing q + 0 (which from now on we simply denote as q 0 ). By substituting η → η 0 and W → A 0 x in (3.11) and performing the integral, we get: where C is a constant of integration. Using (3.27) this expression can be rewritten as: where b is given by: The remaining functions of the metric can be found in a similar way. We get for f and h the following asymptotic expressions: Let us write the above expressions in terms of the original r variable, which can be related to x by integrating the equation: For large x we get: 36) and the functions g, f , and h can be written in terms of r as: where L 0 is given by: In terms of the asymptotic values η 0 and q 0 , L 0 can be written as: Using these results we find that the asymptotic metric takes the form: where we have rescaled the Minkowski coordinates as x µ → L 2 0 x µ . The metric (3.40) corresponds to the product of AdS 4 space with radius L 0 and a squashed CP 3 . The parameter b will play an important role in the following. Its interpretation is rather clear from (3.40): it represents the relative squashing of the CP 3 part of the asymptotic metric with respect to the AdS 4 part.
It is now straightforward to show that in the UV the dilaton reaches a constant value φ 0 , related to q 0 and η 0 as: while the RR four-form approaches the value: where Ω AdS 4 is the volume element of AdS 4 . Interestingly, when the profile function η is constant and equal to η 0 , the metric, dilaton, and forms written above solve the BPS equations not only in the UV, but also in the full domain of the holographic coordinate. Equivalently, W = A 0 x is an exact solution to the master equation (3.10) if η is constant and equal to η 0 and A 0 is given by (3.27). Actually, when η 0 = 1 one can check that q 0 = b = 1 and the asymptotic background becomes the ABJM solution (W = 2x for this case). Moreover, when η = η 0 > 1 the background corresponds 1 to the one found in [23] for the ABJM model with unquenched massless flavors, if one identifies η 0 with 1 + 3N f 4k , where N f is the number of flavors. The main objective of this paper is the construction of solutions which interpolate between the η = 1 ABJM background in the IR and the AdS 4 asymptotics with η 0 > 1 in the UV. Equivalently, we are looking for backgrounds such that the squashing function q(x) varies from the value q = 1 when x → 0 to q = q 0 > 1 for x → ∞. These backgrounds naturally correspond to gravity duals of Chern-Simons matter models with massive unquenched flavors. Indeed, in such models, when the energy scale is well below the quark mass the fundamentals are effectively integrated out and one should recover the unflavored ABJM model. On the contrary, if the energy scale is large enough the quarks can be taken to be massless and the corresponding gravity dual should match the one found in [23]. In the next section we present a one-parameter family of analytic unflavored solutions which coincide with the ABJM background in the deep IR and that have a squashing function q which grows as we move towards the UV. In Sections 5 and 6 we show that these running solutions can be used to construct the gravity duals to massive flavor that we are looking for.

The unflavored system
In this section we will consider the particular case in which the profile is η = 1. In this case dF 2 = 0 and there are no flavor sources. It turns out that one can find a particular analytic solution of the BPS system written in Appendix A. This solution was found in [23] in a power series expansion around the IR. Amazingly, this series can be summed exactly and a closed analytic form can be written for all functions. Let us first write them in the coordinate r. The functions f and g are given by: where c is a constant. For c = 0 this solution is ABJM without flavor (i.e., AdS 4 × CP 3 with fluxes), while for c = 0 it is a running background which reduces to ABJM in the IR, r → 0. The squashing function q can be immediately obtained from (4.1): 1 Notice that the expression for b written in (3.33) is equivalent to the one obtained in [23], namely: In order to check this equivalence it is convenient to use the following relation between q 0 and η 0 : q 0 + η 0 q 0 − η 0 = (q 0 + 1) (η 0 + q 0 )/4 .
For c = 0 the squashing function q interpolates between the ABJM value q = 1 in the IR and the UV value: The warp factor for this solution is: where α is a constant which has to be fixed by adjusting the behavior of the metric in the UV. Finally, the dilaton can be related to the warp factor as: (4.5) Let us now re-express this running analytic solution in terms of the variable x, related to r by (3.7), which in the present case becomes: This equation can be easily integrated: where γ is a constant of integration which parameterizes the freedom from passing to the x variable. By solving (4.7) for r we get: It is straightforward to write the functions f and g in terms of x: while the squashing function is: (4.10) The warp factor h in terms of the x variable is: By choosing appropriately the constant α in (4.11) this running solution behaves as the G 2 -cone in the UV region x → ∞. The dilaton as a function of x is: Working in the variable x, it is very interesting to find the function W (x). For the solution described above, W can be found by plugging the different functions in the definition (3.9). We find: One can readily check that the function written in (4.13) solves the master equation (3.10) for η = 1. For large x, the function W (x) behaves as: which corresponds to an exponent α = 3/2 in (3.16). This is consistent with the asymptotic value q 0 = 2 of the squashing found above.
Let us finally point out that we have checked explicitly that the geometry discussed in this section is free of curvature singularities.

Interpolating solutions
Let us now construct solutions to the BPS equations which interpolate between an IR region in which there are no D6-brane sources (i.e., with η = 1) and a UV region in which η > 1 and, therefore, the Bianchi identity of F 2 is violated. In the r variable the profile η(r) will be such that η(r) = 1 for r ≤ r q , while η(r) > 1 for r > r q . In the region r ≤ r q our interpolating solutions will reduce to the unflavored running solution of Section 4 for some value of the constant c. In order to match this solution with the one in the region r > r q it is convenient to work in the x coordinate (4.7). The point r = r q will correspond to some x = x q . Notice, however, that we have some freedom in performing the r → x change of variables. This freedom is parameterized by the constant γ in (4.7). We will fix this freedom by requiring that x q = 1, i.e., that the transition between the unflavored and flavored region takes place at the point x = 1. Then, (4.7) immediately implies that γ is given in terms of c and r q : We will use (5.1) to eliminate the constant c in favor of r q and γ. Actually, if we defineγ as: then c is given by In this running solution the squashing factor q is equal to one in the deep IR at x = 0. When x > 0 the function q(x) grows monotonically until it reaches a certain valueq at x = 1, which is related to the parameterγ as:γ In the region x ≥ 1 we have to solve the master equation (3.10) with η(x) > 1 and initial conditions given by the values of W and W ′ attained by the unflavored running solution at x = 1. These values depend on the parameterγ. They can be straightforwardly found by taking x = 1 in the function (4.13) and in its derivative. We find: Let us now write the different functions of these interpolating solutions in the two regions x ≤ 1 and x ≥ 1. For x ≤ 1 we have to rewrite (4.9) after eliminating the constant c by using (5.3) (which implies that γ/c = (γ + 1) r q /2). For the functions f , g, and the dilaton φ we get: where h is the function written in (4.11) for c = (γ −1)/(2 r q ). By using the general equations of Section 3, the solution for x ≥ 1 can be written in terms of W (x), which can be obtained by numerical integration of the master equation with initial conditions (5.5). This defines a solution in the full range of x for every γ and r q . Notice that (3.11), (3.12), and (3.13) contain arbitrary multiplicative constants, which we will fix by imposing continuity of f , g, and h at x = 1. We get for f , g and φ: The warp factor h(x) for x ≥ 1 is given by (3.12), where the integration constant β is related to the constant α of (4.11) by the following matching condition at x = 1: For a given profile function η(x), the solution described above depends on the parameter γ, which determines W (x) for x ≤ 1 through (4.13) and sets the initial conditions (5.5) needed to integrate the master equation in the x ≥ 1 region. The solution W (x) obtained numerically in this way grows generically as x 3/2 for large x which, according to our analysis in Section 3.1.1, gives rise to the geometry of the G 2 -cone in the UV. We are, actually, interested in obtaining solutions with the AdS asymptotics discussed in Section 3.1.2, for a set of profiles that correspond to flavor D6-branes with a non-zero quark mass. In order to get these geometries we have to fine-tune the parameter γ to some precise value which depends on the number of flavors. This analysis is presented in the next section.

Massive flavor
We now apply the formalism developed so far to find supergravity backgrounds representing massive flavors in ABJM. These solutions will depend on a deformation parameterǫ, related to the total number of flavors N f and the Chern-Simons level k as: where the factor 3/4 is introduced for convenience and N f /k is just N f N λ with λ = N/k being the 't Hooft coupling. The profile η, which corresponds to a set of smeared flavor D6-branes ending at r = r q , has been found in [23]. The main technique employed in [23] was the comparison between the smeared brane action for the distribution of flavor branes and the action for a fiducial embedding in a background of the type studied here. This fiducial embedding was determined by using kappa symmetry. With our present conventions, 2 assuming as above that r = r q corresponds to x = x q = 1, the function η(x) is given by: where Θ(x) is the Heaviside step function. It follows from (6.2) that the asymptotic value η 0 of the profile is: We want to find interpolating solutions of the type studied in Section 5 which have the AdS asymptotic behavior in the UV corresponding to the value of η 0 written in (6.3). These solutions have an asymptotic squashing (corresponding to q + 0 in (3.30)) given in terms ofǫ as: It follows from (6.4) that the asymptotic squashing q 0 grows with the deformation parameter ǫ. Indeed, q 0 = 1 forǫ = 0, whereas forǫ → ∞ the squashing reaches its maximum value: q 0 → 5/3. By using the relation between b and q 0 (eq. (3.33)) we also conclude that b → 5/4 whenǫ → ∞. To find the solution for x ≥ 1 we have to solve numerically the BPS system in this region. The most efficient way to proceed is by looking at the master equation for W (x) with the initial conditions (5.5). For a generic value of γ the numerical solution either gives rise to negative values of W (x) (which is unphysical for k > 0, see the definition (3.9)) or behaves in the UV as W (x) ∼ x 3 2 , which corresponds to the G 2 -cone asymptotics with q 0 = 2 discussed in Section 3.1.1. Only when γ is fine-tuned to some particular value (which depends onǫ) we get in the UV that W (x) ∼ x and that q 0 is given by (6.4). To determine this critical value of γ we have to perform a numerical shooting for every value ofǫ. In what follows we understand that γ = γ(ǫ) is the function of the deformation parameter which results of this shooting. The function γ(ǫ) is plotted in Fig. 1, where we notice that γ(ǫ = 0) = 0 and, therefore, we recover the unflavored ABJM background when the deformation parameter vanishes. In the opposite limitǫ → ∞, the function γ(ǫ) grows asǫ 2 . Actually, γ(ǫ) can be accurately represented by a function of the type: From the function W (x) we can obtain f and g by performing the integrals in (5.7). The whole metric is determined if h is known. We will compute h from (3.12) with β = 0, which corresponds to requiring that h → 0 in the UV. In the x ≤ 1 region the warp factor h(x) is given by (4.11), with the constant α determined by the matching condition (5.8). The limit on the left-hand side of (5.8) can be determined explicitly from (4.11): Notice that r q (the value of the r coordinate at the tip of the flavor branes) appears as a free parameter in eqs. (5.6) and (5.7). Actually, r q can be easily related to the mass m q of the quarks which deform the geometry. Indeed, by computing the Nambu-Goto action of a fundamental string stretching along the holographic direction between r = 0 and r = r q at fixed value of all the other spacelike coordinates in the geometry (3.1), we get that m q and r q are linearly related as: where α ′ is the Regge slope (which we will take to be equal to one in most of our equations). Equivalently, we can relate m q to the constant c appearing in the solution in the x ≤ 1 region: where γ(ǫ) is the function obtained by the shooting and only depends on the deformation parameterǫ.
We have computed the curvature invariants for the flavored metric and we have checked that the geometry is regular both in the IR (x → 0) and UV (x → ∞). However, the curvature has a finite discontinuity at x = 1, as can be directly concluded by inspecting Einstein's equations (see Appendix A). This "threshold" singularity occurs at the point where the sources are added and could be avoided by smoothing the introduction of brane sources with an additional smearing (see the last article in [19] for a similar analysis in other background).

UV asymptotics
The full background in the x ≥ 1 region must be found by numerical integration and shooting, as described above. However, in the UV region x → ∞ one can solve the master equation (3.10) in power series for large x. Indeed, one can find a solution where W (x) is represented as: where the coefficients A 2i can be obtained recursively. The coefficient A 0 of the leading term was written in (3.27). The next two coefficients are: Notice that a linear behavior of W (x) with x corresponds to a conformal AdS background, whereas the deviations from conformality are encoded in the non-linear corrections. From the result written above one can immediately obtain the asymptotic behavior of the squashing function for large x. Indeed, let us use in the expression of q in terms of W and W ′ (eq. (3.14)) the following large x expansion: We get: where q 0 is the asymptotic value of the squashing (see (6.4)) and q 2 is given by: where b is related to q 0 andǫ by (3.33) and (6.4). Similarly, we can find analytically the first corrections to the UV conformal behavior. The details of these calculations are given in Appendix B. In this section we just present the final results. First of all, let us define the constant κ (depending on the deformationǫ) as: Then the functions g and f can be expanded for large x as: where the coefficients g 2 and f 2 are: Moreover, the UV expansion of the warp factor h and the dilaton is: with the coefficients h 2 and φ 2 given by: It is also interesting to write the previous expansions in terms of the r variable. Again, the details are worked out in Appendix B and the final result is: where the coefficientsg 2 ,f 2 ,h 2 , andφ 2 are related to the ones in (6.17) and (6.19) as: Recalling (see (6.8)) that r q = 2π m q , it is clear from (6.20) that the deviation from con- formality is controlled by the quark mass and that the parameter b determines the power of the first mass corrections. In our holographic context this is quite natural if one takes into account that b determines the dimension ∆ of the quark-antiquark bilinear operator in the theory with dynamical quarks (∆ = 3 − b, see [16,23] and below). The coefficients of these mass corrections depend on the constants g 2 , f 2 , h 2 , and φ 2 (whose analytic expressions we know from eqs. (6.17) and (6.19)), as well as on the constant κ, which must be determined numerically. κ as a function of the deformation parameter is plotted in Fig. 3. From this plot we notice that κ(ǫ) interpolates continuously between κ = 1 forǫ = 0 and some positive constant value at largeǫ.

Holographic entanglement entropy
In a quantum theory the entanglement entropy S A between a spatial region A and its complement is defined as the entropy seen by an observer in A which has no access to the degrees of freedom living in the complement of A. It can be computed as the von Neumann entropy for the reduced density matrix obtained by taking the trace over the degrees of freedom of the complement of A. For quantum field theories admitting a gravity dual, Ryu and Takayanagi proposed in [25] a simple prescription to compute S A from the corresponding supergravity background. The holographic entanglement entropy between A and its complement in the proposal of [25] is obtained by finding the eight-dimensional spatial surface Σ whose boundary coincides with the boundary of A and is such that it minimizes the functional: where the ξ's are a system of eight coordinates of Σ, G 10 is the ten-dimensional Newton constant (G 10 = 8π 6 in our units) and g 8 is the induced metric on Σ in the string frame. The functional S A evaluated on the minimal surface Σ is precisely the entanglement entropy between the region A and its complement. In our case A is a region of the (x 1 , x 2 )-plane. In this section we will study in detail the case in which the region A is a disk with radius R as depicted in Fig. 4 (see Appendix C for the analysis of the entanglement entropy of a strip in the (x 1 , x 2 )-plane). In order to deal with the disk case it is convenient to choose a system of polar coordinates for the plane: We will describe the eight-dimensional fixed time surface Σ by a function ρ = ρ(x) with ρ being the radial coordinate of the boundary plane and x the holographic coordinate of the bulk. The eight-dimensional induced metric is: where ρ ′ denotes the derivative with respect to the holographic coordinate x and the function G(x) is defined as: Let us next define a new function H(x) as: Then, the entanglement entropy as a function of R is given by: where V 6 = 32 π 3 /3 is the volume of the internal manifold and x * is the x coordinate of the turning point of Σ. The Euler-Lagrange equation of motion derived from the entropy functional (7.6) is: Notice that the integrand in (7.6) depends on the independent variable x and we therefore cannot find a first-integral for the second-order differential equation (7.7). Thus, we have to deal directly with (7.7), which must be solved with the following boundary conditions at the tip of Σ: Notice also that the radius R of the disk at the boundary is just the UV limit of ρ: The integral (7.6) for S(R) diverges due to the contribution of the UV region of large x . In order to characterize this divergence and to extract the finite part, let us study the behavior of the integrand in (7.6) as x → ∞. From the definitions of the functions H(x) and G(x) and the UV behavior written in (6.16) and (6.18), it follows that H and G display a power-like behavior as x → ∞, where the coefficients H ∞ and G ∞ are By taking the x → ∞ values of ρ and ρ ′ (ρ = R and ρ ′ = 0) inside the integral in (7.6), as well as the asymptotic form of H(x) and G(x) (eq. (7.10)), we get: where x max is the maximum value of the holographic coordinate x (which acts as a UV regulator). After performing the integral, we obtain: Let us rewrite (7.13) in terms of physically relevant quantities. First of all, we notice that: where F U V (S 3 ) is the free energy 3 of the massless flavored theory on the three-sphere: When the field theory is formulated on a three-sphere, its free energy is defined as: where Z S 3 is the Euclidean path integral. For a CFT whose gravity dual is of the form AdS 4 × M 6 , the holographic calculation of F (S 3 ) gives [31]: where L is the AdS 4 radius and G N is the effective four-dimensional Newton's constant.
where the function ξ N f k is given by: (7.18) In (7.18) η 0 = 1 +ǫ and q 0 is written in (6.4) in terms of the deformation parameter. For the unflavored ABJM theory the free energy is given by (7.17) with ξ = 1. This formula displays the famous N 3 2 behavior. The function ξ(N f /k) encodes the corrections to this behavior due to the smeared massless flavors. It was first computed in [23], where it was shown that it is remarkably close to the value found in [11] for localized embeddings. The function ξ is a monotonic function of the deformation parameter which grows as √ǫ for large values ofǫ. Using (7.14) and the fact that, in the deep UV region of large x, r max = κ r q x 1 b max (see (B.22)), we can rewrite S div (R) as: We notice in (7.19) that S div (R) diverges linearly with r max . The coefficient of this divergent term is linear in the disk radius R and in F U V (S 3 ). The latter is a measure of the effective number of degrees of freedom of the flavored theory in the high-energy UV limit in which the flavors can be considered to be massless. The appearance of The separation between the divergent and finite parts of S(R) has ambiguities. In order to solve these ambiguities, Liu and Mezei proposed in [26] to consider the function F (R), defined as: It was argued in [26] that F (R) is finite and a monotonic function of R which provides a measure of the number of degrees of freedom of a system at a scale R.
In a 3d CFT the entanglement entropy for a disk of radius R has the form: where α is a UV divergent non-universal part and β is finite and independent of R. It was shown in [32] that the finite part β is equal to the free energy of the theory on S 3 . Notice that F = β when S is of the form (7.21). Therefore, for a conformal fixed point the function F is constant and equal to the free energy on the three-sphere of the corresponding CFT. In the next subsection we will obtain the UV and IR values of F and we will show that they coincide with the free energies on S 3 of the massless flavored theory and of the unflavored ABJM model, respectively.
It is interesting to point out that the entanglement entropy of a disk in a (2+1)-dimensional system at large R can also be written in the form (7.21), if one neglects terms which vanish in the R → ∞ limit. In a gapped system, the R-independent part β of the right-hand side of (7.21) is the so-called topological entanglement entropy [34,35] and serves to characterize topologically ordered many-body states which contain non-local entanglement due to nonlocal correlations (examples of such states are the Laughlin states of the fractional quantum Hall effect or the Z 2 fractionalized states). The topological entanglement entropy β is related to the so-called total quantum dimension D of the system as β = log D. In general D > 1 (or β > 0) signals a topological order (for example D = √ q for the quantum Hall system with filling fraction ν = 1/q, with q an odd integer).  Figure 5: Plot of the running free energy F as a function of r q R forǫ = 0 (bottom, blue), 1 (middle, red), and 9 (top, brown).
For our system, we can obtain the embedding function ρ(x) by numerical integration of the differential equation (7.7) and then we can get the functions S(R) and F (R) by using (7.6) and the definition of F (R) in (7.20). The results for the latter are plotted as a function of R r q ∝ R m q in Fig. 5. We notice that F (R) is a monotonically decreasing function that interpolates smoothly between the two limiting values at R m q = 0 and R m q → ∞. The UV limit of F at small R m q equals the free energy (7.17) of the massless flavored theory, whereas for large R m q the function F approaches the free energy of the unflavored ABJM model (i.e., the value in (7.17) with ξ = 1). This behavior is in agreement with the general expectation in [26,28] and corresponds to a smooth decoupling of the massive flavors as their mass m q is increased continuously.
We will study the UV and IR limits of F analytically in the next two subsections. Some details of these calculations are deferred to Appendix C, where we also study the entanglement entropy for the strip geometry.

UV limit
In order to study the UV limit of the disk entanglement entropy, let us write the Euler-Lagrange equation (7.7) when H(x) and G(x) are given by their asymptotic values (7.10): This equation can be solved exactly by the function: which clearly satisfies the initial conditions (7.8), with the following value of the turning point coordinate x * : Therefore, the turning point coordinate x * is large if r q or R are small. In this case it would be justified to use the asymptotic UV values of the functions G and H, since the minimal surface Σ lies entirely in the large x region. Notice also that (7.23) can be written in terms of x * as: In order to calculate the entropy in this UV limit it is very useful to use the following relation satisfied by the function ρ U V (x) written in (7.23): Making use of (7.26) in (7.6), we find the following expression for the entanglement entropy: The divergent part of this integral is due to its upper limit and is just given by (7.19). The finite part of S is: where, in the second step, we used (7.24) to eliminate x * . Notice that the right-hand side of (7.28) is independent of the disk radius R. Moreover, by using (7.11) and (7.14) we find that: Therefore, in this UV limit, the dependence on R of the entanglement entropy takes the form (7.21), where β is just the free energy of the massless flavored theory on the three-sphere. It follows trivially from this form of S U V and the definition (7.20) that F U V = β and therefore: It is also possible to compute analytically the first correction to (7.30) for small values of r q R. The details of this calculation are given in Appendix C. Here we will just present the final result, which can be written as: where c U V is a constant coefficient depending on the deformation parameter (see eq. (C.28)) and ∆ U V = 3 − b is the dimension of the quark-antiquark bilinear operator in the UV flavored theory (this dimension was found in Section 7.3 of [23] from the analysis of the fluctuations of the flavor branes, see also [16]). It is interesting to point out that (7.31) is the behavior expected [26] for a flow caused by a source deformation with a relevant operator of dimension ∆ U V . Moreover, one can verify that c U V is negative for all values of the deformation parameterǫ, which confirms that the UV fixed point is a local maximum of F .

IR limit
Let us now analyze the IR limit of the entanglement entropy S(R) and of the function F (R). This limit occurs when the 8d surface Σ penetrates deeply into the geometry and, therefore, when the coordinate x * of the turning point is small (x * ≪ 1). This happens when either the disk radius R or the quark mass m q = r q /2π are large. Looking at the embedding function ρ(x) obtained by numerical integration of (7.7) one notices that, when x * is small, the function ρ(x) is approximately constant and equal to its asymptotic value ρ = R in the region x ≥ 1 (see Fig. 6). Therefore, the dependence of ρ on the holographic coordinate x is determined by the integral of (7.7) in the region x ≤ 1, where the background is given by the unflavored running solution. Actually, when x * is small it is a good approximation to consider (7.7) for the unflavored background, i.e., when the constants c, γ → 0, with γ/c fixed and given by γ/c = r q . In this limit the different functions of the background are: It follows that G(x) and H(x), as defined in (7.4) and (7.5), are then given by: 33) where, in the last step, we have defined the constants G 0 and H 0 , and φ IR is the constant dilaton corresponding to the unflavored background, Let us now elaborate on the expression (7.6) for the entanglement entropy. We split the integration interval of the variable x as [x * , ∞] = [x * , 1] ∪ [1, ∞] and take into account that one can put ρ = R = constant in the region x ≥ 1. We get: The second term in (7.35) is linear in R and will not contribute to F (R). To evaluate the first integral in (7.35) we must determine the embedding function ρ(x) by integrating (7.7) when G(x) and H(x) are given by their IR values (7.33). The resulting equation is just the same as (7.22) with G ∞ and H ∞ substituted by G 0 and H 0 and b = 1. Then, the function ρ(x) can be written as in (7.23), whereR is a constant. By requiring that ρ IR (x = 1) = R, we get: It follows from (7.36) that the coordinate x * of the turning point is given by: Notice that when r q → ∞ (and G 0 → 0) or R is large one can neglect the G 0 in the denominator of (7.38) and then x 2 which is a small number. Moreover, by using the explicit form (7.36) of ρ(x) in this IR region, we get: 39) and the first integral in (7.35) can be explicitly evaluated: Notice that, at leading orderR ≈ R and, thus, the first term in (7.40) does not contribute to F (R). Then, the IR limit of F is determined by the second contribution in (7.40) and given by: Moreover, from the values of G 0 and H 0 written in (7.33) we get: where F IR (S 3 ) is the free energy on the three-sphere of the unflavored ABJM theory: It follows that the IR limit of the F function is: as expected in the deep IR limit in which the flavors become infinitely massive and can therefore be integrated out. The corrections to the result (7.44) near the IR fixed point could be obtained by applying the techniques recently introduced in [33]. We will not attempt to perform this calculation here.

Wilson loops and the quark-antiquark potential
In this section we evaluate the expectation values of the Wilson loop and the corresponding quark-antiquark potential for our model. We will employ the standard holographic prescription of refs. [36,37], in which one considers a fundamental string hanging from the UV boundary. Then, one computes the regularized Nambu-Goto action for this configuration, from which the qq potential energy can be extracted. In a theory with dynamical flavors this potential energy contains information about the screening of external charges by the virtual quarks popping out from the vacuum. In our case we expect having a non-trivial flow connecting two conformal behaviors as we move from the UV regime of small qq separation (in units of the quark mass m q ) to the IR regime of large qq distance. We will verify below that this expectation is indeed fulfilled by our model. Let us denote by (t, x 1 , x 2 ) the Minkowski coordinates and consider a fundamental string for which we take (t, x 1 ) as its worldvolume coordinates. If the embedding is characterized by a function x = x(x 1 ), with x being the holographic coordinate, the induced metric is: where x ′ denotes the derivative of x with respect to x 1 . The Nambu-Goto Lagrangian density takes the form: As L N G does not depend on x 1 , we have the following conservation law: Therefore, if x * denotes the turning point of the string, we have the first integral of the equations of motion: . Then x ′ is given by: where the two signs correspond to the two branches of the hanging string. The qq separation d in the x 1 direction is: In order to compute the potential energy of the qq pair, let us evaluate the on-shell action. By using the first integral (8.4) it is straightforward to check that the on-shell value of the Nambu-Goto Lagrangian density is: Therefore, the on-shell action becomes: where T = dt. The integral (8.8) is divergent and must be regularized as in [36,37] by subtracting the action of two straight strings stretched between the origin and the UV boundary, which corresponds to subtracting the (infinite) quark masses in the static limit.
After applying this procedure we arrive at the following expression for the regulated on-shell action: from which we get the qq potential energy: where r * is the r coordinate of the turning point: From (8.6) and (8.10) we have computed numerically the potential energy E qq as a function of the qq distance d. The result of this numerical calculation is shown in Fig. 7. As mentioned above, we expect to have a potential energy which interpolates between the two conformal behaviors with E qq ∝ 1/d at the UV and IR. Actually, in the limiting cases in which r q d is small or large the function E qq (d) can be calculated analytically (see Appendix D). In both cases E qq ∝ 1/d, but with different coefficients. Indeed, in the UV limit r q d → 0, the qq potential can be approximated as: 4 where λ = N/k is the 't Hooft coupling and σ is the so-called screening function: Notice that σ encodes all the dependence of the right-hand side of (8.12) on the number of flavors. Actually, the potential (8.12) is just the one corresponding to having massless flavors (which was first computed for this model in [23]), as expected in the high-energy UV regime in which all masses can be effectively neglected. The function σ characterizes the corrections of the static qq potential due to the screening produced by the unquenched massless flavors (σ → 1 for N f → 0, whereas σ decreases as σ ∝ k/N f for N f large). In Fig. 7 we compare the leading UV result (8.12) with the numerical calculation in the small r q d region.
Similarly, one can compute analytically the qq potential in the region where r q d is large. At leading order the result is (see Appendix D): In Fig. 7 we compare the analytic expression (8.14) to the numerical result in the large distance region. Notice that the difference between (8.12) and (8.14) is that the screening function σ is absent in (8.14). Therefore, in the deep IR the flavor effects on the qq potential disappear, which is consistent with the intuition that massive flavors are integrated out at low energies.

Two-point functions of high dimension operators
In this section we study the two-point functions of bulk operators with high dimension. The form of these correlators can be obtained semiclassically by analyzing the geodesics of massive particles in the dual geometry [38][39][40], where m is the mass of the bulk field dual to O. We are assuming that m is large in order to apply a saddle point approximation in the calculation of the correlator. In (9.1) L r (x, y) is a regularized length along a spacetime geodesic connecting the boundary points x and y.
To find these geodesics, let us write the Einstein frame metric of our geometry as: Then, the induced metric for a curve parametrized as x = x(x 1 ) is: with x ′ = dx/dx 1 and G(x) is the function defined in (7.4). Therefore, the action of a particle of mass m whose worldline is the curve x = x(x 1 ) is: The geodesics we are looking for are solutions of the Euler-Lagrange equation derived from the action (9.4). This equation has a first integral which is given by: where φ * ≡ φ(x = x * ) and h * ≡ h(x = x * ), with x * being the x coordinate of the turning point, i.e., the minimum value of x along the geodesic. It follows from (9.5) that: The spatial separation l of the two points in the correlator can be obtained by integrating 1/x ′ . We get: Moreover, the length of the geodesic can be obtained by integrating ds 1 over the worldline, This integral is divergent. In order to regularize it, let us study the UV behavior of the integrand. For large x, the functions h(x) and G(x) behave as in (6.18) and (7.10), respectively. Thus, at leading order for large x, In the UV region x → ∞, the integrand in L behaves approximately as 2 b −1 L 0 e −φ 0 /4 x −1 , which produces a logarithmic UV divergence when it is integrated. In order to tackle this divergence, let us regulate the integral by extending it up to some cutoff x max and renormalize the geodesic length by subtracting the divergent part. Accordingly, we define the renormalized geodesic length as: where C is a constant to be fixed by choosing a suitable normalization condition for the correlator.
Our background interpolates between two limiting AdS 4 geometries, at the UV and IR, with different radii. For an equal-time two-point function O(t, l) O(t, 0) the UV and IR limits should correspond to the cases in which r q l is small or large, respectively. At the two endpoints of the flow, the theory is conformal invariant and the two-point correlator behaves as a power law in l. We can use this fact to fix the normalization constant C in (9.10). Actually, we will assume that the field O is canonically normalized in the short-distance r q l → 0 limit and, therefore, the UV limit of the two-point correlator is: where the √ λ and r q factors have been introduced for convenience. In (9.11) ∆ U V is the conformal dimension of the operator O in the UV CFT, which for the dual of a bulk field of mass m is: where we have taken into account that m is large and that L 0 e − φ 0 4 is the AdS 4 radius of the UV massless flavored geometry in the Einstein frame. It is shown in Appendix E that, indeed, the correlators derived from (9.10) display the canonical form (9.11) if the constant C is chosen appropriately (see (E.8)). In Appendix E we have also computed the first deviation from the conformal UV behavior. In this case the numerator on the right-hand side of (9.11) is not one but a function f ∆ (r q l/ √ λ) such that f ∆ (r q l/ √ λ = 0) = 1. We show in Appendix E that f ∆ (r q l/ √ λ) − 1 ∝ (r q l/ √ λ) 2b for small r q l. The explicit form of the first correction to the non-conformal behavior can be computed analytically from the mass corrections of Section 6.1 and Appendix B (see eqs. (E.19)-(E.21)).  Figure 8: Plot of m L r / 2∆ IR log(r q l/ √ λ) versus the logarithm of r q l/ √ λ. Notice that, according to (9.1), m L r = − log O(t, l) O(t, 0) . The three curves correspond toǫ = 0 (bottom, blue),ǫ = 1 (middle, red), andǫ = 9 (top, brown). In the deep UV (r q l → 0) the curves approach the constant value ∆ U V /∆ IR . The dashed curves correspond to the behavior (9.14), with the normalization constant N given in (E.35).
When the distance r q l is large the theory reaches a new conformal point. Accordingly, the two-point function should behave again as a power law. Notice, however, that the conformal dimension ∆ IR in the IR of an operator dual to a particle of mass m is different from the UV value (9.12). Indeed, in the IR the conformal dimension ∆ IR for an operator O of mass m is the one corresponding to the unflavored ABJM theory, where L ABJM and φ ABJM are given, respectively, in (2.3) and (2.4). Actually, one can check that ∆ U V ≥ ∆ IR and that ∆ U V /∆ IR ∝ǫ 1 16 for large values of the deformation parameterǫ.
The calculation of the two-point function in the IR limit of large r q l is performed in detail in Appendix E, with the result: where N is a constant whose analytic expression is written in (E. 35). Notice that N = 1 due to our choice of the constant C in (9.10), which corresponds to imposing the canonical normalization (9.11) to the two-point function in the UV regime. The results obtained by the numerical evaluation of the integral (9.10) interpolate smoothly between the conformal behaviors (9.11) and (9.14). This is shown in Fig. 8 ) as a function of log(r q l/ √ λ). For small values of r q l/ √ λ the curve asymptotes to the ratio ∆ U V /∆ IR of conformal dimensions, in agreement with (9.11), whereas for large r q l/ √ λ we recover the IR behavior (9.14).

Meson spectrum
Let us now test the flow encoded in our geometry by analyzing the mass spectrum of qq bound states. We will loosely refer to these bound states as mesons, although our background is not confining and quarkonia would be a more appropriate name for them. To carry out our analysis we will introduce additional external quarks, with a mass µ q not necessarily equal to the mass m q of the quarks which backreact on the geometry. To distinguish between the two types of flavors we will call valence quarks to the additional ones, whereas the unquenched N f dynamical flavors of the geometry will be referred to as sea quarks. The ratio µ q /m q of the masses of the two types of quarks will be an important quantity in what follows. Indeed, µ q /m q is the natural parameter for the holographic renormalization group trajectory. When µ q /m q is large (small) we expect to reach a UV (IR) conformal fixed point, whereas for intermediate values of this mass ratio the theory should flow in such a way that the screening effects produced by the sea quarks decrease as we move towards the IR. Within the context of the gauge/gravity duality, the valence quarks can be introduced by adding an additional flavor D6-brane, which will be treated as a probe in the backreacted geometry. The mesonic mass spectrum can be obtained from the normalizable fluctuations of the D6-brane probe. The way in which the D6-brane probe is embedded in the tendimensional geometry preserving the supersymmetry of the background can be determined by using kappa symmetry. For fixed values of the Minkowski and holographic coordinates, the D6-brane extends over a cycle inside the CP 3 which has two directions along the S 4 base and one direction along the S 2 fiber. In order to specify further this configuration, let us parameterize the SU(2) left invariant one-forms ω i of the four-sphere metric (2.7) in terms of three anglesθ,φ andψ, with 0 ≤θ ≤ π, 0 ≤φ < 2π, 0 ≤ψ ≤ 4π. Then, our D6-brane probe will be extended along the Minkowski directions and embedded in the geometry in such a way that the anglesθ and ϕ are constant and that the angle θ of the S 2 fiber depends on the holographic variable x. The pullbacks (denoted by a hat) of the left-invariant SU(2) one-forms (10.1) areω 1 =ω 2 = 0 andω 3 = dψ. The kappa symmetric configurations are those for which the function θ(x) satisfies the first order BPS equation [23]: which can be integrated as: Here x * is the minimum value of the variable x for the embedding, i.e., the value of x for the tip of the brane. This minimum value of the coordinate x for the embedding is related to the mass µ q of the valence quarks introduced by the flavor probe. Indeed, by computing the Nambu-Goto action of a fundamental string stretched in the holographic direction between x = 0 and x = x * we obtain µ q as: x dx . (10.4) In the following we will take the Regge slope α ′ = 1. Moreover, to simplify the description of the embedding, let us introduce the angular coordinate α, defined as follows: and let us define new angles β and ψ as: where ϕ is the angle in (2.9). One can check that the ranges of the new angular variables are 0 ≤ α < π, 0 ≤ β , ψ < 2π. We will take the following set of worldvolume coordinates for the D6-brane: ζ a = (x µ , x, α, β, ψ) .
Then, it is straightforward to verify that the induced metric on the D6-brane worldvolume takes the form: (10.8) We will restrict ourselves to study a particular set of fluctuations of the D6-brane probe, namely the fluctuations of the worldvolume gauge field A a . The equation for these fluctuations is: ∂ a e −φ − det g g ac g bd F cd = 0 , (10.9) where g ab is the induced metric (10.8). More concretely, we will study this equation for the following ansatz for A a : where ξ µ is a constant polarization vector and A i denote the components along the angular directions. These modes are dual to the vector mesons of the theory, with k µ being the momentum of the meson (η µν k µ k ν = −m 2 , with m being the mass of the meson). The non-vanishing components of the field strength F ab are: Moreover, by taking b = µ in (10.9) we arrive at the following differential equation for the radial function R: The mass levels correspond to the values of m for which there are normalizable solutions of (10.13). They can be obtained numerically by the shooting technique. One gets in this way a discrete spectrum depending on a quantization number n (n ∈ Z, n ≥ 0). The numerical results for the first three levels are shown in Fig. 9 as functions of the mass ratio µ q /m q . One notices in these results that the meson masses increase as we move from the IR (µ q /m q → 0) to the UV (µ q /m q → ∞). This non-trivial flow is due to the vacuum polarization effects of the sea quarks, which are enhanced as we move towards the UV and the sea quarks become effectively massless. This is the expected behavior of bound state  Figure 9: Numerical values of the meson masses for the first three levels (n = 0, 1, 2) as a function of the sea quark mass m q for deformation parameterǫ = 9. The solid curves depict the WKB estimate (10.14).
masses for a theory in the Coulomb phase, since the screening effects reduce effectively the strength of the quark-antiquark force. One can get a very accurate description of the flow by applying the WKB approximation. The detailed calculation is presented in Appendix F. The WKB formula for the mass spectrum is: 2 ξ(x * ) (n + 1)(2n + 1) , (n = 0, 1, 2, · · · ) , (10.14) where ξ(x * ) is the following integral: The WKB mass levels (10.14) are compared with those obtained by the shooting technique in Fig. 9. We notice from these plots that the estimate (10.14) describes rather well the numerical results along the flow. Moreover, we can use the UV and IR limits of the functions g and h to obtain the asymptotic form of the WKB spectrum at the endpoints of the flow. This analysis is performed in detail in Appendix F (see eqs. (F.14) and (F.20)). As expected, in the deep IR the mass levels coincide with those of the unflavored ABJM model. In this latter model the mass spectrum of vector mesons can be computed analytically since the fluctuation equation can be solved in terms of hypergeometric functions [12,13]. When µ q /m q → ∞ the meson masses coincide with those obtained for the massless flavored model of [23]. We can use the WKB formulas (F.14) and (F. 20) for the spectrum at the endpoints of the renormalization group flow to estimate the variation generated in the meson masses by changing the sea quark mass m q and switching on and off gradually the screening effects. It is interesting to point out that, within the WKB approximation, the ratio of these masses only depends on the number of flavors, and is given by: where σ is the screening function defined in (8.13). As expected, when N f = 0 the right-hand side of (10.16) is equal to one, i.e., there is no variation of the masses along the flow. On the contrary, when N f > 0 the UV/IR mass ratio in (10.16) is always greater than one, which means that the masses grow as we move towards the UV and the screening effects become more important. In Appendix F we have expanded (10.16) for low values of the deformation parameterǫ (see (F.21)). Moreover, for largeǫ the UV/IR mass ratio grows as √ǫ (see (F.22) for the explicit formula).

Summary and conclusions
In this paper we obtained a holographic dual to Chern-Simons matter theory with unquenched flavor in the strongly-coupled Veneziano limit. The flavor degrees of freedom were added by a set of D6-branes smeared along the internal directions, which backreacted on the geometry by squashing it, while preserving N = 1 supersymmetry. We considered massive flavors and found a non-trivial holographic renormalization group flow connecting two scale-invariant fixed points: the unflavored ABJM theory at the IR and the massless flavored model at the UV.
The quark mass m q played an important role as a control parameter of the solution. By increasing m q our solutions became closer to the unflavored ABJM model and we smoothly connected the unquenched flavored model to the ABJM theory without fundamentals. After this soft introduction of flavor no pathological behavior was found. Indeed, our backgrounds had good IR and UV behaviors, contrary to what happens to other models with unquenched flavor [22]. This made the ABJM model especially adequate to analyze the effects of unquenched fundamental matter in a holographic setup.
We analyzed different flavor effects in our model. In general, the screening effects due to loops of fundamentals were controlled by the relative value of the quark mass m q with respect to the characteristic length scale l of the observable. If m q l was small, which corresponds to the UV regime, the flavor effects were important, whereas they were suppressed if m q l is large, i.e., at the IR. Among the different observables that we analyzed, the holographic entanglement entropy for a disk was specially appropriate since it counts precisely the effective number of degrees of freedom which are relevant at the length scale given by the radius of the disk. By using the refined entanglement entropy F introduced in [26], we explicitly obtained the running of F and verified the reduction of the number of degrees at the IR that was mentioned above. The other observables studied also supported this picture.
We end this paper with a short discussion on the outlook. We are convinced that our model could serve as a starting point to gain new insights on the effects of unquenched flavor in other holographic setups. One possible generalization could be the construction of a black hole for the unquenched massive flavor. Such a background could serve to study the meson melting phase transition which occurs when the tip of the brane approaches the horizon. This system was studied in [16], in the case in which the massive flavors are quenched and the corresponding flavor brane is a probe. Another possibility would be trying to find a gravity dual of a theory in which the sum of the two Chern-Simons levels is non-vanishing. According to [41] we should find a flavored solution of type IIA supergravity with non-zero Romans mass.
Our program is also to converge toward increasingly realistic holographic condensed matter models capable of testable predictions. To make contact with any condensed matter system, one is forced to consider non-vanishing components of the gauge field in the background or at the probe level. A natural flow of ideas taking one to land in the former case typically requires a deep understanding of the probe brane dynamics with worldvolume gauge fields turned on. As an initial step in this direction, we have started exploring what physical phenomena we will encompass by turning on a charge density, magnetic field, and internal flux on the worldvolume of an additional probe D6-brane. The variety of different phenomena seems incredibly rich, much in parallel with recent works on the D3-D7' system [42] and the closely related D2-D8' system [43]. Our findings will be reported elsewhere. FPU12/00481. N. J. is supported also through the Juan de la Cierva program. N. J. and A. V. R. wish to thank Centro de Ciencias de Benasque Pedro Pascual and N. J. the Kavli IPMU for warm hospitalities while this work was in progress.

A BPS equations
In this Appendix we will derive the master equation (3.10), as well as the equations that allow to construct the metric and dilaton from the master function W (x) (i.e., (3.11), (3.12), and (3.13)).
Let us begin by writing the BPS equations that guarantee the preservation of N = 1 SUSY. They can be written in terms of the function Λ introduced in [23], which is defined as the following combination of the dilaton and the warp factor: Then, it was proved in [23] that Λ(r), f (r), and g(r) are solutions to the following system of first-order differential equations: Moreover, the warp factor h(r) can be recovered from Λ, f , and g through: where α is an integration constant. Given h and Λ, the dilaton φ is obtained from (A.1). The function K of the RR four-form can be related to the other functions of the background by using (3.6). Alternatively, K can be obtained from the BPS system as: In terms of the x variable defined in (3.7), the BPS system (A.2) becomes: In order to reduce this system, let us define as in [23] the functions Σ(x) and ∆(x), Then, one can easily show that Σ(x) and ∆(x) satisfy the system: whereas g can be obtained from Σ and ∆ by integrating the equation: Let us next define the master function W (x) as in (3.9). One immediately verifies that, in terms of the functions ∆ and Σ, this definition is equivalent to By computing the derivative of (A.9) and using the BPS system (A.7), one can easily prove that: From (A.10) one immediately finds: where the prime denotes derivative with respect to x. Moreover, from the BPS system we can calculate the derivative of Σ + ∆ and write the result as: Plugging (A.11) into (A.12), we arrive at the following second-order equation for W (x): x d dx which can be straightforwardly shown to be equivalent to the master equation (3.10).
Let us see now how one can reconstruct the full solution from the knowledge of the function W (x). First of all, we notice that from the expression of W in (A.9), we get: (A.14) By combining this expression with (A.11) we obtain Σ and ∆: .
(A. 15) By noticing that e 2∆ = q we arrive at the representation of the squashing function q written in (3.14). Moreover, by using this result in (A.8), we obtain the differential equation satisfied by g: which allows to obtain g(x) once W (x) is known. The result of this integration is just the expression written in (3.11). Moreover, taking into account the expression of the squashing factor q we get precisely the expression of f written in (3.11).
Let us now compute Λ by using Λ = Σ + ∆ + g and (A.11). We get: and, after using (3.11), we arrive at: By using this result and (3.11) in (A.3), we get that the warp factor can be written as in (3.12). The expression (3.13) for the dilaton is just a consequence of the definition of Λ in (A.1) and of (A.18).

A.1 Equations of motion
Let us now verify that the first-order BPS system (A.2) implies the second-order equations of motion for the different fields. Let us work in Einstein frame and write the total action as: where the action of type IIA supergravity is given by: .20) and the source contribution is the DBI+WZ action for the set of smeared D6-branes. Let us write this last action as in [23]. First of all, we introduce a charge distribution three-form Ω. Then, the DBI+WZ action is given by: where the DBI term has been written in terms of the so-called calibration form (denoted by K), whose pullback to the worldvolume is equal to the induced volume form for the supersymmetric embeddings. The expression of K has been written in [23]. Let us reproduce it here for completeness: where the e i 's are the one-forms of the basis corresponding to the forms (2.11) and (2.14) (see [23] for further details). Notice that the equation of motion for C 7 derived from (A.19) is just dF 2 = 2π Ω. Therefore, the Ω for our ansatz can be read from the right-hand side of (3.5). The Maxwell equations for the forms F 2 and The radial derivative of the dilaton can be put in terms of the derivative of Λ and h by using (A.1): where T sources µν is the stress-energy tensor for the flavor branes, which is defined as: In order to write the explicit expression for T sources µν derived from the definition (A.28), let us introduce the following operation for any two p-forms ω (p) and λ (p) : Then, by computing explicitly the derivative of the action (A.21) with respect to the metric, one can check that: It is now straightforward to compute explicitly the different components of this tensor. Written in flat components in the basis in which the calibration form has the form (A.22), we get: 5 , By using these values one can verify that, the Einstein equations (A.27) are indeed satisfied as a consequence of the first-order system (A.2). Notice that, for the massive flavored background of Section 6, dη/dr (and, therefore, T µν ) has a finite discontinuity at r = r q . It follows from (A.27) that the Ricci tensor R µν has also a finite jump at this point.

B Mass corrections in the UV
In this appendix we show how to obtain the first corrections to the conformal behavior of the metric and the dilaton from the UV asymptotic expansion of the master function W (x) written in (6.10). First we notice that the ratio η 0 /A 0 can be written in terms of b as: leading to a useful identity: Inserting the unit written in this way in the integral appearing in the expression of e g in (5.7), we get: where J is the following constant integral (depending onǫ): and G(x) is the function Using these results, we can immediately write: The function G(x) can be easily expanded for large x. At first non-trivial order we get: Using that, for large x, we easily arrive at the expansion of e g written in (6.16). Let us next find the large x expansion for f , which can be obtained by using the expansion of g and q in the relation e f = √ q e g .
We get: where the coefficient f 2 is given by: It is straightforward to demonstrate that (B.10) coincides with the value of f 2 written in (6.17).
Let us now analyze the UV expansion of h. First, we evaluate the integral appearing in (3.12): By combining this result with the expansions (6.16) and (6.12) we get: where L 0 is the UV AdS radius (3.39) and the coefficient h 2 is: One can readily check that the prefactor in (B.12) coincides with the one in (6.18) and that the coefficient h 2 written in (B.13) is the same as the one in (6.19).
Let us now obtain the UV expansion of the function Λ(x), defined in (A.1), which can be written as: 14) The UV expansion of this expression is: where the coefficient Λ 2 can be written in terms of g 2 as: By using the value of g 2 in (6.17), we can find a more explicit expression for Λ 2 : Let us next obtain the expansion of e φ by using e φ = e Λ h 1 4 . We get: where the prefactor is and can be shown to be the same as the asymptotic value of the dilaton written in (3.41).
The coefficient φ 2 is: which has been written more explicitly in (6.19).
Finally, let us write these results in terms of the r variable. At next-to-leading order, the relation between the coordinates x and r is obtained by integrating the differential equation: We get: This relation can be inverted as: By plugging the expansion (B.23) into (6.16) and (6.18) one easily arrives at (6.20) and (6.21).

C More on the entanglement entropy
Let us study analytically the entanglement entropy and the F function near the UV fixed point at R = 0. We shall represent S(R) in terms of a local functional L = L(H(x), G(x), ρ(x)) as: We will use this representation to compute the first-order variation of S when the background functions H and G and the embedding function ρ are varied around their UV values: is the function written in (7.23). At first order the corrections δH(x) and δG(x) can be parameterized as: The constants H 2 and G 2 in (C.3) can be related to the ones characterizing the behavior of the background at the UV: It is useful to define the following combination of H 2 and G 2 : The perturbation of the profile δρ(x) is the correction, at first-order, of the solution of (7.7) which satisfies δρ(x → ∞) = 0. By computing the first variation of (7.7) we find that δρ(x) is solution to the following second-order inhomogeneous differential equation: More explicitly, after using the UV values of H, G, and ρ, the differential equation (C.6) can be written as: Let us next calculate the first-order variation of the entanglement entropy. From (C.1) we get: where Π U V (x) is defined as: . (C.9) Using the explicit form of ρ U V we get: and it follows that: Therefore, the lower limit contribution to the last term in (C.8) vanishes. Let us now study the contribution to this term of the x → ∞ limit. We will now check that δρ(x) decreases as x → ∞ in such a way that: Thus, the upper limit of the last term in (C.8) also vanishes. To prove (C.12) we have to integrate the differential equation (C.7) and extract the large x behavior of the solution.
Amazingly, (C.7) can be integrated analytically. Its general solution can be written as the sum of two terms: where δρ P is a particular solution of the equation and δρ G is a general solution of the homogeneous part of (C.7). We have found a particular solution δρ P , which can be written in terms of hypergeometric functions and is given by: (C.14) We are only interested in the behavior of δρ P for large x. It is straightforward to prove that, for large x, δρ P can be approximated as: The general solution of the homogeneous differential equation which vanishes as x → ∞ is: where C is an arbitrary constant. For large x the function in (C.17) behaves as: Then, eq. (C.12) holds and the last term in (C.8) does not contribute to δS, as claimed above. Let us now calculate the other two contributions. First of all, the term due to the variation of H is given by: which, after using (7.24) to eliminate x * , becomes: which, after eliminating x * , gives: Thus, δS can be written as: where H 2 is the constant defined in (C.5). Using the definition (7.20), it follows that the change in the function F is: Taking into account (7.11) and (7.14) this expression can be written as: Therefore, we can write the F function near the UV fixed point as: where the constant coefficient c U V can be read from (C.26), Eq. (C.27) coincides with (7.31) when the former is written in terms of the dimension ∆ U V = 3 − b of the quark-antiquark bilinear operator in the massless flavored theory. Moreover, as F U V (S 3 ) is always positive, the sign of c U V depends on the sign of the coefficient H 2 . We have checked that H 2 is negative for all values of the deformation parameterǫ, in agreement with the expectation that F is maximized at the UV fixed point.

C.1 Entanglement entropy on the strip
Let us now evaluate the entanglement entropy in the case in which the region A is the strip − l 2 ≤ x 1 ≤ + l 2 of width l in the (x 1 , x 2 )-plane (see Fig. 10). In this case we consider a constant time surface Σ, whose embedding in the ten-dimensional space is determined by a function x = x(x 1 ). The induced metric on Σ is, where x ′ denotes the derivative of the holographic coordinate x with respect to the cartesian coordinate x 1 and the function G(x) has been defined in (7.4). The entropy functional for the strip of width l is given by: where V 6 is the volume of the internal manifold, whose value was given after (7.6), H(x) is the function defined in (7.5), and L 2 = dx 2 is the length of the strip. As the integrand in S(l) does not depend on the coordinate x 1 , the corresponding Euler-Lagrange equation admits the following first integral: x ′ L strip ∂x ′ − L strip = constant , (C.31) Figure 10: The surface Σ ends on a strip of length l at the boundary.
or, more explicitly, where H * = H(x = x * ) and x * is the holographic coordinate of the turning point. It follows from this last expression that x ′ is given by: Therefore, the width of the strip is given by the integral: 34) and the entropy S(l) is: The integral (C.35) for S(l) is divergent in the UV. Indeed, from the behavior of the functions H(x) and G(x) at large x (eqs. (7.10) and (7.11)), it follows that, when x → ∞, the integrand in (C.35) behaves as: and S(l) therefore diverges as x 1 b max , where x max is the upper limit of the integral, which can be regarded as a UV cutoff. More explicitly, one can rearrange the integral in (C.35) as: where the divergent term is explicitly shown. We now define the rescaled functionsĜ(x) andĤ(x) as:Ĝ Then, the finite term in S(l) can be written as: whereas the divergent term takes the form: We have evaluated numerically the right-hand side of (C.39) as a function of the strip width l. The result is displayed in Fig. 11. One notices from these results that S f inite (l)/L 2 is negative and diverges as −c/l for small l, where c is a constant. Actually, for small l the surface Σ is in the UV region of the geometry and one can evaluate the entropy analytically in this limit. This is the subject of the next subsection, where we show that the constant c is proportional to the free energy on the three-sphere of the massless flavored theory (see eq. (C.51) below). The analytical UV calculation is compared to the numerical results in Fig.  11.

C.1.1 UV limit
In order to study the UV limit it is convenient to change variables in the integrals (C.34) and (C.39). Let us introduce a new variable z, related to x as x = x * z. Then, l can be represented as: while S f inite (l) is given by: Let us now obtain the expressions of l and S f inite (l) in the limit in which x * → ∞. In this case the argument of the functionsĜ andĤ in the integrals is always large and one can takê We get: where I 1 is the following integral: Using (C.44) in (C.43), we get: Similarly, the finite part of the entropy is: where I 2 is the integral: (C.47) Plugging (C.47) into (C.46), we arrive at: Eliminating x * in favor of l, we get: However, from the definition of G ∞ and H ∞ in (7.11) it follows that: where F U V (S 3 ) is the free energy of the massless flavored theory on the three sphere (see (7.17)). It follows that: which is the result we were looking for. From the comparison of Fig. 11 to the numerical results we conclude that (C.51) is a good description of the entropy in the small l region. Notice that the coefficient of 1/l on the right-hand side of (C.51) is determined by the free energy of the massless flavored theory. This is expected on general grounds in the UV region, since the masses of the quarks can be neglected in the high-energy regime.

C.1.2 IR limit
Let us now evaluate S f inite (l) in the regime in which l is large. In this case the surface Σ penetrates deeply in the bulk and its tip is near the origin (i.e., x * is small). We will proceed as in Section 7. with x a ≫ x * , H * ∝ x 4 * ≪ H(x) ∝ x 4 and one can neglect the terms containing H * in the integrals. Then, the strip width l can be approximated as: When, x a ≫ x * and x * → 0 we can extend the first integral in (C.52) up to ∞ and neglect the contribution of the second integral (which is proportional to √ H * ∝ x 2 * ). Then, l can be approximately written as: and, as expected, a small value of x * leads to a large value of l. Similarly, we can perform the same type of manipulations in the expression (C.39) of S f inite . We find: Performing the first integral in (C.54) in the limit x * → 0 for fixed x a , we arrive at: where S ∞ is the constant: (C.56) We can now eliminate in (C.55) the turning point x * in favor of l by using (C.53). We find: where we have written the result in terms of F IR (S 3 ). The first term in (C.57) is just the strip entanglement entropy for the unflavored theory. The constant S ∞ represents the asymptotic value of S f inite (l)/L 2 as l → ∞. One can approximate this constant by taking x a → 0 in (C.56). After using (7.14) to relate G ∞ and H ∞ to F U V (S 3 ), we arrive at: Using this result, we get the following dependence of E qq with the distance d where we have represented this relation in terms of the rescaled quantities E qq /r q and r q d/L 2 0 . Notice that the leading term (the first term on the right-hand side of (D.14)) is given by the potential of the massless flavored background, as expected.

D.2 IR potential
We now estimate the qq potential for large separations. In this case we will content ourselves to compute the leading order contribution. For large qq separations the hanging fundamental string penetrates deeply in the geometry, which is equivalent to saying that x * is small. We will follow an approach similar to the one in Sections 7.2 and C.1.2 and we will split the [x * , ∞] interval of the integrals (8.6) and (8.10) as [x * , ∞] = [x * , x a ] ∪ [x a , ∞] with x a < 1 being small. We will assume that x a is small enough so that the background functions are well approximated by (7.32) in the interval [x * , x a ] . Then, we can estimate the integral (8.6) for d as: where L ABJM is the unflavored AdS radius (2.3) and we have used the fact that for x * ≪ 1 we have that h * ≫ h(x) when x ∈ [x a , ∞] and, therefore, in this interval we can neglect h(x) in the square root of the denominator of the integrand in (8.6). The first integral in (D.15) can be done analytically, yielding the result: For small x * , at leading order, we get from (D.16): 17) which confirms that r q d is large when x * is small. We can make similar approximations in the integral (8.10) with the result: Performing explicitly the integral, we get: 19) which, at leading order, becomes: Eliminating x * by using (D.17) we arrive at the estimate (8.14).

E Asymptotics of the two-point functions
In this appendix we study the renormalized geodesic distance, and the corresponding twopoint function of bulk operators with large conformal dimensions, in the UV limit of small separation l and in the IR regime in which r q l → ∞. In the former case, our study will serve to fix the normalization constant C of (9.10), as well as the analytic form of the correlator near the UV fixed point. This is the case we analyze first in the next subsection.

E.1 UV behavior
We will now obtain the form of the correlator in the UV limit in which the turning point x * is large and the geodesic does not penetrate much into the bulk of the geometry. In order to study this limiting case, it is convenient to perform a change of variables in the integral (9.7) and write it as: (E.1) Similarly, the renormalized geodesic length can be represented as: where z max = x max /x * . When x * is large the argument of the functions in the integrals (E.1) and (E.2) is large and one can use the UV asymptotic expressions (6.18) and (7.10). Therefore, we can approximate l as: where the function f ∆ parameterizes the deviation from the conformal behavior near the UV fixed point and is given, at first order, by where the new constant c ∆ is just:

E.2 IR behavior
Let us now analyze the two-point functions in the IR limit, which corresponds to taking x * → 0. In this case the geodesic penetrates deeply in the bulk. We will proceed similarly to what we did in Section D.2 for the qq potential in the IR. Accordingly, we split the interval [x * , x max ] as [x * , x max ] = [x * , x a ] ∪ [x a , x max ] for some x a < 1 small. Near the turning point, i.e., when x * ≤ x ≤ x a , we can approximate the functions h(x) and g(x) as: where L IR = L ABJM is the AdS radius of the unflavored ABJM solution and γ and c are the parameters of the unflavored running solution. It follows that, in this IR region, the function G(x) behaves as: Moreover, the dilaton φ IR is constant and given by the ABJM value φ ABJM written in (2.4). Suppose that we have chosen some x a ≫ x * . Let us then split the integral (9.7) for l as: where x max should be sent to +∞ at the end of the calculation. Let us approximate the first integral in (E.24) by taking the functions in the deep IR, as in (E.22). We obtain: (E.25) For small x * /x a we can approximate this integral as: To evaluate l U V approximately we notice that h * ∝ x −4 * and, therefore, it is large for x * → 0. Then, we can neglect the one inside the square root and approximate l U V as: h(x) 1 4 G(x) . (E. 27) Since the integral in (E.27) converges and is independent of x * , it follows that l U V ∼ x * and, therefore, it can be neglected with respect to the large value of l IR ∼ 1/x * . Thus, we take Notice that it follows from this equation that r q l is large if x * is small, as it should be in the IR regime. Let us next perform a similar analysis for the renormalized geodesic length L r . First, we consider the IR integral: To evaluate the UV integral we proceed similarly to what we did for the integral for l and, in this case, we only keep the one inside the square root. Then: which is just (9.13). Let us next define the following quantities: Then, after using the relation (E.28) to eliminate x * in favor of l, we get: Then, it follows that the IR limit r q l → ∞ of the two-point correlator is as in (9.14), where N is the normalization constant given by: (E. 35) It turns out that, in the expression of the integrals in (E.33) we can take the limits x a → 0 (in I IR ) and x max → ∞ (in I U V ). Actually, it can be easily proved that, after taking these limits, I IR and I U V can be recast as: Notice that the form of the correlator is consistent with the fact that the conformal symmetry is recovered in the IR limit r q l → ∞. The non-canonical normalization factor N is due to the fact that we chose to renormalize L in such a way that the correlator is canonically normalized in the opposite UV limit r q l → 0.

F WKB mass levels
Consider the following differential equation for the function R(x): where x takes values in the range x * ≤ x ≤ ∞, m is the mass parameter and P (x) and Q(x) are two arbitrary functions that are independent of m. We will assume that near x ≈ x * , ∞ these functions behave as: where P i , Q i , s i , and r i are constants. By a of suitable change of variables, the differential equation (F.1) can be converted into a Schrödinger equation, which only admits a discrete set of normalizable solutions for a set of values of m parameterized by a quantum number n ≥ 0. The mass levels for large values of n can be evaluated in the WKB approximation [44], with the result: where ξ = ξ(x * ) is the integral: In our case, the functions P (x) and Q(x) which correspond to the fluctuation equation (10.13) are: and it is immediate to extract the exponents s 1 , s 2 , r 1 , and r 2 : Using these values one immediately finds that the WKB mass levels are given by (10.14).

F.1 Asymptotic spectra
Let us now evaluate analytically the meson spectrum in the two limiting cases in which the sea quark mass m q = r q /2π is small or large compared to the valence quark mass µ q . Notice that e g ∝ m q (see eqs. (5.6) and (5.7)). Therefore, (10.4) can be regarded as giving the relation µ q /m q as a function of x * . Let us consider first the case in which m q is small. Clearly, when m q → 0 for fixed µ q one necessarily must have x * large. Actually, integrating (10.4) by using the asymptotic expression of g written in (6.16) at leading order, we get: where the constant κ has been defined in (6.15). Thus, when m q → 0 for fixed µ q the x coordinate of the tip of the flavor brane increases as: Then, we have: µ q ≈γ + 1 2 m q x * , (F. 16) and it follows that, for fixed µ q and large m q , the coordinate of the tip of the flavor brane behaves as: Thus, for large sea quark mass x * → 0 and thus the dynamics of the fluctuating flavor brane is dominated by the IR, where the solution corresponds to the running solution of the unflavored system. In this case, we have at leading order near x ∼ 0: where, in the last step, we used that c/γ ≈ x * /r * . Thus, and, after evaluating the integral (F.9) and writing the result in terms of the 't Hooft coupling λ and the valence quark mass µ q , we have: which is exactly the mass spectrum of vector mesons in the unflavored ABJM model [13]. By combining (F.14) and (F.20) we obtain the UV/IR mass relation written in (10.16). It is now straightforward to obtain this ratio for small and large values of the deformation parameter. Indeed, for smallǫ one can expand m