Holographic Entanglement Entropy in Anisotropic Background with Confinement-Deconfinement Phase Transition

We discuss a general five-dimensional completely anisotropic holographic model with three different spatial scale factors, characterized by a Van der Waals-like phase transition between small and large black holes. A peculiar feature of the model is the relation between anisotropy of the background and anisotropy of the colliding heavy ions geometry. We calculate the holographic entanglement entropy (HEE) of the slab-shaped region, the orientation of which relatively to the beams line and the impact parameter is characterized by the Euler angles. We study the dependences of the HEE and its density on the thermodynamic (temperature, chemical potential) and geometric (parameters of anisotropy, thickness, and orientation of entangled regions) parameters. As a particular case the model with two equal transversal scaling factors is considered. This model is supported by the dilaton and two Maxwell fields. In this case we discuss the HEE and its density in detail: interesting features of this model are jumps of the entanglement entropy and its density near the line of the small/large black hole phase transition. These jumps depend on the anisotropy parameter, chemical potential, and orientation. We also discuss different definitions and behavior of c-functions in this model. The c-function calculated in the Einstein frame decreases while increasing $\ell$ for all $\ell$ in the isotropic case (in regions of $(\mu,T)$-plane far away from the line of the phase transition). We find the non-monotonicity of the c-functions for several anisotropic configurations, which however does not contradict with any of the existing c-theorems since they all base on Lorentz invariance.


Contents
Fundamental questions addressed in studies of high energy heavy ions collisions (HIC) at RHIC and LHC, and future experiments NICA and FAIR, concern understanding of quark-gluon plasma (QGP) formation, i.e. thermalization of media produced in HIC, thermodynamic entropy production, and its characteristics such as quantum entanglement, decoherence etc. Most of our knowledge on the formation and properties of QGP resulting from HIC is obtained from measurements of the yields and spectra of particles in the final state of colliding heavy ions and their thermo/hydrodynamic interpretation. According to our common understanding, within a short time of order 1 − 2 fm/c collision systems reach a state that can be approximated by a thermal medium located in an expanding ball. This medium is characterized by local thermodynamic parameters including temperature and entropy density.
HIC experiments at RHIC and LHC have provided strong evidence that this medium is a QGP at large temperatures and densities. There is also a strong evidence for the existence of the confinement/deconfinement phase transition in the (µ, T )-plane, i.e. temperature on chemical potential. The experimental search for the QCD phase transition is nowadays one of the central goals for current and future collider facilities [1].The experimental search is mainly related to the measurement of fluctuations of net-proton or net-charge multiplicity [2][3][4][5] which are expected to exhibit non-monotonic behavior near the phase transition. The proper understanding of the experimental results requires careful theoretical analysis of the dynamical processes taking place near the phase transition lines, especially near the critical endpoint (CEP) 1 [7]. There are several theoretical approaches to searches of the QCD phase transitions. One of the promising theoretical approaches is based on lattice calculations, but this approach encounters difficulties with non-zero chemical potential [8,9].
Holographic duality provides an alternative approach to the study of the QCD phase transitions [10][11][12][13][14] (and for review see [15][16][17][18][19][20] and refs therein). Holography allows to define thermal entropy and free energy as functions of temperature and chemical potential of the gravitational background. With this approach, the thermalization process is dual to the black hole formation, linking the study of thermalization in the QCD with the study of the black hole formation in 5-dimensional gravity. The thermalization process is conveniently monitored, at least theoretically, by tracking the evolution of entanglement entropy of a selected region 2 . The entan-glement entropy of the selected domain A is defined as the von Neumann entropy of the reduced density matrix obtained by tracing out the degrees of freedom located out of domain A. The entanglement entropy is hard to compute in QCD, but it is suitable to compute it in the lattice [23][24][25], see also [26]. The entanglement entropy can be also evaluated at the gravity side. It is called the holographic entanglement entropy (HEE) and is defined as the area of a minimal surface extending from some predefined surface A on the boundary into the bulk [27][28][29]. The HEE during thermalization usually evolves to the thermal entanglement entropy [30][31][32][33][34][35][36][37][38]. With this approach, there is a natural possibility of studying the evolution of entropy in HIC (thermalization) and phase transitions for the obtained thermal media in the framework of the same holographic model.
Starting from the Landau thermodynamic approach to high energy collision [39] one connects the entropy of the ball of QGP produced in HIC with the total multiplicity of particle production with HIC. The ideal hydrodynamics preserves the total entropy, but dissipative one does not. One can relate the entanglement entropy associated with a given region of the ball with the multiplicity of particles produced in this region during the HIC [40]. After thermalization the entanglement entropy of the area depends on whether this area belongs to Glauber's participant area of ions collisions or does not. It occurs that the dependence of the entanglement entropy on geometrical size of entangled areas is closely related to energy loss and jet quenching.
We start with the most general anisotropic holographic model and consider the general orientation of the slab-shaped entangled region with respect to the geometry of HIC. The natural coordinate system defined by the HIC is such that the first axis (the longitudinal axis) is directed along the line of collision and the second (the transversal axis) is determined by the direction of the impact parameter. The orientation of the entangled region can be set by Euler angles. Then, we calculate the HEE for an anisotropic model with symmetry in transversal directions studied in [60,61]. The choice of this anisotropic model is motivated by its previous use for the holographic description of HIC [60]. It is this model that for special parameter of anisotropy ν = 4.5 gives the dependence of the produced entropy on energy in accordance with the experimental data for the energy dependence of the total multiplicity of particles produced in heavy ion collisions [62]. Isotropic holographic models had not been able to reproduce the experimental multiplicity dependence on energy [63][64][65][66][67][68][69][70][71][72][73][74][75]. As shown in [76], the model [61] describes smeared confinement/deconfinement phase transitions. As we will see, the model indicates the relations of the fluctuations of the multiplicity, i.e. the entanglement entropy, with the phase transitions.
There are more reasons to use the anisotropic models in context of HIC [77] and QCD itself [78]. Other anisotropic holographic models have been actively studied in recent years [79][80][81][82][83][84], and the main motivation for these studies was also the anisotropic nature of HIC.
We take an entangled slab-shaped area that has a finite extent in one direction and infinite extent in the other two directions 3 . We implicitly assume that the slab of the interest is located inside the overlapping area of two ions. In this case of the spacetimes with the fully anisotropic metric (see (2.1) below) the problem of finding the extremal area functional effectively reduces to finding the geodesics in some auxiliary 2-dimensional Euclidean space. For the model [61] we find that varying the angle between the axis of collisions and the direction of the smallest side of the slab-shaped entanglement area changes the slab HEE. Increasing this angle we enhance the HEE, and according to conjecture [40], this means the enhancement of the multiplicity of particle production from this area. This enhancement depends on geometrical (length), and thermodynamical (temperature T and chemical potential µ) parameters. The HEE density [86][87][88][89] is a more convenient object for study since it does not suffer from ultraviolet divergencies. The HEE and its density undergo jumps and these jumps increase while increasing angle. Moreover, the values of the HEE density and its jumps increase with the anisotropy increasing.
We also discuss various definitions and behavior of c-functions in this model. In the isotropic case we use different frames, the Einstein frame or the string frame, as well as different types of renormalizations. Here we mainly use the geometric renormalization scheme, which consists in subtracting the disconnected configuration from the connected one [41,45]. For the isotropic case in regions of (µ, T )-plane far away from the phase transition line, the c-function calculated in the Einstein frame decreases while increasing for all . We find the non-monotonicity of cfunctions in the string frame, that is related to dynamics of the dilaton field in UV in our model. There are several proposal in the literature how to define the holographic c-function in the anisotropic backgrounds [90][91][92][93][94][95][96][97]. Here we use prescription [95] adapted to our renormalization schemes. We find the non-monotonicity of cfunctions for several anisotropic configurations and discuss their origin. Note, that generally speaking, the non-monotonicity of c-functions does not contradict any of the existing c-theorems [98][99][100][101][102][103][104][105]107] since they all base on Lorentz invariance.
The paper is organized as follows. In Sect. 2.1 we briefly describe anisotropic holographic models. In Sect. 2.1.1 we present the most general anisotropic holographic model. We present the action and the ansatz that solves the EOM for the anisotropic model with symmetry in transversal directions in Sect. 2.1.2 and thermodynamics of the background in Sect. 2.2. In Sect. 3 we present an expression for the HEE for a slab-shaped entangling region oriented differently with respect to the HIC axes. In Sect. 3.4.1 and Sect. 3.4.2 we present its special forms corresponding to transversal and longitudinal orientations in respect to the collision axis, respectively, and discuss the regularization procedure used for the transversal and longitudinal orientations. Sect. 4 is devoted to entanglement entropy density and definitions to c-functions. In Sect. 5 we display and discuss our main numerical results. In Sect. 5.1 we demonstrate the dependence of the HEE on the geometrical parameters -length , anisotropy ν and orientations, and on thermodynamical parameterstemperature T and chemical potential µ. In Sect. 5.2 we show dependencies of the HEE density on the thickness, anisotropy and orientation of the slab as well as on temperature and chemical potential. In Sect. 5.3 we discuss the scaling behavior of the modified c-function. In Sect. 5.4 we discuss possible origins of a non-monotonic behavior of c-functions. In Sect. 5.5 we present behavior of different c-functions near the background phase transition. In Sect. 5.6 we present a table of the dependencies of various c-functions on . In Sect. 5.7 we compare the position of the phase transition for HEE with the positions of phase transitions related to the background instability.
Finally, we end the paper with the conclusion and discussion of future directions of research on the subject. We start with a general anisotropic holographic model .
is the AdS deformation factor (in the string frame in the presence of the dilaton), g(z) is a blackening function and g i (z) are anisotropy factors.

