Gravity dual of a multilayer system

We construct a gravity dual to a system with multiple (2+1)-dimensional layers in a (3+1)-dimensional ambient theory. Following a top-down approach, we generate a geometry corresponding to the intersection of D3- and D5-branes along 2+1 dimensions. The D5-branes create a codimension one defect in the worldvolume of the D3-branes and are homogeneously distributed along the directions orthogonal to the defect. We solve the fully backreacted ten-dimensional supergravity equations of motion with smeared D5-brane sources. The solution is supersymmetric, has an intrinsic mass scale, and exhibits anisotropy at short distances in the gauge theory directions. We illustrate the running behavior in several observables, such as Wilson loops, entanglement entropy, and within thermodynamics of probe branes.


Introduction
The holographic AdS/CFT correspondence relates strongly interacting quantum field theories (QFT) to classical gravity theories in higher dimensions [1]. Apart from being a conceptual breakthrough, this duality has become an important and versatile tool in the study of the possible states of matter (for reviews see [2]).
There are basically two types of approaches to implement the holographic idea. In the socalled bottom-up approach, in order to model a d-dimensional QFT, one considers a gravity system in d + 1 dimensions. This is somehow the minimal version of the correspondence and has been very successful in describing many interesting phenomena with models that evade the rigid constraints of string theory. On the contrary, the top-down models use the full machinery of string theory. They employ ten-dimensional gravity solutions of the corresponding equations of motion of type II supergravity. These ten-dimensional solutions are typically more difficult to obtain. However, the extra internal directions encode precious information which is lost in many phenomenological bottom-up approaches. Moreover, knowing the brane configuration which gives rise to the supergravity background allows to determine its field theory dual. Actually, one can engineer such holographic duals from the corresponding brane setup.  In this paper we will employ string theory techniques to find gravity duals of systems composed by multiple parallel (2 + 1)-dimensional layers in a (3 + 1)-dimensional ambient theory (see Fig. 1). The ambient 4d theory will be N = 4 super Yang-Mills theory with gauge group SU (N c ), which is realized in a stack of coincident N c color branes. The multiple layers are obtained by adding parallel D5-branes sharing two spatial directions with the stack of color branes, according to the array: The field theory corresponding to the array (1.1) is well-known [3][4][5][6]. It consists of a supersymmetric theory with matter hypermultiplets (flavors) living on the (2 + 1)-dimensional defect and coupled to the ambient N = 4 theory. To construct a multilayer structure of the type shown in Fig. 1 we consider a continuous distribution of D5-branes along the third direction in (1.1). In the absence of D5-branes, the geometry generated by the D3-branes is AdS 5 × S 5 . We want to obtain a supergravity solution which includes the backreaction of the flavor D5-branes. This backreaction can be obtained by using the techniques developed to study unquenched holographic flavor (see [7] for a review and references). In some cases the solutions found with these techniques are analytic and preserve some amount of supersymmetry. In this approach the D5-branes of the array (1.1) are smeared both in the cartesian direction orthogonal to the defect and on the internal directions. The dual gravity background can be obtained by solving the supergravity equations of motion with D5-brane sources. These backgrounds are dual to anisotropic systems, since there is a distinct field theory direction, the direction orthogonal to the defect. Let us also comment on the novelty of our approach. There are several previous holographic works which provide anisotropy, see, e.g., [8][9][10][11][12][13][14][15][16][17][18][19]. The main difference between our model and other anisotropic models in the literature is that in our case the anisotropy is produced by the presence of dynamical objects (the D5-branes of the multiple layers) and not by fluxes or fields depending anisotropically on the coordinates. The D3-and D5-branes of the array (1.1) can be separated in the directions 789 transverse to both types of branes. When this separation is zero the mass of the hypermultiplets living on the defect vanishes, i.e., we have massless flavors. This D3-D5 massless flavor case was considered in [20], where an analytic supersymmetric solution was found that displays a Lifshitz-like anisotropic scaling invariance. The non-zero temperature generalization of the scaling solution was found and studied in [21]. In this paper, we study the massive flavor case of the D3-D5 intersection (1.1). The solutions we present here preserve the same supersymmetry as the massless flavor case, but do not possess the scaling invariance of the latter.
The gravity duals of theories with massive flavors are running solutions which naturally represent a renormalization group flow. This flow is generated by changing the quark mass (see [22] for an example of these massive flavored backgrounds in the ABJM theory). When the mass of the quarks is very large, the flavors decouple and we expect to recover the unflavored solution (AdS 5 × S 5 in our case). On the other hand, when the flavors are massless we obtain the anisotropic scaling solution of [20]. For a finite non-vanishing value of the quark mass we expect to get a background interpolating between these two solutions: Unflavored in the IR and massless flavored in the UV. Once we have the background at our disposal we can study the effects of the flow on different observables. In general, we expect to obtain the results corresponding to the isotropic AdS 5 × S 5 solution in the IR and to be able to tune the amount of anisotropy by changing some of the parameters of our solutions. In our analysis of several observables we will find that, indeed, the UV behavior is determined by the scaling solution of [20], whereas the long distance IR behavior depends on a free parameter of our geometry.
The rest of this paper is organized as follows. In Sec. 2 we present our ansatz for the metric, for the dilaton, and for the forms. All functions of the model depend on a master function, which in turn satisfies a second order differential equation. This equation follows from supersymmetry analysis detailed in App. A. A crucial ingredient entering the equations is the so-called profile function, which encodes the distribution of sources along the holographic coordinate. To determine this function for D5-brane sources one needs to analyze in detail the embeddings of the D5-branes and their kappa symmetry. This analysis is deferred to App. B.
In Sec. 3 we tackle the problem of integrating the master equation. This task requires the redefinition of some of the functions and a change of variables. In its final form our solution depends on a constant parameter which characterizes the IR deformation of the metric. In Sec. 4 we begin our study of the observables in our background. In this section we study the Wilson loops and the potentials for quark-antiquark pairs, when these particles are in the same layer or separated in the direction orthogonal to the layers. In Sec. 5 we do a similar analysis for the entanglement entropy of slabs.
In Sec. 6 we explore our supergravity solution with a probe D5-brane with a worldvolume gauge field dual to a chemical potential. We analyze the zero temperature thermodynamics of the probe and, in particular, the UV-IR flow of the speed of sound. This is not the only interesting configuration of probe branes with non-zero worldvolume gauge field. An interested reader is invited to App. D where we consider D5-branes with worldvolume flux along the internal directions. This internal flux induces a bending of the probe brane along the direction x 3 in the array (1.1), which can be interpreted as a recombination of the flavor D5-branes and the color D3-branes, realizing the Higgs branch of the theory [23,24]. Finally, in Sec. 7 we summarize our results and discuss some research lines for the future.