Anisotropic Model with Symmetry in Transversal Directions
Taking in the previous formula g 1 (z) = 1 and g 2 (z) = g 3 (z) = g(z) we get 3) with a special form of anisotropy factor g(z) as in [35] and a particular case of b(z) describing the holographic model for heavy quarks [61] present a special interest for us: where L is the characteristic length scale of the geometry and b(z) = e cz 2 /2 . In all our calculations we set c = −1. Here the metric is in the Einstein frame. In the next section to perform calculation of the HEE we also switch to the string frame adding an extra dilaton-dependent exponential prefactor [27,41]. The metric (2.4) is supported by the Einstein-Dilaton-two-Maxwell action with special potential V for the dilaton field φ and strength potentials f 1 and f 2 for two Maxwell fields: where F (1)2 and F (2) 2 are the squares of the Maxwell fields. The ansatz with the Maxwell fields F µ = A t (z)δ 0 µ , and F (2) = q dy 1 ∧ dy 2 , φ = φ(z) and metric (2.4) satisfies the equations of motion under relations between the warp factor b(z), dilaton and Maxwell potentials V (φ), f 1 (φ), f 2 (φ) [61].
The ansatz (2.4) breaks isotropy while preserves translation and (t, x)-boost invariances. The advantage of this model is that it allows to find potentials V , f 1 and f 2 and the blackening function g explicitly. Note that in this case the dilaton potential can be approximated by the sum of two exponents: Note, that in [108] an explicit isotropic solution for the dilaton potential as a sum of two exponents and zero chemical potential has been constructed. The isotropic version of the model has beed considered earlier in [109]. There is another holographic anisotropic model that has symmetry in transversal directions, but breaks the boost invariance in (t, x 1 ) plane, keeping the (t, x 2 ) and (t, x 3 ) planes. It is based on the Mateos-Trancanelli metric [82,83] This metric is supported by the Einstein-Axion-Dilaton action.

Thermodynamics of the Background (2.4)
The thermodynamical properties of the anisotropic holographic model were studied in [61]. The temperature of the model is given as T = |g (z h )| /4π, where g(z h ) is given by the incomplete gamma function (see (2.37) in [61]). The expressions for thermal entropy and the free energy are the following: where z h 2 is the second horizon, i.e. at this point T (z h 2 ) = 0. For the considered blackening function the function T (z h ) is three-valued i.e. has the Van der Waals type of behavior at 0 < µ < µ cr (ν) (Fig.1.A)) and free energy shows the swallowtailed dependence on the temperature in this range of chemical potentials ( Fig.1.B)). The loop of the swallow-tailed shape disappears at (µ, T ) = (µ cr (ν), T cr (ν)). For µ ≥ µ cr , the curve of free energy increases smoothly from higher to lower values of temperature. The lowest values of free energy correspond to the thermodynamically stable phases. The line of free energy intersects itself at T = T BB (ν, µ), and here a small black hole transits to a large one. Note that for the non-zero values of chemical potential, the Hawking-Page (HP) transition occurs at the temperatures larger than the temperature of black hole to black hole transition (BB) (Fig.1.B)). So the black hole solution is always dominant with respect to thermal AdS for µ = 0. The position of BB phase transition line determines the phase diagram of the model [61,76]. The entropy function s(T ) is multivalued for µ < µ cr and becomes one-to-one for chemical potentials µ ≥ µ cr . The value of the µ cr depends on the anisotropic parameter ν, µ cr (ν = 1) = 0.117, µ cr (ν = 1.5) = 0.189, µ cr (ν = 2) = 0.270, µ cr (ν = 3) = 0.299, µ cr (ν = 4.5) = 0.349. In Fig.2.A) the BB phase transitions in the (µ, T )-plane are presented for the isotropic and anisotropic backgrounds with different ν. At the temperature values T = T BB (ν, µ) the thermal entropy undergoes a significant jump (see Fig.2.B)), where the entropy is presented in the logarithmic scale. On these plots we see that the jumps disappear for µ ≥ µ cr (µ cr = 0.117 for ν = 1 and µ cr = 0.349 for ν = 4.5).

General Framework
The entanglement entropy is used to probe correlations in quantum systems. If the system is divided into two spatially disjoint parts A andĀ, the entanglement entropy S(A) gives an estimation of the amount of information loss corresponding to the restriction of an A. It is not to be simple to calculate the entanglement entropy from the strongly coupled system side, in particular in QCD. However, one can compute its holographic dual. For some boundary region A the HEE is obtained by extremizing the 3-surface functional that ends on the boundary surface A. In the dual field theory the entanglement entropy of a subsystem A is given by the formula [27][28][29] In what follows we set G 5 = 1. From (3.1) we see that the entanglement entropy depends on the geometry of the area A. It is difficult to do calculations for arbitrary A, compare for example with [85]. We will do the calculations for the areas having the shape of parallelepipeds, two sides of which are long and one short, i.e. for parallelepipeds which look as slabs.
The orientation of these parallelepipeds can be specified by the Euler angles (φ, θ, ψ) relative to the axes x 1 , x 2 , x 3 specifying the HIC, see Fig.4. The axis x 1 is taken along the collision line, the axis x 2 is taken in the transversal direction along the direction of the impact parameter b, and the axis x 3 is taken along the emerging magnetic field, see Fig.5. In Fig.4 the initial slab is oriented along the axes specified by the HIC geometry and is shown in green. The rotated slab shown in pink defines the entanglement area. We are going to calculate the HEE for the rotated pink parallelepiped.
The entangled slabs are supposed to be in the 3-dimension regions of overlapping of two nuclei regions, so called Glauber regions, see Fig.3. The overlapping region of two nuclei depends on time. At a fixed moment in time it is a three-dimensional body in which the cross section perpendicular to the axis of collisions has a shape bounded by arcs of two circles shifted relative to each other by an impact parameter. The area and the shape of this area depend on the impact parameter b. The overlapping area for a peripheral collision is approximated by the parallelepiped with sizes L 1 , L 2 and = D − b. We assume L 1 , L 2 >> , each ion is presented as a disk with radius D/2 and the impact parameter b is essentially less than D, << D. The overlapping area of two ions is almost cylinder and we can consider the case of very short time after collision, << D, L 1 , L 2 ≈ D. There are two limited cases: peripheral collision and central one, see Fig.3.
Let us show that for the spacetimes with the metric (2.1) the problem of finding the extremal area functional (3.1) for the slab with an arbitrary orientation effectively reduces to finding geodesics in some auxiliary 2-dimensional Euclidean space. To show this we consider the embedding in the static gauge and assume that orientation of the slab in respect to the HIC axes is given by the Euler angles, see Fig.4. Since there is no dependence of the integrand on ⇠ 2 , ⇠ 3 we can perform the integration over these variable that give the sizes of the parallelepiped in second and third directions S = L 2 L 3 f actor?
Here`, L 2 , L3 are lengths of the parallelepiped in the first, second and third direc- The Nambu-Gato action is Since there is no dependence of the integrand on ⇠ 2 , ⇠ 3 we can perform the integration over these variable that give the sizes of the parallelepiped in second and third directions S = L 2 L 3 f actor?
Here`, L 2 , L3 are lengths of the parallelepiped in the first, second and third directions. The action (3.17) is a particular case of the Born-Infield action  The Nambu-Gato action is We have 3) x 4 (ξ) = z(ξ 1 ), Or in term of elements of the Euler matrix we get: (a 21 a 23 a 22 a 31 ) 2 g 1 g 2 (3.13) + (a 31 a 32 a 21 a 33 ) 2 g 1 g 3 + (a 23 a 32 a 22 a 33 ) 2 g 2 g 3 ⌘ z 0 2 g (3.14) CHECK (see FulAniz-NG-action.nb ) g 22ḡ33 ḡ 2 23 = (g 1 a 2 21 + g 2 a 2 22 + g 3 a 2 23 )(g 1 a 2 31 + g 2 a 2 32 + g 3 a 2 33 ) (3.15) (g 1 a 21 a 31 + g 2 a 22 a 32 + g 3 a 23 a 33 ) 2 = (a 21 a 23 a 22 a 31 ) 2 g 1 g 2 + (a 31 a 32 a 21 a 33 ) 2 g 1 g 3 + (a 23 a 32 a 22 a 33 ) 2 g 2 g 3 Figure 4. The entangling subsystem is taken as a green slab. Rotating the green slab on the Euler angles (φ, θ, ψ) we get the pink slab that is oriented along the axes (x 1 , x 2 , x 3 ) associated with the HIC geometry and shown in Fig.5. x i are spatial coordinates in (2.1), a ij (φ, θ, ψ) are entries of the rotation matrix and Here φ is the angle between the ξ 1 axis and the node line (N), shown in Fig.4 in black, θ is the angle between the ξ 3 and x 3 axis, ψ is the angle between the node line N and the x 1 axis. We write the line element for the induced metric as and substitute the differentials dx M following from the embedding relations (3.3) in the RHS of (2.1) : We have 13 + g 2 a 2 23 + g 3 a 2 33 , g 12 (z, φ, θ, ψ) = g 1 a 11 a 12 + g 2 a 21 a 22 + g 3 a 13 a 32 , g 13 (z, φ, θ, ψ) = g 1 a 11 a 13 + g 2 a 21 a 23 + g 3 a 31 a 33 , g 23 (z, φ, θ, ψ) = g 1 a 12 a 13 + g 2 a 22 a 23 + g 3 a 32 a 33 , g 21 =ḡ 12 , g 32 =ḡ 23 ,ḡ 32 =ḡ 23 . (3.9) The determinant of the induced metric is and the Nambu-Goto action is where g, g 1 , g 2 , g 3 are functions of z andḡ 22 ,ḡ 33 ,ḡ 23 are functions of z and the Euler angles. Since there is no dependence of the integrand on ξ 2 , ξ 3 we can perform the integration over these variable that give the sizes of the parallelepiped in second and third directions: Here , L 1 , L 2 are the lengths of the parallelepiped in the first, second and third directions. can be fixed by boundary conditions, see below (3.17). In what follows we take L 1 = L 2 = 1. The action (3.12) is a particular case of the BI action This action defines the dynamical system with dynamical variable z = z(ξ) and time ξ. An effective potential is (3.14) This system has the first integral: From (3.15) we can find the "top" point z * (the closed position of the minimal surface to the horizon), where z (ξ) = 0: Finding z from (3.15) one gets representations for the length and the action S (3.13), that defines the HEE S (3.2) up to the factor 1/4: .
• The form of the effective potential V(z) does not depend on the slab orientation. Therefore, the location of the dynamical wall, defined by location of the minimum of the the effective potential V(z), is the same for all orientation for fixed z h and ν.