Supergravity ansatz
In this section we review the supergravity setup of [20], corresponding to the array (1.1) of D3-and D5-branes. More details are given in App. A. The D3-branes are color branes which generate an AdS 5 × S 5 space, whereas the flavor D5-branes create a codimension one defect in the (3 + 1)-dimensional gauge theory and, when the backreaction is included, the original AdS 5 × S 5 metric gets deformed. The D5-branes are homogeneously distributed in the internal space in such a way that some amount of supersymmetry is preserved. When the S 5 space is represented as a U (1) bundle over CP 2 , the deformation of the five-sphere depends on a single radial function, which measures the relative squashing between the fiber and the base of the deformed S 5 . Choosing a convenient radial coordinate ζ (with boundary corresponding to ζ = ∞), the ten-dimensional Einstein frame metric takes the form: where φ is the dilaton of type IIB supergravity, h is the warp factor, and f is the squashing function. These functions are assumed to depend only on ζ. Moreover, A is a one-form on CP 2 which implements the non-trivial U (1) bundle. The Minkowski directions x 1 and x 2 are parallel to the defect, whereas x 3 is orthogonal to it. Besides the metric and the dilaton, the type IIB supergravity solution contains a RR fiveform F 5 and a RR three-form F 3 . The former is self-dual and given by the standard ansatz in terms of the dilaton φ and warp factor h: Clearly, d F 5 = 0, since the D3-branes have been replaced by a flux in the supergravity solution. On the contrary, the D5-branes are dynamical and are governed by the standard DBI+WZ action which, in particular, contains the term: where 1/T 5 = (2π) 5 g s α 3 and C 6 is the six-form potential for F 7 = −e φ * F 3 (the hat over C 6 denotes its pullback to the D5-brane worldvolume M 6 ). Therefore, S W Z contributes to the equation of motion of C 6 or, equivalently, to the Bianchi identity of F 3 . Indeed, let us write the six-dimensional integral in (2.3) as a ten-dimensional integral: where Ξ is a four-form with support on the worldvolume of the D5's and with legs along the directions orthogonal to M 6 , which is just the RR D5-brane charge distribution. The equation of motion for C 6 is: where 2 κ 2 10 = (2π) 7 g 2 s α 4 . In the smearing approach Ξ does not contain Dirac δ-function singularities. Its form can be obtained once the ansatz of F 3 compatible with supersymmetry is fixed. This has been done in [20], a result which we now review.
(2.6) Then, F 3 is given by: where Q f is a constant and p(ζ) is an arbitrary function of the holographic coordinate ζ. By computing the exterior derivative of F 3 we get its modified Bianchi identity: Comparing (2.8) and (2.5) we can extract the D5-brane charge distribution Ξ which, in what follows, we will refer to as the smearing form. Clearly, Ξ does not depend on the x 3 coordinate, although it contains dx 3 in its expression. This means that we are continuously distributing our D5-branes along x 3 , giving rise to a system of multiple (2 + 1)-dimensional parallel layers. Moreover, the function p(ζ) introduces a profile of the charge distribution in the holographic coordinate. Notice that the D3-and D5-branes in the array (1.1) can be separated in the 789 directions. In principle, we could have an arbitrary distribution of D5-branes in these transverse coordinates, which is reflected in the fact that the profile function p(ζ) is arbitrary. However, for a stack of flavor D5-branes with the same quark mass the function p(ζ) has a well-defined form (see App. B) and Q f is related to the density of smeared branes along the direction x 3 . As shown in App. B, if we distribute N f D5-branes along a distance L 3 in the third cartesian direction, then (see (B.59)). The preservation of two supercharges for our ansatz imposes a system of first-order differential equations in the radial variable. These BPS equations are reviewed in App. A, where it is shown that they can be reduced to a single second-order differential equation for a master function W (ζ). This equation is: From W we can reconstruct the full solution. The squashing function f (ζ) is given by: while the dilaton is: A nice way of measuring the deformation of the metric (2.1) with respect to the AdS 5 × S 5 geometry is obtained by considering a squashing factor, which we define as follows It is clear from our ansatz (2.1) that q represents the relative size of the U (1) fiber with respect to the CP 2 base. In terms of the master function W , q is given by: (2.13) Once f and φ are known, the warp factor h can be obtained by integrating the following first-order differential equation: where Q c is related to the number N c of D3-branes as Q c = 16 πg s α 2 N c . In what follows we will study several solutions of the master equation (2.9).

Unflavored solution
Let us consider the master equation in the case in which there are no flavor brane sources. It turns out that the general solution of (2.9) can be analytically found in this case. Indeed, when Q f = 0 the master equation (2.9) can be trivially integrated once as: A further integration gives: where C and b are constants. Plugging this expression of W into the right-hand-side of (2.11) one readily verifies that the dilaton is constant and given by: Moreover, the function f and the squashing function q are given by: 18) and the warp factor h can be obtained by integrating (2.14). This solution coincides with the general unflavored one found in [20]. For b = 0 this geometry is just AdS 5 × S 5 . When b = 0 the solution approaches AdS 5 × S 5 in the UV. If we take b to be real and positive, then the minimal value of ζ is ζ = b and the metric has a blown-up CP 2 cycle at ζ = b. It was argued in [25] that this b = 0 background is dual to the superconformal N = 4 field theory deformed by the VEV of a dimension 6 operator.

Massless flavored solution
Let us now consider the massless flavor case with Q f = 0 and p = 1. In this case it is possible to find a special solution of (2.9): Indeed, by plugging this ansatz into the master equation we readily get that the exponent α is given by: Similarly, we can obtain the value of the constant A. The final formula for W is: Using this expression in (2.10) and (2.11) we arrive at the following values of f and φ: It is straightforward to verify that this solution coincides with the one found in [20] for massless flavors. 1 The warp factor for this solution is: whereR is the same as in [20], Since the profile function p is constant, this solution represents massless smeared flavors extending all the way down to ζ = 0. As shown in [20], the background corresponding to the master function (2.21) is invariant under a set of anisotropic scale transformations in which the x 3 coordinate transforms with an anomalous scaling dimension. Moreover, the squashing factor (2.12) q is constant and equal to 3 2 √ 2 ≈ 1.06, see (2.22). The purpose of this paper is to find solutions corresponding to massive flavors, which should interpolate between the unflavored solution of subsection 2.1 at the IR and the scaling solution studied in this subsection at the UV. We start to discuss these solutions in the next subsection. 1 The relation between ζ and the radial variable r used in [20] and in the ansatz (A.1) is ζ = 3 2 √ 2 r.

D3
⇣ q miércoles, 7 de noviembre de 18 ⇣ q miércoles, 7 de noviembre de 18 Figure 2: On the left we depict a localized embedding of a D5-brane with a separation ζ q from the stack of color D3-branes. On the right, several D5-branes with different orientations and the same distance ζ q generate a cavity ζ ≤ ζ q which does not contain flavor sources.

Massive flavors
In the holographic approach the mass of the fundamentals is related to the distance between the color and flavor branes (in our case D3's and D5's, respectively). When these two sets of branes are separated, the fundamentals are massive and there is an IR region of the geometry which is not occupied by the the flavor brane and, as a consequence, the D5-brane charge density vanishes there. In the smearing setup we have many D5-branes with different orientations in the internal space which, if they correspond to flavors with the same mass, should have the same separation from the D3-branes. As illustrated in Fig. 2, the sourceless region has a ζ coordinate less or equal to some value ζ q , a region we will call cavity in the following. The profile function p(ζ) vanishes inside the cavity and should approach the value appropriate for massless flavors, i.e., p → 1 as ζ → ∞. To determine the explicit form of the function p(ζ) one has to specify the set of source D5-branes of our smeared distribution. This is done in detail in App. B. The final result for the profile function is: Notice that p(ζ) is continuous at ζ = ζ q and asymptotically p(ζ → ∞) = 1. The master function W for the profile (2.25) must be obtained numerically. However, we can expand (2.9) in a series expansion about any radial coordinate ζ. Let us start in the neighborhood of ζ = ∞. Indeed, the profile function (2.25) can be expanded in powers of ζq ζ as: Plugging this expansion into the master equation (2.9) and integrating order by order, we get the following solution for the master function W (ζ): The asymptotic UV expansions of the different functions of the background are easily obtained from (2.27). We have collected these expansions in App. A. Another useful expansion is the neighborhood where the sources set in, i.e., at the edge of the cavity. We have hence perturbatively solved the master equation for ζ ≥ ζ q and ζ − ζ q small; details are relegated in App. A.
Recall that inside the cavity, i.e., for ζ ≤ ζ q , we have the analytic solution for the master function. This solution was written in (2.16) and depends on two parameters C and b. We want to extend this solution for ζ ≥ ζ q and determine the values of C and b for which the function W behaves as in (2.27) for ζ → ∞. This matching has to be done numerically, for which we implemented the shooting technique accompanied with a convenient change of variables. We will discuss this in the next section.