Geometric Renormalization
It is convenient to perform renormalization of the entanglement entropy by the subtraction the "disconnected" surface contribution from the "connected" one [41,45,53]. The contribution of the "disconnected" surfaces is given by the doubled area of the surface hanging from the boundary to the horizon or on the dynamical wall (if it exists at the considered parameters). The difference between the "connected" and "disconnected" parts S CD ≡ S conn − S diccon for arbitrary oriented slab is where z D is the minimum of the two values z DW and z h (z DW is the position of the dynamical wall, at this point V (z DW ) = 0).
Here we used the fact that the determinant of 3×3 induced matrix corresponding to the 3-dim body handing along the z-axis is where m 13 , m 23 and m 13 are given by (3.24). In (3.25) z = z D is the position of the horizon, or the dynamical wall. Here we consider the dynamical wall for HEE, that is defined as a position of the minimum of the effective potential V(z), It is interesting to note that V does not depend on the Euler angles. The HEE dynamical wall does not coincide with Wilson loop dynamical walls, which depend on the orientation of the loop in respect to the collision axes [76].
3.3 HEE for g 1 = 1, g 2 = g 3 = g and arbitrary orientation In this section we consider the metric (2.3). From the general formula (3.12) we get and a ij are given by (3.6) and satisfy the relation (a 21 (φ, θ, ψ)) 2 + (a 31 (φ, θ, ψ)) 2 + (a 11 (φ, θ, ψ)) 2 = 1. (3.29) So we can use a new parametrization a 11 (φ, θ, ψ) = cos ϕ, (3.30) and in this case Note that ϕ is nothing but the angle between the ξ 1 and x 1 axes. In this parametrization we have  This answer is equivalent to the HEE of the slab that is obtained by rotation of the initial slab oriented along the natural HIC axes about the ξ 3 axis by the angle ϕ, see Fig. 6.

The Holographic Entanglement Entropy for
We get the following representations for the character length and the action. For (2.4) we have determines the degree of divergence of S ϕ . Near z = 0 it has a different behavior for ϕ = 0 and ϕ = 0. However, we can use the universal renormalization (3.25).
( 3.43) The HEE is a UV divergent quantity, so the HEE needs the renormalization for z ∼ 0. Note that we can analyze the behavior of the integral at the upper limit z * . If [76]. Now we can determine the power of the integrand singularity at z = 0. It is defined by M behavior near z = 0: Here κ 0 is defined by the asymptotic of M and V at z = 0 and depends on the orientation, so we use subscripts to specify the different orientations. Taking into account the dilaton field φ(z, z h , c, ν) the asymptotic is given by (eq. (2.58) in [61]): where k(z h , ν, c) does not depend on z. Therefore we have the following asymptotic of the functions b s (z, ν) and M (z) at z → 0: The expression for B s (ν, c, z h ): The potential V asymptotic at z ∼ 0: The integrand (3.44) for the yXY case has a larger degree of divergence than for xY Y case, so Note that κ 0 > −1 for ν ≥ 1 (see Fig.7) and this fact means that the HEE is finite after the geometric renormalization for ν ≥ 1.

Minimal renormalization for ϕ = 0
Here we consider (3.38) in a particular case ϕ = 0, i.e. the smallest size of the slab is oriented along the longitudinal direction. In notation (3.13) we have: The UV divergencies are defined by behavior of M xY Y (z) at z ∼ 0 (3.48). We see that M xY Y (z) has an integrable singularity at z = 0 for ν ≥ 1.67.
For ν = 1 we have to perform a renormalization and for the renormalized HEE we get: where the last term means the indefinite integral over M as (z) at z = z * . For xY Y case and ν > 1.67 we have an integrable singularity. For xY Y case and 1 < ν < 1.67 we have to perform renormalization: where we use the expression for the indefinite integral of M :

Minimal renormalization for ϕ = 0
Here we consider a configuration for ϕ = 0. As we can see from (3.38), the type of UV asymptotic is the same for all ϕ = 0 and further we consider the case ϕ = π/2 for simplicity. From the invariance in the transversal directions we can take θ = 0, ψ = 0. This corresponds to parametrization (3.3) with φ = π/2, θ = 0, ψ = 0 and in this configuration the small subsystem orientation is delineated along the transversal direction y 1 -direction. We call this configuration the transversal and denote with the subscript yXY . In notation (3.13) we have: UV divergences are now defined by the asymptotic of M yXY (z) at z ∼ 0 (3.48). The plots of functions κ yXY (ν) and κ xY Y (ν) are presented in Fig.7. We see that M yXY (z) has a nonintegrable singularity for all ν values at z = 0.
For yXY case and ν > 1 there is a nonintegrable singularity and we have to perform a renormalization:

Entanglement Entropy Density and c-functions
It is also instructive to consider the entanglement entropy density that is defined as [86,87] (compare with [88,89]). The advantage of dealing with the HEE density is that it has no divergences. For fixed temperature T (i.e. fixed z h ) the entanglement entropy density can be obtained from equation (3.18). We see that it can be expressed in terms of the value of the effective potential V, The potential V(z) in the Born-Infield action (3.12) is the same for different orientations of the entangling domain. From (3.27) we get However, the angular dependence of the density exists due to the fact that the length given by (3.43), depends on the orientation. Definition (4.3) requires a comment. Despite the fact that the density of entropy does not contain divergences, its definition may depend on a finite renormalization. In particular, adopting the geometric renormalization (3.25) we have The anisotropic c-function can be defined as Here m is the scaling power. It depends on the model and we make few comments about m in the next subsection 5.3. Generally speaking, unlike conformal theories [98][99][100] and especially the holographic theories with conformal invariance [101][102][103][104][105][106], m ϕ is different in UV and IR, can depend on the orientation of the slab as well ( see considerations of candidates for c-functions in theories with Lorentz violation [90][91][92][93], and especially, [95][96][97]).
We can also consider

Numerical Results
In this section we display and discuss our main results using numerical calculations.
In what follows we set L = 1. We present all the plots according the color scheme of the lines indicated in Table 1.
various Magenta Brown Cyan various Gray Table 1. The color scheme of the lines used in the plots in this article

Entanglement Entropy near the Background Phase Transition
In this section we present plots of the entanglement entropy dependence on the geometric characteristics of the entangling region (orientation and thickness of the slab) and the thermodynamic characteristics of the medium (temperature and chemical potential) for the model (2.4). We find the slab HEE dependence on the smaller size numerically performing integration of (3.43) and (3.38), and then excluding the dependence on z * solving equation (3.43) for a given .
In this construction the location of the domain walls and the horizon plays a special role. In Fig.8 we show the appearance of the domain wall. The dynamical wall appears when the effective potential V gets a saddle point.
In Fig.8 we also illustrate the location of z DW for various choices of the horizon size z h for ν = 1 A) and ν = 4.5 B). Solid lines show V and dashed lines show V . We see that the location of the dynamical wall does not depend on the horizon. In Fig.8.C) we show the location of the dynamical wall for arbitrary 1 ≤ ν ≤ 4.5.
In Fig.9 we see that becomes larger when z * approaches z DW or the horizon. The location of the domain wall does not depend on the orientation and also does not depend on the horizon, see  We present in Fig.9.A) the dependence on z * for ν = 1, z h = 0.5, 1, 1.5 and µ = 0. The red dashed line shows the position of the dynamical wall. The plots in Fig.9.B) show the dependence on z * for ν = 4.5, T = 0.25 and µ = 0.2 for different angles ϕ = 0, π/6, π/4, π/2. The value of z * does not exceed z DW (T, µ), the position of the dynamical wall for given T and µ. Here z DW = 0.349.
The divergencies in (3.38) we tried as the minimum subtraction of the UV asymptotic in accordance with the formulas above in 3.4.1 and Sect. 3.4.2.