Integration of the master equation
Let us introduce a new holographic variable x, related to ζ as follows: Clearly, for our interpolating solutions x ≥ 0. In this new variable the cavity (i.e., the region without flavor branes) corresponds to 0 ≤ x ≤ x q , where x q is the edge of the cavity given by: Notice that x q depends on the ratio between the deformation parameter b in the sourceless region and the location of the boundary of this unflavored region in the ζ variable. Moreover, as b ≤ ζ q , we have that: Let us write the background in terms of the new variables. It is convenient to absorb the constant C appearing in the master function (2.16) inside the cavity. Accordingly, we definê W as: In the region without flavor sources this function is given by:  Z(x q ) Figure 3: Plot of the function Z(x q ) = C 3 2 ζ q /Q f , obtained by the shooting method of the master equation (3.8). This function has as asymptotics Z(x q = 0) ≈ 1.38 and Z(x q = 1) ≈ 1.26. Notice that Z(x q ) is almost constant in the interval 0.5 ≤ x q ≤ 1.
This function vanishes at the origin x = 0 andŴ and its derivative take the following values at the edge of the cavity: Notice that, for x q = 1, i.e., in the case in which the cavity is largest,Ŵ (x) = 1 inside the cavity and the solution becomes (isotropic) AdS in this region. This implies that x q measures the amount of anisotropy of our solution.
The profile function in terms of the x variable is: Moreover, the master equation forŴ becomes: This equation is solved analytically by (3.5) inside the cavity, where 0 ≤ x ≤ x q . For x ≥ x q we can solve (3.8) numerically by imposing the initial conditions (3.6) and by requiring that as x → ∞ it asymptotically behaves as the massless flavored solution, Notice that (3.8) depends on two parameters x q and C 3 2 ζ q /Q f . In order to have a solution with the asymptotic behavior (3.9) these two parameters must satisfy a relation, which can be determined numerically. This relation is where the function Z(x q ) is determinend numerically for 0 ≤ x q ≤ 1. Our results for this function have been plotted in Fig. 3. We observe that Z(x q ) is a decreasing function of x q and, for given values of Q f and the constant C, ζ q reaches its maximum at x q = 0. Notice also that x q gives the size of the cavity in the x variable. This size should be related to the quark mass m q . In our holographic setup m q can be determined by evaluating the Nambu-Goto action of a fundamental string extended along the holographic direction, from the origin of the space to the tip of the flavor brane. In the ζ variable we have: where the e φ 2 factor is due to the fact that we are working in the Einstein frame. By using the metric ansatz (2.1), as well as (2.10) and (2.11), we obtain m q in terms of an integral involving the master function W : (3.12) Putting α = 1 from now on and writing the result in terms of x q , we get: In Fig. 4 we plot C 2 m q /Q f as a function of x q . We observe in this plot that m q grows with x q and that m q = 0 for x q = 0, i.e., when the cavity has zero size. Moreover, we can increase the quark mass by decreasing C (C → 0 corresponds to m q → ∞ for x q = 0). The precise relation between m q and the constant C depends on the value of x q . For x q ∼ 0, 1 one can expand the right-hand side of (3.13) and find: (3.14) All the functions of the background can be written in terms of Z andŴ . For example, the squashing factor q is: whereas f , g, and the dilaton φ are: where, in the second step, we have written these functions in terms of q(x).
In Fig. 5 we plot the functionŴ for different values of x q . We notice that, as x q → 1, W is almost constant and equal to 1 in the interior of the cavity, where it is given by (3.5). Moreover, outside the cavity, i.e., for x ≥ x q it rapidly approaches the asymptotic expression (2.27), which can be rewritten in terms of x and x q as: The UV expansion of f , φ, h, and q in the x variable can be obtained from The comparison of the numerical values for x ≥ x q and those given by (3.17) shows that the agreement with the UV expansion is better if x q is small, since in this case the background is closer to the massless scaling solutions (for x q = 0 there is no cavity). On the other hand, for x q close to one the flavor effects are smaller in the IR, sinceŴ is almost constant and equal to one inside the cavity. Indeed, from (3.2) we get that x q ∼ 1 implies that the IR deformation parameter b is small. This is consistent with the fact that, for a given value of C, the quark mass is maximal in this case (see figure 4).
Let us now write the metric in the x variable. First of all, we define the rescaled warp factor and cartesian coordinates as: Then, in terms of these hatted variables, we have: We have only one free parameter, x q , corresponding to a family of geometries. Actually, in our unquenched flavored background one would expect the geometry to depend on two quantities, the amount of flavors and their mass. We have rescaled out the quark mass with the definition of the hatted quantities in (3.18) and therefore x q will serve us to parametrize the amount of flavors. It is interesting to write down the relation between the hatted and unhatted cartesian coordinates in terms of the quark mass m q . Taking into account that C ∼ m − 1 2 q , we get for the longitudinal (x 1 x 2 ) and transverse (x 3 ) directions: (3.20) Therefore, for fixed distances x and x ⊥ , taking m q → 0 is equivalent to considering the UV x ,x ⊥ → 0 region in the hatted variables, whereas taking large m q amounts to zooming in the IR region of large x and x ⊥ . For fixed quark mass, x q is the parameter that controls the flavor effects in the IR. To illustrate how the flavor branes deform the metric as we move in the holographic direction, we have plotted in Fig. 6 the squashing function q(x) for different values of x q . For x q ∼ 1 the function q(x) is nearly constant and equal to one in the sourceless region x ≤ x q and grows monotonically for x ≥ x q , until it reaches its asymptotic UV value q = 3/2 √ 2. It follows from these results and those represented in Fig. 5 that the geometry becomes more and more isotropic inside the cavity as we increase the parameter x q .
In the following sections we will study our gravity dual by computing several observables with the purpose of exploring their change as we move from the UV to the IR as we vary the size of the cavity x q .

Wilson loops
In holography, the potential energy between a "quark" and an "antiquark" is obtained from the solution of the equation of motion of a fundamental string hanging from the UV boundary q(x) Figure 6: The squashing function q(x) for x q = 0.005 (black), x q = 0.5 (blue) and x q = 0.9 (brown). As x → 0 the squashing factor goes to zero for all the values of x q , except for the case x q = 1. In that particular case q(x) remains finite and equal to one inside the cavity. [26,27]. These equations are derived from the Nambu-Goto action: where g 2 is the induced metric on the worldsheet of the string. We will calculate the qq potential in two cases. First we will consider the intralayer case, in which the quark and the antiquark are in the same layer and have the same value of the coordinate x 3 . After this we will consider the interlayer configuration, where the quarks are separated in the anisotropic direction.