Entanglement entropy dependence on
In Fig.10 we present the dependence of the slab HEE on its smallest length at various values of the temperature. These temperatures are taken to be equidistant around the temperature of the BB phase transition at µ = 0.05 in the isotropic case. These temperatures are indicated on the (z h , T )-curve by the points (A, B, C) above the T BB and the points (D, E, F) below T BB for µ = 0.05 (see Fig.10.A)). For these temperatures and µ = 0.05 we depict the dependence of the HEE on . Here is taken less then c , the length at which the turning point of the connected entangling surface moves closer towards the domain wall z DW (see Fig.9). We see that the slab HEE undergoes a jump when the temperature crosses the phase transition line at the point T BB (µ). We also see that these jumps depend smoothly on at least for The similar picture takes place for the anisotropic background. The HEE dependencies on for the transversal and longitudinal orientations at equidistant values of the temperature around the BB phase transition for µ = 0.2 and ν = 4.5 are presented in Fig.11.A) and Fig.11.B), respectively. We see that decreasing the angle ϕ from ϕ = π/2 to ϕ = 0 we increase the HEE at fixed and fixed thermodynamic parameters (with our adopted method to eliminate divergences).  Therefore, we have demonstrated that the HEE depends on the regularization scheme.
• For the minimal regularisation scheme both the longitudinal and the transversal HEE -increase linearly at large and C yXY (T, µ) > C xY Y (T, µ). The linearity region for longitudinal entanglement entropy begins at lower compared to the transverse entanglement entropy.
they have different behavior for → 0 (compare with results of [35]).
• For the CD regularisation scheme both the longitudinal and the transversal HEE go to constants at large and a yXY (T, µ) > a xY Y (T, µ);

Entanglement entropy dependence on temperature
It is instructive to depict the HEE dependence on the temperature for a fixed value of near the phase transition. For this purpose we first find the dependence of z * on the length (at fixed horizon) and then calculate the HEE from (3.58) and (3.59).
The dependence of z * on z h for the fixed value of = 1 and different values of the angle ϕ at the chemical potentials µ = 0 are shown by blue lines in Fig.12 . For comparison, at this plot the dependence of z * on z h is shown for the isotropic case by green lines also for µ = 0. We notice a significant dependence on z h only on small values of z h (large BHs). Substituting the dependence of z * on z h to equations (3.58) and (3.59) we get the HEE dependence on the temperature at the fixed value of = 1 at varieties of the chemical potential µ for A) longitudinal and B) transversal cases. These dependences are shown by blue lines in Fig.13 A) and Fig.13 B) for various chemical potentials µ near the critical value µ cr = 0.349 for the longitudinal and transversal cases, respectively. For comparison, the dependences of the HEE on temperature ln z * Figure 12. The dependence of z * on z h at fixed = 1 and varieties of angles for the anisotropic case (blue lines for angles ϕ = 0, π/4, π/3, π/2 from bottom to top) and for isotropic case (green lines) at µ = 0).  for the same slab are shown by green lines for different chemical potentials µ near µ cr = 0.117 for the isotropic case, ν = 1 in both graphs. These plots Fig.13 show that in the vicinity of the phase transition temperature T = T BB (ν, µ) the HEE, like the thermal entropy (see Fig.2.B) undergoes significant jumps. From these plots we also see that for µ > µ cr (µ cr = 0.117 for ν = 1 and µ cr = 0.349 for ν = 4.5) the jumps disappear, similar to the thermal entropy behavior.
To summarise this subsection, we note that, • the HEE of the slab with constant thickness increase with the temperature increasing and fixed chemical potential; • the HEE jumps near the phase transition do not very dependent on the type of renormalization (geometric and minimal renormalization), see Fig.15.

Entanglement Entropy Density
The entanglement entropy density, defined in the previous Sect. 5.2 as (4.3), in our case takes the form In this section we study the dependence of the entanglement entropy density (5.5) on temperature and chemical potential near background phase transition. Also we compare the results in minimal and geometric renormalizations, (5.5) and (5.6).

Entanglement entropy density dependence on
In Fig.16 we show the dependence of the HEE density (5.5) on at fixed values of the temperature and chemical potential near the background phase transition line for the model (2.4). We see that the density function decreases monotonically depending on the length. Here, similarly to the considerations in the previous Sect.5.1, we consider the density of the entanglement entropy for different T and µ near the phase transition point (T BB = 0.2457, µ BB = 0.2), namely for the points A, B, C and D, E, F indicated in Fig.1.A. We see that the densities change significantly when we cross the phase transition point (T BB , µ BB ).
In Fig. 17 we show the dependence of the HEE density (5.6), on at the same temperatures and chemical potential as in Fig. 16. We get the similar behavior in both cases. For comparison we present the dependence of the densities (5.5) and (5.6) on in the isotropic case in Fig. 18.

Entanglement entropy density dependence on temperature
Now, we fix the length l = 1 and present in Fig.19 the slab HEE density dependence on the temperature for different angles, ϕ = 0 and ϕ = π/2, and different chemical potentials µ, below and above the critical values µ (anis) cr = 0.349 corresponding to ν = 4.5. For comparison, we present here the similar plots for isotropic cases where µ (iso) cr = 0.117. We see a dependence on orientation. In Fig.21 the angular dependence of the log of entanglement entropy density for angle values ϕ = 0, π/6, π/4, π/2 (thickness increases with increasing angle) at the chemical potential µ = 0 and µ = 0.2 are shown in A and B, respectively. We see that the shapes of the curves shown in Fig.21.A) are the same for different orientation angles, and the curves are just shifted with increasing angle ϕ. The same applies to the curves shown in Fig.21.B).
In Fig.22 we show the angular dependence of the entanglement entropy density for different values of (thickness increases with increasing a slab width) at the chemical potential µ = 0.2 and T = 0.25 in minimal A) and B) geometrical regularization.
Comparing the plots in Fig.21 and Fig.22 we see that the transversal HEE density weakly depends on the regularization type while there is a substantial difference for the longitudinal density in different regularization schemes.
In Fig.23 we show the angular dependence of the entanglement entropy density for different values of (thickness increases with increasing a slab width) at the chemical potential µ = 0.2 and T = 0.25 in minimal A) and B) geometrical regularization. Since the magnitude of the HEE density above the critical temperatures is  substantially larger than the values of HEE density below the phase transition, the plot of the jumps of HEE density Fig.24 is similar to the Fig.23. We see that the HEE density undergoes jumps near the BB phase transition and these jumps depend on the anisotropy, slab orientation and chemical potential. To summarise this subsection, note that similarly to the plots presented for the entaglement entropy in the previous subsection 5.1, in Fig.16 -Fig.23, we present dependence of the HEE densities (5.5) and (5.6) on the slab width l for equidistant values of temperature near background phase transition. We observe that the HEE density undergoes a jump when the temperature crosses the phase transition line. We also observe a non-substantially dependence on the regularization scheme.

c-functions
Now we consider the behavior of c-function numerically where m ϕ,F is a scaling power. It depends on the model, in our particular case on the anisotropy parameter ν, on the orientation (index ϕ) and the frame in which calculations are performed (index F ). In isotropic conformal invariant case in our dimension, m AdS 5 = 3. In our model in the EF we have This metric is invariant under the rescaling Using the suggestion of [95], see also [96,97], we get where n 1 = 1, d 1 = 1 and n 2 = 1/ν, d 2 = 2. Taking into account the behavior of the dilaton in UV [61] we get in the SF y j → = Λ K U V /2−1+1/ν y j , i = 1, j = 1, 2 (5.15) and n 1 = K U V /2, (5.16) In the SF we take In the particular case, ν = 4.5, we have where F is the index indicated the frame, F = EF, SF and We see that the UV scaling factor m EF U V,x,y approximates the function m EF x,y (z, z h , ν) relatively well up to holographic coordinates about a half of the horizon position,. The same concerns also the approximation m SF U V,y to m SF y (z, z h , ν), meanwhile the UV scaling factor m SF U V,x cannot be used as an approximation for m EF y (z, z h , ν). Moreover, all m F x,y (z, z h , ν) exhibit the singular behavior at large . Let us first consider the case of zero chemical potential. The dependence of temperature on the horizon for µ = 0 is presented in Fig.26. For large BH's the temperature increases when z h → 0, and for small BH's the temperature increases also for z h → ∞. At z h = z HP the HP phase transition takes place and all small BH's are unstable.

The c-function in the isotropic case
For the isotropic case in the EF, i.e. with b-factor without a dilaton, the c-function is defined as in the conformal invariant case In Fig.27 we depict the c ν=1,F vs iso by brown lines for z h = 0.8, 1, 1.5 (thickness increases with increasing decreasing z h ). In all these cases the disconnected parts hang up to the horizon (there is no dynamic wall in the EF) and in (5.27) we take For comparison, in Fig.27 the dependences of c ν=1,SF on , i.e. for calculations that are done in the SF with b S (z), are shown by green lines. In this case there is the dymamical wall and if it is not covered by the horizon, z DW = 1.413 < z h , the disconnected parts hang up to the dynamical wall. For the cases depicted on Fig.27 this corresponds to the cases z h = 0.8 and z h = 1. Note that z h = 1.3 corresponders to the thermodynamically unstable phase.
We observe a completely different behavior of green and brown curves. There are saddle points on the green lines, but not on the brown ones. These saddle points are related with existence of solution of the equation Notice, that by definition the DW corresponds to V z=z DW = 0, but → ∞ at z → z DW . In the second term the difference is zero at z = z DW , but → ∞ and we have uncertainty at z → z DW . Calculations show that c SF → 0 for at z → z DW .

The c-function in the anisotropic case
In anisotropic cases the definition of the c-function is modified [95][96][97] Here the power m ϕ depends on the orientation. For the particular cases of the transversal and longitudinal orientations these powers, m y and m x , are defined according to different scaling in transversal and longitudinal directions. We use here the scaling factors m ϕ defined according to the behavior of the metric in the UV region. Note that since the behavior of b s (z) and b(z) are different in the UV region, we get different factors m ϕ at different frames, and we put a subscript m ϕ,F to indicate this; F = EF for the Einstein frame and F = SF for the string frame (we have the same for the isotropic case). The factor m ϕ,F is defined according the formula (5.18) and (5.19) and for ν = 4.5 we have (5.20) and (5.21).
In Fig.28 and Fig.29 we present dependences of c on for ν = 4.5, various z h and µ = 0, 0.5 calculated in the EF and SF for transversal and longitudinal orientations. We see that we get rather different behavor in EF and SF .
• We see in Fig.28.A that in the transversal case the c-function has rather nontrivial behavior in the EF for µ = 0 already. Namely we see that there is a region of z h , here we see multivalued c-function as function of . This takes place not only at nonstable points (one of them is shown by orange), but also at z h,HP = 1.138, and this multivalued phenomena disappears at z h = 1.13.
We also see that c-function has local minimum and maximum at the points min and max correspondingly, and there is a region of min < < max where c increases when increases).
• We see at Fig.28.B that this multi-validity is preserved at µ = 0.5 and the region of z h where this multi-validity takes place, is more wide as compare with µ = 0 case.
• The c-function for the longitudinal orientation for µ = 0 is monotonically decreases, • While at µ = 0.5 and rather large z h (small BH), the c-function exhibits a nonmonotonic behavior around some cr (see the plot in the inset in Fig.29.D).
• In the SF the c-function for the transversal orientation decreases monotonically.
• For the longitudinal orientation there is a saddle point 0 (depending on z h and µ) and the c-function starts to decrease only for > 0 .

Origin of non-monotonic behavior of c-functions
Let us make a few comments about the origin of the behavior of the c-functions presented in Fig.28 and Fig.29. To understand the origin of such an untypical behavior we first consider the dependence of c-function on z * for the cases presented here, and then incorporate the dependence of on z * . η CD (ν) over the (z * ,ν) plane for the EF at different z h and µ, and different orientations (top line for the transversal case and bottom line for the longitudinal one). All these plots demonstrate that the c-function as a function of z * monotonically decreases while increasing z * . One can also see that the c-function essentially depends on the orientation and c-function is larger for the transversal case for the same thermodynamical parameters and the same z * . Also the c-function falls with temperature faster in the transversal case. Note that the c(z * ) does not change substantially when we change the chemical potential. This is due to the fact that the c-function depends mainly on the effective potential, which does not depend on the blackening function. There is a weak dependence of the form of the function c = c(z * ) on z h and µ, compare the graphs on the same lines in Fig.30. In the top line of Fig.32 contour plots of m F (ν) φ η CD (ν) over z * , ν plane for the SF at z h = 1 and ϕ = 0 (longitudinal case), and for different values of the chemical potential are presented. Here in : Fig.32.A) µ = 0 and in Fig.32.B) µ = 0.5. We see that two contour plots are practically identical. the bottom line of Fig.32 m F (ν) φ η CD (ν) as functions of z * for discrete values of ν are presented. For comparison of c-functions for µ = 0 (magenta lines) and µ = 0.5 (darker magenta lines) for discrete values of ν are presented on the same plot in Fig.32 .D). We see that two sets of lines are almost coincide. η CD (ν) vs z * (horizontal axis) and ν (vertical axis) for F=SF at z h = 1 and ϕ = 0 (longitudinal case), and for different values of the chemical potential: A) µ = 0 and B) µ = 0.5. We see that two contour plots are practically identical. Bottom line. C) The same as in A) for discrete values of ν; D) comparison of c-functions for µ = 0 (magenta lines) and µ = 0.5 (darker magenta lines) for discrete values of ν. We see that two sets of lines are almost coincide.