Intralayer potential
We take (t,x 1 ) as worldvolume coordinates of a fundamental string and we will consider an ansatz with x = x(x 1 ) with the other cartesian coordinates being constant. The induced metric on the two-dimensional worldsheet is: where the prime denotes derivative with respect tox 1 and we have denoted Z q ≡ Z(x q ). The Nambu-Goto action of the string is: whereT = dx 0 and the Lagrangian density L is: As L does not depend explicitly on the holographic coordinate x, one has the following first integral: from which we get: where q 0 ,ĥ 0 , andŴ 0 are the values of the functions q,ĥ, andŴ , respectively, at the turning point x = x 0 . From this relation we obtain x as: which can be straightforwardly integrated to give: The (hatted) quark-antiquark distance at the boundary is: Let us now calculate the on-shell action of the fundamental string. Plugging the solution into the action, we get: . (4.10) As usual, this on-shell action is divergent and must be regulated. We do it by subtracting the action of two straight fundamental strings stretched from the origin x = 0 to x = x max : The quark-antiquark potential is then given by: where we have used the relation between T andT : More explicitly: q d for x q = 0.1 (red), x q = 0.5 (blue), and x q = 0.9 (brown). The dashed curve is the UV limit (4.18). Notice that the maximal separation forqq increases with x q .
We have numerically evaluated the potential V qq as a function ofd . The results are presented in Fig. 7 for different values of x q . We notice that all curves become coincident for smalld . Asd ∼ m 3 4 q d (see eq. (3.20)), one expects to recover the massless scaling solution in the UV domaind → 0. Indeed, we prove below that V qq ∼d − 4 3 in the UV region of small d . This behavior matches the one obtained numerically, as shown in Fig. 7. As we move towards the IR by decreasing the turning point coordinate x 0 and increasingd , we obtain that there is a maximal value ofd (corresponding to a minimal value x min 0 < x q of the turning point coordinate). For x 0 < x min 0 the dominant configuration is a disconnected one, in which the two ends of the string go straight from the boundary to the origin. This behavior has been obtained previously in backgrounds dual to unquenched flavors [28][29][30][31] (in other types of backgrounds, see also [32][33][34][35][36]). Indeed, dynamical quarks produce string breaking and a maximal length due to pair creation. In our case this breaking is not produced in the scaling solution with m q = 0. Moreover, the critical distance at which the string breaks grows with x q , as is also evident in Fig. 7, and becomes very large when x q ∼ 1. This is easy to understand since the breaking occurs when the string penetrates deeply in the cavity, whose size is maximal when x q ∼ 1, and the integrals (4.9) and (4.14) get their main contribution from the sourceless region inside the cavity. For large enough values ofd the dominant configuration is the disconnected one with zero energy, which means that the external quarks are completely screened by dynamical quarks popping out from the vacuum. A recent interesting work [19], in a seemingly unrelated context of holographic QCD, parallels our findings. The authors of [19] demonstrated that large amounts of anisotropy will completely screen the interactions between quarks and anti-quarks, while in the absence of anisotropy the model would otherwise be confining.

UV limit
Let us now evaluate the potential in the UV limit in which d is small (large x 0 ) and our embedding is close to the boundary. In this limit we can use that, at leading order in the UV, the squashing function q is constant and thatŴ ∼ 20), and (A.21)). We will use these values to calculate the integrals ford and V qq in (4.9) and (4.14). At leading order we get the following relation betweend and x 0 : Similarly, we approximate the potential V qq by the following integral: Performing the integral, we get: By using the relation (4.15) we can eliminate x 0 in favor of the qq distanced . After some calculation we get: which coincides with the result found in [21] for the massless scaling background. In Fig. 7 we show that the potential (4.18) does indeed coincide with the numerical results in the UV domaind → 0 for all values of the parameter x q .

Interlayer potential
We now consider a Wilson loop that extends in the x 3 direction with x 1 and x 2 constant, which corresponds to two fundamentals located at different layers. Accordingly, we take (t,x 3 ) as worldvolume coordinates and consider an ansatz in which x = x(x 3 ). The twodimensional induced metric is now: where the prime now denotes derivative with respect tox 3 . The Nambu-Goto action is: The first integral derived from L is: where again q 0 ,ĥ 0 , andŴ 0 are the values of the functions q,ĥ, andŴ , respectively, at the turning point x = x 0 . Then: Integrating this equation we get: (4.24) Therefore, the quark-antiquark distance along the direction transverse to the layers is: (4.25) The unregulated on-shell action for this case is: Proceeding as in the intralayer case to regulate this action, we arrive at the following quarkantiquark potential: q d ⊥ for x q = 0.1 (red), x q = 0.5 (blue), and x q = 0.9 (brown). The dashed curve is the UV potential (4.31). On the right we compare the intralayer and interlayer potentials for x q = 0.1 (red), x q = 0.5 (blue), and x q = 0.9 (brown). The continuous (dashed) curves correspond to the intralayer (interlayer) potentials.
The numerical results for the interlayer potential have been plotted in Fig. 8. They are qualitatively similar to the intralayer case of Fig. 7. In the UV regiond ⊥ ∝ m ⊥ , in agreement with our analytic calculation of Sec. 4.2.1. In this case there is also a maximal length which increases with x q . We have also compared the intralayer and interlayer potentials for the same value of x q . These potentials are plotted together on the right panel of Fig. 8, where we notice that they have very different behavior in the UV but they become very similar in the IR (cf. also [19]). This IR similarity increases as x q approaches its maximal value x q = 1, which is consistent with the fact that for large values of x q the isotropic unflavored limit is rapidly attained in the IR.

UV limit
Let us consider the UV limit in whichd ⊥ is small and we have the following approximate relation between d ⊥ and x 0 :d (4.28) In this limit the potential can be approximated as: The integral can be computed analytically and yields the following result for the potential: In terms ofd ⊥ we get: which is equivalent to the result obtained in [21] for the massless background. As illustrated in Fig. 8, the potential (4.31) nicely matches the numerical results.

Entanglement entropy
The entanglement entropy of a region and its complement is a good measure of the quantum correlations of the system. In holography the entanglement entropy is obtained by minimizing an area functional for an eight-dimensional surface embedded in ten-dimensional spacetime [37,38]. Let A be a spatial region in the gauge theory. The holographic entanglement entropy between A and its complement is: where Σ is the eight-dimensional spatial surface whose boundary is A and minimizes S A , 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 Einstein frame. In this section we will apply this prescription when A is a slab of infinite extent in the two cartesian directions and having a finite width in the remaining cartesian coordinate. Clearly, there are two cases to study, namely parallel and transverse slabs, which we analyze separately. We note that again we find striking similarity with the results in [19].

Parallel slab
First we consider the case in which A is an infinite slab with a finite width parallel to the layers, namely: We will parameterize Σ by a function x = x(x 1 ). The eight-dimensional induced metric on Σ is: where the prime now denotes derivative with respect tox 1 . If we integrate over all the coordinates except x, we get: The first integral derived from S is: where the subscript nought denotes that the corresponding quantity is evaluated at the minimal value x 0 of x. From this last equation we get: which can be integrated to givel : The entanglement entropy for this configuration is given by the divergent integral: We will regularize S by subtracting the entropy of a configuration that we call the flat surface in the following. To conform with the homology constraint, the surface is not disconnected, but it is connected at the bottom x = 0. The full flat surface consists of three constant pieces: two straightx 1 = ±l /2 ones and a horizontal one x = const. = 0, each individually being solutions to the equation of motion. It turns out that the contribution of the horizontal surface to the area integral vanishes. The entropy for the flat embedding is thus Thus, we define the finite entanglement entropy as:  After some calculations, we get: The numerical results for S versusl are presented in Fig. 9 for several values of the parameter x q . In this plot we notice that S becomes positive whenl is large enough. According to our regularization procedure (5.10) this means that the disconnected surface is dominant with respect to the connected one whenl >l c , wherel c is a critical length which grows with x q . On the contrary, for smalll the entropy behaves as S ∼l − 4 3 , with a coefficient independent of x q . As pointed out in [21], this universal behavior is the one corresponding to an effective D2-brane and can be obtained analytically, as we show in the next subsection. A similar behavior of the EE has previously been obtained in backgrounds dual to confining theories and unquenched flavors, see [39][40][41].

UV limit
Whenl is small, the minimal value x 0 of x is large and we can use the expansion of the functions of the background valid for large x. At leading order, we get: Similarly, the regulated entanglement entropy for the parallel slab is: which can be integrated analytically with the result: (5.14) If we eliminate x 0 and write the entanglement entropy in terms ofl , we get: which is the same result as in [21] and, as shown in Fig. 9, matches perfectly the numerical results whenl is small.

Transverse slab
We now consider a slab with finite width in the direction ofx 3 . The region A in this case is: 16) and the induced metric on Σ becomes: where now x = x(x 3 ). The transverse entropy functional is: Now the first integral is: from which it follows that: Therefore, the transverse lengthl ⊥ is: The entropy functional evaluated on the minimal surface for this configuration is: whose divergent part is: We define S f inite ⊥ as: After some calculations one can demonstrate that: as functions of their corresponding rescaled width. In both panels curves with the same color correspond to the same value of x q : red for x q = 0.1, blue for x q = 0.5, and brown for x q = 0.9.
In Fig. 10 we plot S ⊥ as a function ofl ⊥ obtained by the numerical computation of the integrals in (5.25) and (5.21). For smalll ⊥ the curves for different values of x q coincide and behave as S ⊥ ∼l −6 ⊥ . This behavior is found analytically in the next subsection. For largê l ⊥ the dominant configuration is the disconnected one and S ⊥ becomes positive. We have also compared in Fig. 10 the parallel and transverse entanglement entropies. In the UV the difference is significant, but at larger distances the two entanglement entropies become very similar. This similarity is more and more pronounced as x q → 1.

UV limit
We now evaluate the entropy in the limit in which x 0 is very large andl ⊥ is very small. Using the UV expansion at leading order of the different functions of the background, we get: is: which, after performing the integral becomes: Eliminating x 0 in favor ofl ⊥ we reproduce the result of [21]: . (5.29)

Flow of mutual information
The mutual information of two entangling regions A 1 and A 2 is a measure of the information shared by these two domains and is defined as: We will analyze the evolution with the intrinsic scale of the background of I(A 1 , A 2 ) when A 1 and A 2 are two slabs of equal length l parallel to each other which are separated by a distance s. In holography there are two possible surfaces contributing to the entanglement entropy of two slabs [42,43]. One of these configurations, which has zero mutual information, dominates when the slab separation is large. Below some critical separation the second configuration is dominant and I(A 1 , A 2 ) becomes positive. In reduced units, the critical separationŝ for two strips of lengthl is determined by the vanishing of the mutual information: 2 S(l) = S(2l +ŝ) + S(ŝ) . where a = 4 3 (a = 6) for parallel (transverse) slabs. The solution of (5.32) in these two cases is: ŝ Let us next study the critical point in the long distance IR region, where we expect to approach the conformal isotropic behavior of AdS 5 × S 5 . In this last case the critical separation is also determined by (5.32) with a = 2 and, thus: Interestingly, viewed as an inverse it equals the metallic meanl/ŝ = σ p=1,q=1/2 = √ 3+1 2 [44,45]. Outside the UV region the criticalŝ/l depends onl. We have numerically verified that this flow behaves as expected. For parallel (transverse) slabsŝ/l grows (decreases) from its UV value (5.33) asl is increased and it actually becomes very close to the conformal isotropic value (5.34) for largel and x q close to one (see [45] for another recent example of a holographic flow of mutual information). We note that for very large distancesl, the dominant phase is that for flat embeddings and the full phase diagram resembles that of [43].

Thermodynamics of a massless probe brane
In this section we test our background with a probe D5-brane, embedded as in the array (1.1), in which we switch on a worldvolume gauge field A = A 0 dx 0 dual to a chemical potential µ. Our goal is to study the zero-temperature thermodynamics of the probe and its evolution as we move from the UV to the IR. For simplicity we will consider massless embeddings in which the probe reaches the IR end of the space. These type of embeddings have been analyzed in detail in App. C, where we check that the equations of motion of the probe are satisfied if the worldvolume gauge potential A 0 satisfies (C.21). In this section we will work directly in the radial coordinate ζ, for which the Lagrangian density takes the form:L where L is the Lagrangian density (C. 19). In (6.1) T is the constant defined in (C.20) and we have writtenL in terms of the master function W . The chemical potential is just the value of A 0 at the UV boundary. In our variables: where d is the integration constant of (C.21) which, as we check below, is proportional to the charge density. The grand potential Ω is given by minus the on-shell action, which in our case is finite and there is no need of regularizing it. This is due to the cancelation of the divergences between the DBI and WZ terms. Removing the Minkowski volume factor, we get: The charge density ρ can be written as: From (6.2) and (6.3) we can compute the derivatives with respect to d that are needed to calculate ρ: Clearly, one has: and we get that, indeed, ρ is related to d as expected: The energy density is given by: = Ω + µρ . (6.8) Plugging the values of Ω and ρ into (6.8), we get: Therefore: as expected. Taking into account that p = −Ω, the speed of sound u s can be obtained as: Thus, we get the following expression of u 2 s : More explicitly, plugging into (6.12) the expression (6.2) of µ, we obtain that u 2 s is given by the following ratio of two integrals over the holographic coordinate: In the unflavored (W = 1, b = 0) and massless flavored (W ∝ ζ − 2 3 , b = 0) cases we get the following d-independent results: (6.14) The unflavored result u 2 s = 1/2 is the one expected for a conformal worldvolume theory in 2 + 1 dimensions. In general, one should get a value depending on d, which interpolates between these two values. In order to facilitate the numerical calculations, let us rewrite these results in x coordinate. Recall that ζ and x are related as: It turns out that Q f and C can be scaled out from our formulas. First of all, we define the rescaled density and chemical potential as: Then, we have:μ whereŴ was defined in (3.4). Moreover, Ω and can be recast as: whereN is defined as:N The speed of sound can either be written as: 20) or if we define the integrals I 1 and I 2 as: then u s is obtained from the ratio between I 1 and I 2 : Recall that we can relate the quark mass of the background m q to the constant C, namely Using this relation we get thatd ∼ d/( Q f m 3 2 q ). Therefore, the UV limit m q → 0 corresponds to takingd → ∞. Accordingly, we should recover in this larged limit the massless flavored result of (6.14): , (UV limit). (6.23) The UV limit written above can also be analytically verified. Indeed, let us introduce a new integration variable y, related to x as: The minimal value of y, corresponding to x = 0, is y = y q , where: It is now straightforward to write I 1 and I 2 as: Whend → ∞ the lower limit of the integrals becomes y q = 0. Moreover, the argument of the functionŴ is large and one can use its UV asymptotic expression: W (x(y)) ≈ 9 8 Then for larged: ∞ 0 dy y 13 3 (1 + y 4 ) and therefore we get the expected UV result: Let us now explore the smalld limit. By takingd = 0 in (6.21) we see that the integrals I 1 and I 2 are proportional (I 2 (d = 0) = Z 4 q I 2 (d = 0)). Therefore, we get: This result is not the one we naively expect sinced → 0 corresponds to m q → ∞ and, as the quarks are not dynamical in this large mass limit, it would seem that we should get the conformal result u 2 s = 1/2 in this IR limit, at least when x q is close to one. In order to clarify the situation, let us taked = 0 in the integrals (6.21) and examine the behavior of the integrand near x = 0. For x q = 1 we get: where we have used (3.5) to obtain the behavior of the integrand near x = 0. This x q = 1 integral is convergent and the two integrals I 1 and I 2 are well defined. On the contrary, for x q = 1 we haveŴ = 1 inside the cavity and the integrand is divergent whend = 0 and thus we cannot take directlyd = 0 in the integrals. Actually, in the calculation of u s , the limitŝ d → 0 and x q → 1 do not commute. By taking x q → 1 first we indeed get the conformal result u 2 s = 1/2 independently of the value ofd. The numerical values of u 2 s obtained from (6.22) as a function ofd for different values of x q have been plotted in Fig. 11. In this plot we notice that the UV asymptotic result (6.29) is satisfied for all values of x q and for x q ∼ 1 there is a minimum at lowd, in which u 2 s approaches the conformal value u 2 s = 1/2. This is consistent with the behavior found in the analysis of other observables: the UV behavior is independent of x q and given by the scaling solution, whereas the IR is controlled by x q . By taking x q close to one, the long distance behavior of our system becomes more isotropic.