5.4.2
as a function of z * as a function of z * has a rather nontrivial form for some values of thermodynamic parameters, especially in the longitudinal case. In Fig.33.A) the functions = (z * , z h , µ) for F=EF for transversal case at z h = 1, 1.06, 1.138, 1.2 and µ = 0 are presented. It is interesting to note that already on stable backgrounds, i.e. z h,ip = 1.108 < z h < z HP = 1.138, the functions = (z * , z h , µ) are non-monotonic. Here z h,ip = 1.108 corresponds to the line, that has an inflection point (see the inset in Fig.33.A)). This means that the same can be realized at the different z * , or 3 different values of entanglement entropy, entropy density and c-function corresponds to the same length. Fig.33.B)- Fig.33.D) show what happens when we change the chemical potential. In Fig.33. B) µ = 0.2 and z h = 1.285 (an unstable point, the corresponding curves are displayed in orange), z h = 1.1266 (the point of the BB phase transition depicted in red), z h = 1.1 (the corresponding curve has an inflection point) and z h = 1. In Fig.33.C) µ = 0.2 and points with z h > 2.76 correspond to the stable backgrounds. An unstable point with z h > 2.6 is in orange. In Fig.33.D)   case) at z h = 0.3, 0.325, 0.35, 0.375 and 0.4. The right curve corresponds to surfaces touching a dynamic wall located at z DW = 0.354 and for the rest, z * is located near the horizons. At Fig.34. D) the same plots as in C) for ϕ = 0 (longitudinal case) are shown.
In Fig.31 contour plots for m F (ν) ϕ η CD (ν) over (z * , ν)-plane for the SF at z h = 1 and ϕ = π/2 (transversal case), and for different values of the chemical potential are shown: A) µ = 0.2 and B) µ = 0.5. We see that two contour plots are almost identical. In the bottom line of Fig.31. C) the plots similar to plots in A) for various discrete values of ν are presented; also the similar plots are presented in Fig.31. D) for various discrete values for ν. We see that there is a small quantitive difference in the right inserts of both plots: the coordinates of the saddle points for µ = 0.2 are z * | ν=1.5 = 0.7146, c| ν=1.5 = 0.4590 and z * | ν=1 = 0.9361, c| ν=1 = 0.1170 (A), and for µ = 0.5 are z * | ν=1.5 = 0.71748, c| ν=1.5 = 0.46135 and z * | ν=1 = 0.8489, c| ν=1 = 0.14450 (B). The multi-valued dependency on z * in holographic models was previously observed in [54]. The authors established a new type of the phase transition associated with the swallow-tail like structure for S HEE as the function of . For completeness, we present the dependence of the entanglement entropy in the EF on in Fig.36 and Fig.37. For the isotropic case the Van-der-Waals behavior of the entanglement entropy on a slab width takes place only for small black holes. For the anisotropic there is a small region of temperatures where S HEE has a multi-valued also for large black holes. c-function undergoes a jump at crit obtained from the self-intersection of (S HEE , ) diagram. lef t,L right,L in Fig.28.D) correspond to the turning points.

The c-function near the Background Phase Transition
In this subsection we present dependences of c on in a neighbourhood of the background phase transitions. First we consider c-functions defined by (5.30) in the SF. In Fig.38 we present this c-function for ν = 1 A) and ν = 4.5 for transversal B) and longitudinal C) orientations. In the main panels plots for the temperature below the phase transitions are shown, while the insets show graphs at temperatures above phase transitions. We see that the c-functions below the phase transition line are negligibly small as compared with c-functions above this line. Note that for isotropic case c-functions have saddle points, but there are no saddle points in the anisotropic case for transversal configuration. In Fig.39 c-fuctions versus in the EF are shown. For ν = 1 we calculate the c-function with η given by (5.5), Fig.39 A) and with (5.6), Fig.39. B. Here see curves corresponding to various values of the temperature below and above the phase transition point T = 0.3445, µ = 0.05: T = 0.33, 0.335, 0.34 (below) and T = 0.345, 0.35, 0.355 (above) (with increasing thickness with increasing temperature) and µ = 0.05. We also present the plots of c-functions for ν = 4.5 and transversal and longitudinal orientation and values of temperature below and above the phase transition, see Fig.39. C) and Fig.39. D). Here we see the prints of the length phase transition similar to [52]. This phase transition especially clearly is observed for transversal configurations at temperature above the phase transition line.  Table 2. A summary of the c-function behavior considered in the present work in the EF and the SF based on definition (5.7). "ISO" denotes isotropic solution. "L" and "T" denote the longitudinal and transversal orientations of the slap in the anisotropic metric with ν = 4.5. "PT" corresponds to configurations near the background phase transition.

Various c-functions as functions of
Let us summarize and comment on the results obtained in the previous Sect.5.3 and Sect.5.5. The summary of these results is presented also in Table 2.
• We have found, that c-functions in the EF decrease with increasing for not too large (see the left part of the Table 2).
-For the isotropic case in regions of (µ, T )-plane far away from the line of the phase transition the c-function decreases with increasing for all .
-For the anisotropic case the c-function decreases while increasing for all transversal and longitudinal orientations for small in all regions of (µ, T )-plane. * For the transversal orientation different behavior exhibits in the different regions. In the region of small µ in the transversal case the c-function has a local minimum at min after which it increases up to a local maximum at max and only then decreases. Increasing µ we get even more interesting behavior. The c-function as a function of becomes multivalued in some interval of , lef t,T < < right,T . For > right,T , crossing the line of the phase transition, the c-function has a jump. * For longitudinal orientation the c-function for small µ just decreases with increasing . But increasing µ we see that the c-function also exhibits multivalued behavior for lef t,L < < right,L and for > right,T it has jumps when crossing the phase transition line.
• The c-functions calculated in the SF also exhibit non-monotonic behavior in some cases, but not multivalued behavior (see the right part of Table 2).
-In isotropic case in the SF c-function increases while increasing in UV up to = max (T, µ). This non-monotonic behavior is related with dilaton behavior near z = 0.
-There is different behavior in anisotropic case in the SF. * The c-function exhibits monotonic behavior for the transversal orientation. It decreases with increasing for all . * The c-function increases with increasing in UV up to = max,SF (T, µ) and when decreases for the longitudinal orientation.
There are two reasons why we do not need to worry about all these.
• First of all, as has been mentioned in the text, a non-monotonicity in the anisotropic case is not in contradiction with any of the existing c-theorems as all of them are based on Lorentz invariance.
• The saddle points max,SF (T, µ) as well as the regions of multi-validity of the c-functions, l,L < < r,L or l,T < < r,T are located in the regions with large enough values of , where the definition of the c-function using the UV asymptotics of the solutions can be violated.