Summary and conclusions
The goal of this paper has been the construction of a gravity dual to a system containing multiple (2 + 1)-dimensional layers in a (3 + 1)-dimensional ambient theory. In order to deal with this problem we adopted a top-down holographic approach and used brane engineering to generate the corresponding gravity dual. We considered the setup (1.1), in which D3-and D5-branes intersect along 2 + 1 dimensions and a codimension-one defect is created along the worldvolume of the D3-branes. Moreover, the D5-branes are distributed homogeneously along the gauge theory directions orthogonal to the defect, giving rise in this way to a multilayer system. To find the corresponding supergravity solution, we regarded the D5branes as flavor branes and used the techniques developed to find the backgrounds dual to unquenched smeared flavor.
The background found is supersymmetric and solves the supergravity equations of motion with D5-brane sources. It is given in terms of the master function W , which can be obtained by solving the master equation (2.9). This master equation contains the profile function which can, in principle, be arbitrarily chosen and depends on the particular distribution of the D5-brane charge along the holographic coordinate. The solutions corresponding to a constant profile function p were obtained in [20]. In the flavor language, the solutions of [20] are dual to models with massless flavors living on the defect (the corresponding black hole was constructed and analyzed in [21]). The corresponding field theories display anisotropy in the third direction which, by construction, is produced by the multiple (2+1)-dimensional layers.
Here we generalized the scaling solution found in [20] to the case in which the quarks are massive and there is a region in the bulk, which we called the cavity, in which the flavor sources vanish. The size x q of this cavity provides us with a parameter which determines the long distance behavior of the model. Indeed, as we tuned x q to its maximal value x q = 1, the IR behavior of several observables we analyzed approaches the one of the un-layered theory, whereas the short distance behavior is independent of x q and given by the scaling solution. Thus, as x q → 1 the theory in the IR becomes more isotropic and effectively retains its (3 + 1)-dimensional character.
We studied the running anisotropic behavior described above in several quantities, but it is clear that we have not exhausted the list of observables to analyze. Let us mention some possible extensions of our work. In Sec. 6 we explored our background with probe D5-branes with a worldvolume gauge field dual to a chemical potential. For simplicity we considered embeddings of the probe with trivial embedding function χ, corresponding to massless flavors. The more general case of massive embeddings can be readily obtained by considering a general function χ(r). It would be interesting to see how the quantum phase transition studied in [46][47][48] between Minkowski and black hole embeddings is modified by the backreaction as has been seen in other (2 + 1)-dimensional systems [49].
It is clearly an interesting task to find the condensed matter system for which our holographic model could offer a framework for concrete calculations. Each layer effectively mimics graphene, when some of the D5-branes blow up in to D7'-branes [50][51][52][53][54][55][56], so the full multilayer system could perhaps be understood as a (gapped) holographic graphite. More study is needed to make this identification precise. Nevertheless, some of the recent multilayered systems are known to have strong coupling dynamics [57], so in an ideal scenario, we would love to engineer the holographic geometry to mimic the physics in these settings. This is, however, a highly non-linear task which requires novel ideas.
Fortunately, we have a rather straightforward avenue ahead of us to find a killer app. In [58] it was demonstrated that the original D3-brane background (the dual to N = 4 SYM theory), supplemented with flavor D7-brane probes (introducing quenched quark matter) at finite chemical potential, provides a realistic equation of state for cold and dense quark matter. This realization paved the road for further investigations [59][60][61][62].
How much then can holography help in understanding the deep cores of neutron stars where deconfined quark matter could reside? In the current context we would like to ask the following question. Black holes are known to forget about their past and are characterized by a handful of parameters (mass, rotation, charges). Neutron stars, on the other hand, seem to depend on a variety of parameters through their complicated internal composition. Surprisingly, however, certain (dimensionless) macroscopic properties (tidal deformability, quadrupole moment, moment of inertia) were found to obey universal relations to quite high degree of accuracy [63] independent of the underlying equation of state. 2 A recent interesting paper [65] drew attention to the resemblance of the neutron star universal relations and the no-hair relations of black holes. In order to formally study the limit of strong gravity and the black hole formation, the equation of state has to be anisotropic due to Buchdahl bound [65]. We believe that the rather theoretical construction laid out in this paper (smeared D5-brane providing the necessary anisotropy) is a key step towards this direction and we hope to report progress on this front in the near future.
left-invariant SU (2) one-forms satisfying dω i = 1 2 ijk ω j ∧ ω k . Then, we can write ds 2 10 as: where the fiber τ takes values in the range 0 ≤ τ ≤ 2π. Notice that we are using in (A.1) a radial variable r which is different from the one in (2.1) (see below for the relation between r and ζ). The one-form A in (2.1) is given by: The one-forms ω i can be represented in terms of three angles (θ, ϕ, ψ) as: where 0 ≤ θ ≤ π, 0 ≤ ϕ < 2π, and 0 ≤ ψ < 4π. The coordinate representation of the canonical vielbein basis of C P 2 is: By plugging these expressions into (2.6) and (2.7) we get the value of the RR three-form F 3 written in (2.7), where we already used the new radial variable ζ defined in (A.10) below. Our solution preserves some amount of supersymmetry. Indeed, let us choose the following vielbein basis for the ten-dimensional metric (A.1): Then, in this basis of one-forms, the Killing spinor of the background can be written as: where η is a constant spinor satisfying some projection conditions. If we represent the Killing spinors of type IIB supergravity as a Majorana-Weyl doublet, then η satisfies the projections: which means that our background preserves two supercharges. Actually, the different functions of our ansatz must satisfy the following system of first-order BPS equations [20]: where p(r) is the profile function of (2.7), the prime denotes derivative with respect to the radial coordinate r and Q f and Q c are related to the number of flavors N f and colors N c as: where L 3 = dx 3 . Let us now write the BPS system in terms of a new radial variable ζ, related to r as: Then, the equations for φ, g, and f take the form: In this new variable ζ, the BPS equation for g in (A.11) can be immediately integrated: and we can rewrite the remaining BPS equations as: Notice that our ansatz for the metric in the ζ variable is precisely the one written in (2.1). We now introduce the new master function W as: (A.14) One can easily show that the BPS system (A.13) implies that W satisfies the following first-order differential equation: Moreover, one can demonstrate that the equation for φ in (A.13) can be written as: One can combine these last two equations in the following second-order master equation for W : which coincides with (2.9). It is also straightforward to verify that the function f and the dilaton φ are given in terms of W as in (2.10) and (2.11), respectively. As mentioned in Subsec. 2.3 the master equation can be solved in powers of ζ q /ζ for large ζ. The result for W was written in (2.27). The function f is readily computed using (2.10): The dilaton in this expansion is: where C 1 is a constant of integration. The squashing factor q = e f /ζ for this solution can be obtained from (2.13) and takes the form: We can also solve the master equation (2.9) perturbatively close to the cavity for ζ ≥ ζ q . In order to do that, let us substitute in the master equation an expansion of the form: Identifying the above approximate solution (A.22) and the analytic solution (2.16) at the edge of the cavity ζ = ζ q we get the following two constraints (A.23) The first non-vanishing coefficients are The function f corresponding to (A.22) is: 26) and the dilaton can be expanded as: Finally, plugging these expansions into (2.14) we can solve for the warp factor h, which is given by The integration constant C 1 is connected with the integration constant coming from solving the warp factor differential equation inside the cavity.

B Probe analysis and profile function
In this appendix we study the embeddings of a D5-brane probe which preserves the supersymmetry of a background given by our ansatz. Once these embeddings are characterized we will be able to obtain the corresponding profile function p(ζ) which, in turn, is a necessary ingredient to solve the BPS equations and determine completely the different functions of the ansatz. The supersymmetric embeddings we are looking for must satisfy the kappa symmetry condition: where is a Killing spinor of the background. For a D5-brane without any worldvolume gauge field, Γ κ is given by: where g 6 is the determinant of the induced worldvolume metric, γ α 1 ···α 6 is the antisymmetrized product of induced Dirac matrices and we are representing the Killing spinors of type IIB supergravity as a Majorana-Weyl doublet. We will take ξ α = (x 0 , x 1 , x 2 , r, θ, ψ) as our set of worldvolume coordinates. In this case the kappa symmetry matrix takes the form: Recall that the Killing spinor of the background can be written as in (A.6) in terms of a constant spinor η which satisfies the projections (A.7). We can cast the kappa symmetry condition (B.1) as the following condition on η: whereΓ κ is defined as the matrix: We will consider an embedding in which x 3 and ϕ are constant, while The induced γ-matrices on the worldvolume for this embedding ansatz are: where we have denoted and the Γ's are constant ten-dimensional Dirac matrices. From these induced matrices we can compute the antisymmetrized product γ rθψ and get an expression of the type: where the different coefficients are given by This expression can be rewritten as: where the coefficients d i are: In order to compute the form of Γ κ , we use that We find Let us now calculate the matrixΓ κ . As and, therefore: We now study how the different terms ofΓ κ act on η. We begin by considering the action of the first term on the right-hand-side of (B.17). First of all, we notice that, using the ten-dimensional chirality condition satisfied by and the projections (A.7), one can show that η satisfies [20]: Using (B.18) and the first projection in (A.7), we get: Finally, the last projection in (A.7) allows us to write: Let us next consider the second term inΓ κ . We first use: from which it immediately follows that: Similarly, since Γ r15 = −Γ 35 Γ r13 , we get that the third term of (B.17) acts on η as: But, according to the last set of projections in (A.7), we have Γ 35 η = −Γ r4 η. Therefore: Finally, as Γ 415 = −Γ r4 Γ r15 , we get: Thus, collecting all these results, we have: AsΓ κ should act as (plus or minus) the identity matrix on η, we must have: The first equation fixes that ψ + 3τ = n π , (B.28) with n ∈ Z. This implies thatτ In order to have τ having values within the range [0, 2π], we choose n = 4 in (B.28). Therefore, the function τ = τ (ψ) is: and the values of τ for the embedding range from τ = 0 to τ = 4π 3 . Moreover, the second equation in (B.27) implies the following first-order equation for χ: (B.31) By making use of (B.29), this equation can be written as: The induced metric on the D5-brane worldvolume for our ansatz is: which shows that our brane configuration preserves the same amount of SUSY as the background. This equation can be integrated as:

B.1 General integration of the BPS equation
where c is an integration constant. By inspecting (B.40), it is clear that the coordinate ζ of the brane must be greater or equal than some minimal value. Let us denote by ζ q this minimal value of ζ. Clearly, ζ q and the constant c are related as: or, equivalently: 3 In terms of ζ q the embedding angle χ is given by: It follows from (B.44) that χ vanishes at the tip of the brane ζ = ζ q and grows as we move towards the UV, until it reaches an asymptotic value χ = χ * , with χ * given by: Notice that χ = χ * is a constant angle solution of the BPS equation (B.39). For this χ = χ * embedding, the tip of the brane is at ζ q = 0 (see (B.44)). Therefore, the brane reaches the origin of the geometry, as it corresponds to having massless quarks. 3 We are implicitly assuming that c 2 is positive. If we take c 2 < 0 in (B.40) the minimal value ζ q of the coordinate ζ is achieved when cos χ = −1 and is given by ζ 2 q = − 3 4 c 2 . Thus, the embedding function in this second branch can be written as: For these embeddings χ = π at the tip of the brane and decreases when we move towards the UV region ζ → ∞, reaching the value χ = χ * asymptotically.

B.2 The profile function
To determine the function p(r) for a set of branes ending at the same distance from the origin, i.e., with the same mass, we compare the smeared and localized brane actions. Let us start with the WZ term, for which the smeared distribution is: where C 6 is the RR six-form potential given by: with K (6) being the calibration six-form written in [20]. In terms of the one-forms E a of our vielbein basis (A.5), K (6) is given by: where Ω is the following complex three-form: In (B.46) Ξ is the smearing four-form for the massive case, dF 3 = 2κ 2 10 Ξ, which, according to (2.8), is Therefore: Integrating in (B.51) over the angular variables and x 3 for − L 3 2 ≤ x 3 ≤ L 3 2 , we get: where the smeared WZ Lagrangian density is: We now compare (B.53) with the WZ action of a localized single D5-brane multiplied by N f . As all the flavor branes that form the set of our smeared distribution are identical, these two quantities should be equal. The WZ term of the action of a D5-brane without worldvolume gauge field is: where the hat denotes pullback to the worldvolume and C 6 has been written in (B.47) in terms of the calibration form of (B.48). Let us suppose that we take ζ α = (x 0 , x 1 , x 2 , r, θ, ψ) as worldvolume coordinates and consider an embedding of the type (B.6) with x 3 and ϕ constant. Then, the pullback of C 6 is: 16 cos Let us next assume that τ (ψ) is given by (B.30) and let us integrate the action over θ and ψ. The resulting Lagrangian density (to be integrated over (x 0 , x 1 , x 2 , r)) multiplied by N f is: Let us now compare the terms in (B.53) and (B.56) without e f . They are equal provided the profile p(r) is given by: When r → ∞ the angle χ → χ * , where χ * is the value written in (B.45). Since cos χ * 2 sin χ * = 4 3 √ 3 , we have: As argued on general grounds in Sec. 2.3, we should have p(r → ∞) = 1, i.e., the profile should approach the one corresponding to massless quarks in the far UV. This only occurs when Q f and N f are related as: This same conclusion can also be reached by comparing (B.53) for p = 1 and (B.56) for χ = χ * . It follows that the function p(r) is related to the function χ(r) for our fiducial embedding as: Moreover, by equating the terms proportional to e f in (B.53) and (B.56) we conclude that p(r) should satisfy the differential equation: By computing the derivative of the right-hand side of (B.60) one can easily show that indeed p satisfies (B.61).
Let us now check the equality of the DBI terms in the smeared and microscopic approach. The smeared DBI action is: where the smeared DBI Lagrangian density L smeared DBI is: Let us compare the action (B.62) to the DBI action of the fiducial brane multiplied by N f , which is equal to: Taking into account (B.36), we get that (B.63) should be equal to: ( which can be shown to be equivalent to (B.61) after using (B.32). Let us finally obtain the explicit form of the profile in terms of the radial variable ζ. First of all, we write p as: Using (B.44), we get: 4 Notice that p(ζ) is continuous at ζ = ζ q , as it should. 4 For the second branch of embeddings with fiducial embedding function given by (B.42), the profile p(ζ) can be obtained by substituting (B.42) into (B.67), with the result:

C Equations of motion of the probe
The DBI action, in Einstein frame, of a D5-brane with a worldvolume gauge field F turned on is: We are interested in having a gauge potential A dual to a U (1) charge density. Therefore, we take for which F = dA is given by: Let us choose ξ α = (x 0 , x 1 , x 2 , r, θ, ψ) as our set of worldvolume coordinates and let us assume an embedding ansatz as in (B.6), i.e., with χ = χ(r) and τ = τ (ψ). The Lagrangian density L DBI corresponding to the action (C.1) is: (C.4) The WZ term of the action has been written in (B.54). Taking into account (B.55), we can easily demonstrate that L W Z takes the form: cos χ 2 cos(3τ + ψ) 2 sin χ + e f (1 + cos χ + 4τ ) χ .

(C.5)
Let us now analyze the equations of motion of χ(r) and τ (ψ) derived from the total Lagrangian L = L DBI + L W Z . The equation of motion of χ(r) derived from the total action is: The derivatives with respect to χ appearing in (C.6) are: Moreover, if we define the auxiliary functions J DBI (χ) ≡ cos χ 2 sin 2 χ + e 2f −2g (1 + cos χ + 4τ ) 2 J W Z (χ, χ ) ≡ cos χ 2 2 sin χ + e f (1 + cos χ + 4τ ) χ , (C.8) the derivatives with respect to χ we need in (C.6) are: Let us analyze the dependence on the angular variable ψ of the different terms in (C.6). The only terms that depend on ψ in (C.6) are those that containτ and 3τ +ψ. By inspecting how these terms enter into (C.7) and (C.9), one easily concludes that the equation of motion of χ can only be satisfied ifτ and 3τ + ψ are constant. This last condition implies thaṫ τ = −1/3, as in (B.29). Let us now look at the equation of motion for τ (ψ): The remaining τ derivative that we have to compute is: cos χ 2 sin(3τ + ψ) 2 sin χ + e f (1 + cos χ + 4τ ) χ , (C.14) which vanishes independent of χ if sin(3τ + ψ) = 0, i.e., when τ depends on ψ as in (B.28). As in the supersymmetric solution, we will take n = 4 in this equation. The conclusion of this analysis is that, for the type of ansatz we are considering, the function τ (ψ) must be given by (B.30) in order to satisfy the equations of motion of the probe D5-brane. We will now study these equations separately in two different cases.

C.1 The BPS solution
Let us consider now the BPS configuration in which τ (ψ) is given by (B.30), A 0 vanishes and χ(r) satisfies the first-order differential equation (B.32) dictated by kappa symmetry. By using (B.35) to evaluate the different derivatives of (C.7), (C.9), and (C.14) one can easily show that: which implies the fullfilment of the equations of motion for the BPS configuration.

C.2 The massless solution
Let us consider a massless embedding of the probe brane with non-zero chemical potential. Therefore, we will try to solve the equations of motion with χ constant. It is clear from the first equation in (C.7) that ∂L DBI /∂χ vanishes if χ is zero. Moreover, from the second equation in (C.7) we conclude that ∂L W Z /∂χ is zero if χ = χ * , where χ * is the angle defined in (B.45). Moreover, since 16) it follows that χ = χ * solves the equations of motion of the probe. Let us next define J * as: It follows that: It remains to satisfy the equation of motion for A 0 . The Lagrangian density for the gauge field is: where the effective tension T is given by: Since A 0 is a cyclic variable in the Lagrangian density (C. 19), the equation of motion for A 0 can be integrated once, giving: where d is a constant proportional to the charge density. From this equation we get: The chemical potential µ is just the value of A 0 at the boundary, and is given by the following integral: (C.23)