Entanglement Entropy Phase Transition
Let us remind that the criterion for the confinement/deconfinement phase transition in QCD is the behavior of the potential between quarks or the behavior of the temporal Wilson loops. The Wilson loops can also be computed in HQCD. It happens that location of the confinement/deconfinement line in the (µ, T ) plane can be close to the background phase transition [54,61,76], but not necessary fit it. For special models the phase transition of the HEE can be used as an indication of the HQCD phase transition [41,55]. The position of the background phase transition depends on the particular holographic model, see [110,111] and refs therein. It can be also located in the rightbottom part of the (µ, T ) plane starting from a point (µ cr , T cr ) and going down with increasing µ till the zero temperature. It also can be located in the left part of the plane, starting from T 0 at µ = 0 and going down with increasing µ till the point (µ cr , T cr ).
To find the location of the HEE phase transition on the (µ, T ) plane one has to find the location of points where the free energy F A corresponding to the reduced ρ A -matrix has a multi-valued behavior. It turns out that the effective free energy corresponding to the entangled region A dF A = −S A dT, (5.31) has a behavior similar to the behavior of the free thermal energy. We can define the density of the entanglement effective free energy as The density of the effective free energy as a function of temperature also has the swallow-tail behavior. This is due to the fact that, on the one hand, slabs that extend infinitely in two directions and being thick enough can take into account enough degrees of freedom to repeat the characteristic form of thermal entropy, and on the other hand, the threevalued behavior of the function z h = z h (T ) is an intrinsic feature of the background and inevitably leads to swallow-tail behavior for both ordinary and effective free energy even for thin slabs. But since the slab does not include into account all degrees of freedom, the thermal entropy, and the effective entanglement entropy do not exactly coincide, as well as the lines of the thermal and entanglement phase transitions are not coincide.
Let us summarize what we find studding behavior of the effective entanglement free energy in different scheme of regularizations and frames.
• We check the behavior of the effective entanglement free energy density as a function of temperature in the SF using different renormalization schemes, see Fig. 40.A. We see that there is no essential dependence on used regularization schemes in the SF. The transition points for the entanglement free energy in both regularization schemes almost coincide, also these points are very closed to the transition point obtained from the thermal free energy.
• We also compare behavior of the effective entanglement free energy and thermal free energy densities as the functions of temperature in the EF, see Fig. 40.B.
Here we see that the temperatures of the phase transition points almost coincide, while the values of the free energy densities at these points do not coincide.
• We also analise behavior of the effective entanglement free energy defined with the minimal renormalization scheme in EF and find the swallow-tail looks more flattened, see Fig.40.C. The same in the EF (the entanglement free energy density with CD in the EF is shown by gray solid line). C) The entanglement free energy density with η M in in the EF. The inset in C) shows that the form of T-dependence of the entanglement free energy density presented in the main panel is in fact of a swallow-tailed form. All these plots are done for the anisotropic case with ν = 4.5 and µ = 0.2.
• We have checked that the locations of the critical points extracted from the HEE density and the HEE itself almost coincide (numerical calculations have been done for different and orientations).
In Fig.41 phase diagrams for the thermal entropy (solid lines) and entanglement entropy densities (dashed lines) for various anisotropy parameters ν, ν = 1 (green lines), ν = 3 (khaki lines) and ν = 4.5 (blue lines) obtained by numerical calculations based on (3.43), (3.38) and (5.5) are presented. Generally speaking, the location of the phase transition for the HEE is rather close to the background phase transition, see Fig.41. Note that in contrast to Wilson loops behavior in anisotropic background [61], there is no visible orientation dependence of HEE phase transition line in anisotropic cases.

Conclusion
We have considered the most general anisotropic holographic model and found the expression for the HEE in terms of the Euler angles defined by the orientation of the slab-shape area in respect to HIC axes. In a particular case, we have considered the HEE for the model invariant in the transversal directions with the unique anisotropy scaling factor supported by the Einstein-Dilaton-two-Maxwell action [61]. The choice of the model [60] is motivated by agreement of the energy dependence of the produced entropy with the experimental data for the energy dependence of the total multiplicity of particles produced in HIC [62] in the anisotropic metric. This model describes multiplicity and quark confinement (for heavy quarks), predicts crossover transition line between confinement/deconfinement phases, anisotropy in hadron spectrum (for a short time after collisions) [76,111] and phase transition for the spatial Wilson loops [19].
We have calculated the HEE and its density in the holographic anisotropic model [61]. We have shown that the HEE and its density have significant fluctuations near the BB phase transition line in (T, µ)-plane for all values of the anisotropy parameter. The lines of thermal and entanglement entropies phase transition in (µ, T ) plane are different in the anisotropic cases but do not depend on the orientation of the entangling area. Note that for isotropic case ν = 1 the differences between these lines in the phase diagram plane are not visible. We have discussed an application of enormous increasing of the HEE as an indicator of the background phase transition.
We have studied the dependence of the c-function on the thickness of the entanglement slab. We have found saddle points of c-function as a function on as well as its multivalued behaviour. The obtained results are schematically presented in Table 2. The c-functions in the EF decrease while increasing for not too large . Moreover, in the isotropic case in regions of (µ, T )-plane remote from the line of the phase transition, the c-function decreases while increasing for all . In the anisotropic case, the c-function decrease with increasing for all transversal and longitudinal orientations for small in all regions of (µ, T )-plane. The c-functions calculated in the SF in some cases exhibit non-monotonic behavior, but there is no multivalued behavior here. In isotropic case in the SF c-function increases with increasing in UV up to = max (T, µ). This non-monotonic behavior is related with dilaton behavior near z = 0. There is different behavior in anisotropic case in the SF. As has been mentioned in the text, a non-monotonicity in an anisotropic case is not in contradiction with any of the existing c-theorems as all of them are based on Lorentz invariance. The saddle points as well as the regions of multi-validity of the c-functions are located in the regions where the definition of the c-function using the UV asymptotics of the solutions can be violated.
As to further development, we suppose to study modifications of the model [61] to include light quarks following [112], incorporate the chiral phase transition [113], and also perform the numerical calculations in full anisotropic case to incorporate the magnetic field, as has been done in [114].
We hope that the results presented in this paper, their interpretations and their further possible adjustment to the phenomenology data can be of interest for experiments at the future facilities of FAIR [115], NICA [116], for RHIC's BES II program and CERN, III run.