D Probes in the Higgs branch
We now study embeddings of a probe D5-brane in which the coordinate x 3 is not constant but instead bends as the holographic coordinate changes. This bending can be interpreted as a recombination between the D3-and D5-branes, realizing the Higgs branch of the dual theory, see, e.g., [66]. In order to find this configuration let us consider the following set of worldvolume coordinates ξ α = (x 0 , x 1 , x 2 , r, θ, ψ) and the following embedding ansatz: In addition, the probe D5-brane will have a worldvolume flux, given by: with Q being a constant. We show below that, if certain BPS equations are satisfied, this configuration preserves the supersymmetries of the background. The induced metric for the ansatz (D.1) is: 3) The DBI action in Einstein frame is given by (C.1). In the present case the DBI determinant is: The WZ action now takes the form: where the WZ Lagrangian density is given by: and J W Z (χ, χ ) has been defined in (C.8).
Let us now study the equations of motion for z(r). By inspecting L DBI and L W Z we immediately conclude that they do not depend on z (only on z ) and, therefore, z(r) is a cyclic variable. Thus, the equation of motion for z(r) can be integrated once as: We will consider the case in which the constant on the right-hand side of (D.8) is zero which, as we will verify below, is the supersymmetric configuration. Therefore, we must have: Moreover, the derivatives of L DBI and L W Z with respect to z are: A useful relation that can be derived from this last equation is: Let us now derive the equations of motion for the embedding function χ(r). First of all, we calculate the derivatives of the DBI Lagrangian density with respect to χ and χ, which are given by: 1 + e 2g 4 (χ ) 2 + h −1 e −2φ (z ) 2 J 2 DBI (χ) + 64 h −1 e −4g−φ Q 2 J DBI (χ) ∂J DBI (χ) ∂χ . (D.13) Let us now use (D.12) to compute the square root of the right-hand-side of (D.13) involving χ . We get: J DBI (χ) 1 + e 2g 4 (χ ) 2 1 + e 2g 4 (χ ) 2 ∂J DBI (χ) ∂χ , (D.14) which are just the derivatives corresponding to the case Q = z = 0 with A 0 = 0 (see (C.7) and (C.9)). As the derivatives of L W Z with respect to χ and χ are the same as in the Q = z = 0 case, we immediately conclude that the equation of χ is satisfied if τ (ψ) and χ(r) fulfill (B.28) and (B.32) respectively. Thus, the addition of internal flux, related to the bending of the brane in the x 3 direction as in (D.11), does not modify the BPS solution for χ(r) and τ (ψ).
Let us now write an explicit equation for the bending function z(r). This equation can be obtained by plugging into (D.11) the relation: 5 1 + e 2g 4 (χ ) 2 = J DBI sin χ cos χ 2 .

D.1 Kappa symmetry
In this subsection we will derive the first-order equations for τ (ψ), χ(r), and z(r) from the kappa symmetry condition of the probe. We begin by computing the induced Dirac matrices 5 Another useful relation is 15) for the ansatz (D.1): The kappa symmetry matrix with F θψ flux in our conventions is: with F θψ = Q. In the second term in (D. 19) we need to compute γ θψ γ θψ = g θθ g ψψ (γ θψ ) 2 , (D.20) with g θθ and g ψψ being elements of the inverse induced metric (D.3). In order to calculate the square of γ θψ , let us write this matrix as: γ θψ = A e −ψ Γ 12 B Γ 13 + C Γ 15 , (D.21) with A, B, and C being given by: where, in the last step, we have used that (Γ 13 ) 2 = (Γ 15 ) 2 = −1 and that {Γ 13 , Γ 15 } = 0. Moreover, by inspecting (D.3) we can write g θθ and g ψψ as: It follows that: γ θψ γ θψ = −1 .
(D. 25) Thus, Γ κ can be written as the sum of two terms: where Γ (1) κ and Γ (2) κ are given by: Let us now proceed as in App. B and compute the antisymmetrized product γ rθψ by using (D.18). We get: where the coefficients d 1 , d 2 , d 3 , and d 4 are given by the same expression as in (B.12) and the new coefficients d 5 and d 6 are: Let us next define the rotated kappa symmetry matrixΓ κ as in (B.5). Then, the kappa symmetry condition is the one written in (B.4). Moreover, we can writeΓ κ =Γ (1) κ +Γ (2) κ , withΓ The rotated matrixΓ (1) κ can be written as: whereasΓ (2) κ is given by: Let us now calculateΓ (1) κ η. This calculation is similar to the one performed to derive (B.26). In order to compute the additional terms, we need to use the projections: Using these results, we arrive at: (D.34) Next, we computeΓ (2) κ η. We need the projections: and we obtain: and one can easily show that it is a consequence of the BPS equations for χ and z. Therefore, we can write: