Chiral Separation effect in non-homogeneous systems

We discuss chiral separation effect in the systems with spatial non - homogeneity. It may be caused by non - uniform electric potential or by another reasons, which do not, however, break chiral symmetry of an effective low energy theory. Such low energy effective theory describes quasiparticles close to the Fermi surfaces. In the presence of constant external magnetic field the non - dissipative axial current appears. It appears that its response to chemical potential and magnetic field (the CSE conductivity) is universal. It is robust to smooth modifications of the system and is expressed through an integral over a surface in momentum space that surrounds all singularities of the Green function. In itself this expression represents an extension of the topological invariant protecting Fermi points to the case of inhomogeneous systems.

In the recent years the non -dissipative transport effects attract attention both in the framework of condensed matter physics and in the high energy physics [1][2][3][4][5][6][7][8]. An important arena for the experimental observation of these effects is given by Dirac and Weyl semimetals [9][10][11][12][13][14][15][16]. These materials also represent the bridge between high energy theory and condensed matter physics because the physics of fermionic quasiparticles inside them models physics of elementary particles. The chiral separation effect (CSE) is one of the non -dissipative transport effects. It has been proposed by M. Metlitski and A.Zhitnitsky [17]. The essence of this effect is appearance of a non -dissipative axial current in the direction of external magnetic field. The original calculations of this effect have been performed in the system of continuum Dirac fermions. Without external magnetic field this system is homogeneous. It has been found that in the chiral limit (when Dirac fermions are massless) the axial current is proportional to the external magnetic field strength F ij and the ordinary chemical potential µ counted from the Fermi point (the point in momentum space where fermion energy levels cross each other): The theoretical prediction of this effect has been followed by the number of papers discussing the possibility to observe it during heavy ion collisions [18][19][20][21]. The fireballs appeared during the heavy ion collisions are widely believed to contain the new state of matter called quark -gluon plasma [22][23][24][25][26][27][28][29][30][31][32]. In this state quarks are free and almost massless. Their masses may actually be neglected completely and we may speak of the system of true chiral fermions. External magnetic field appears here because the ion beams carry electric current. During the non -central collisions the two colliding ions, therefore, produce strong magnetic field orthogonal to their trajectories. The CSE results in the appearance of axial current within the fireball. After decay of the fireball the asymmetry in the distribution of outgoing particles carries the signature of the chiral separation effect. It is worth mentioning, that the CSE as well as its cousin -the chiral vortical effect (CVE) may also be relevant for the description of the quark matter under extreme conditions in the rotated neutron stars [33]. Actually, the CSE represents a certain incarnation of chiral anomaly. This relation has been discussed in a number of papers (see, for example, [14]). The similar conjecture has been proposed for the so -called Chiral Magnetic Effect (CME) [18,[34][35][36][37]. The important difference between the two is, however, that the true equilibrium theory does not admit the presence of the CME [5][6][7][8][40][41][42][43] while the CSE exists as an equilibrium phenomenon [44]. Out of equilibrium, however, the CME is back [39], which, possibly, manifests itself in experiments with Dirac semimetals [38].
We would like to notice several recent works on the CSE. In the framework of continuum quantum field theory the CSE was discussed, for example, in [3]. The lattice regularization has been used in [8] and in [44]. It has been shown that if the model is considered at small but finite temperatures, then the lattice regularization gives the conventional result for the CSE. It appears that in the framework of continuum theory the order of summation over Matsubara frequencies and integration over momenta is important. The uncertainty appears if integration over the 3 -momenta is performed first. This explains the importance of lattice regularization for the investigation of the CSE. [44] reports conventional expression for the conductivity of CSE in the lattice models, which describe chiral fermions at low energies. In the presence of the finite mass of the fermions the expression for the CSE is changed. On the formal level the theory with massless fermions suffers from various infrared divergencies [46][47][48]. For this reason the finite fermion masses are worth to be introduced even if the limit of small masses is assumed. Notice, that in [34] the neutral particles were discussed, and for them there are no infrared divergencies related to radiation of photons. Interaction corrections to CSE have been considered, for example, in [45]. It has been argued that the higher orders of perturbation theory do give corrections to the version of the CSE of massive fermions.
To the best of our knowledge the CSE has not been considered in sufficient details for the inhomogeneous systems. At the same time in any real situation, when the CSE takes place, it exists in essentially inhomogeneous systems. In case of the heavy ion collisions the chiral quarks in the fireballs exist in the presence of non -uniform environment.
The effective action of the chiral quarks here cannot have the form of a homogeneous Dirac action. Instead it should depend on various background fields that depend on coordinates. The external magnetic field in this problem as well as the quark chemical potential is also coordinate dependent. In case of the CSE in Dirac/Weyl semimetals the more or less uniform external magnetic field can be provided as well as the uniform chemical potential. However, the effective action for electron quasi -particles in itself is always not homogeneous even in the absence of external magnetic field. There are always impurities, and various sources of elastic deformations -both internal (dislocations and disclinations) and external (mechanical stress caused by external forces). In the present paper we concentrate on the latter situation -uniform external magnetic field and uniform chemical potential, but non -uniform fermionic action. In order to deal with the non -homogeneous systems we use the Weyl -Wigner formalism.
Weyl-Wigner formalism [49,50] as an alternative formulation of non-relativistic quantum mechanics has been developed by H. Groenewold [51] and J. Moyal [52]. Later it was adopted in some form both for the quantum field theory and for the condensed matter physics. It is based on the notions of the Weyl symbol of operator and the Wigner distribution function. The Wigner -Weyl formalism in quantum mechanics is often referred to as the phase space formulation. It is defined in phase space that is composed of both coordinates and momenta, while, the conventional formulation uses either coordinate space or momentum space representations. In the phase space formulation the quantum state is described by the Wigner distribution (instead of a wave function), while the product of operators is replaced by the Moyal product of functions defined in phase space.
Lattice field theories were proposed as a mathematical tool to deal with divergences in quantum field theories calculations in high energy physics. On the other hand, in addition to the "traditional", approach of quantum mechanics in solid state physics [63], there is a lot of activity of using quantum field theory ideas in condensed matter systems. The attempts to apply the numerical lattice QFT methods in condensed matter started with Monte Carlo simulations of graphene [64].
Although there were attempts to construct an exact phase space formalism for lattice theories and finite state quantum systems (Schwinger [53], Buot [54][55][56], Wootters [57], Leonhardt [58], Kasperkovitz [59] and Ligabo [60]), until recently such an approach has not been proposed. It has been developed recently in [61]. However, the present paper is based on the simplified version of lattice Wigner -Weyl calculus valid for the case when an inhomogeneity is sufficiently weak. This method is an approximation in case of condensed matter physics, where the lattice describes a real material, but is exact in case of high energy physics where the lattice is a mathematical tool. In the case of condensed matter systems, it is shown that this approximation holds for any physically reasonable fields from the experimental point of view. In this approach, Weyl-Wigner phase space formalism is used to calculate Dirac operators and Green's functions. These techniques are widely used in recent research [65][66][67][68][69] dealing with linear response to electromagnetic fields which are shown to be topological invariants, as quantum Hall conductance for example.
We consider the lattice tight -binding models of a rather general type. In these models the fermions are placed in four -component Dirac spinors. Extra internal indices of these spinors are also admitted (valley indices, real spin indices etc). As a result action for the fermionic quasiparticles contains the structure of 4 × 4 matrix to be expressed through Dirac matrices γ k for k = 1, 2, 3, 4, 5, and their derivatives σ kj = 1 4i [γ k , γ j ]. We are interested in the situation, when low energy effective theory of such models obeys chiral symmetry. For the homogeneous models this would mean that matrix γ 5 commutes or anti -commutes with the one -particle Hamiltonian in a small vicinity of Fermi surfaces/Fermi points, where the low energy effective theory arises. The Fermi surface manifests itself in the two -point Green functionĜ =Q −1 as the position of its singularities in momentum space. HereQ is the so -called lattice Dirac operator. For the non -homogeneous systems with sufficiently weak inhomogeneity operatorQ is not diagonal in momentum, and the notion of ordinary Fermi surface/Fermi point may be replaced by the coordinatedependent Fermi surface/Fermi point [62]. The space dependent Fermi point is known also as emergent gauge field. Mathematically the weakness of inhomogeneity means that the poles of the Wigner transformed two -point Green function G W (p, x) at each value of x are given by zeros of the Weyl symbol Q W (p, x) of operatorQ. (The precise definitions of Wigner transformation and Weyl symbol will be given in the next Section of the present paper.) For the case when inhomogeneity is more strong the zeros of Q W do not coincide with the poles of G W . An extension of the notion of Fermi surface to this case may be given by the hyper -surface in phase space (space of both coordinates x and momenta p). We may choose the hyper -surface, where Q W (p, x) = 0. The other possible definition is localization of the singularities of G W (p, x). One may also consider position of the singularities of a certain combination of Q W and G W entering expression for the CSE conductivity (to be specified below), and this is our way for the definition of hypersurface Ξ in phase space, which is the extension of the notion of the Fermi surface. For the non -homogeneous systems the low energy physics appears in a certain vicinity of Ξ, and we require that Q W and G W commute or anti -commute with γ 5 in this vicinity. Recall that the precise chiral symmetry cannot be maintained within the whole phase space for the lattice tight -binding models except for the marginal cases. An example of the marginal case is the system of naive lattice Dirac fermions. In this case the 16 fermion doublers appear, and the sum of their contributions to the CSE conductivity vanishes [44].
We show that under the above conditions (chiral symmetry of low energy effective theory) the axial current of CSE in the non -homogeneous system of general type is still proportional to external magnetic field. Being averaged over the whole volume of the system it may be expressed as where N is a topological invariant expressed through G W and Q W . (Here µ is counted from the level, whereJ 5 = 0.) This expression contains an integral over a surface Σ 3 in momentum space surrounding the singularities of an expression standing inside the integral (that is the hyper -surface Ξ mentioned above). Being the topological invariant N is robust to an arbitrary smooth modification of the given system as long as this modification does not break chiral symmetry of quasiparticles existing close to Ξ. All inhomogeneous lattice systems with chiral fermions (those with the chiral symmetry at low energies) may be subdivided into the homotopic classes. Within each class the operatorsQ are connected to each other by smooth deformation. The values of N are constant within each homotopic class. For the homogeneous representative of each homotopic class the value of N may be calculated easily. It is given by the number of the species of chiral Dirac fermions in the corresponding low energy effective theory. This allows to calculate easily the CSE conductivity N /(2π 2 ) for any inhomogeneous system.
In the present paper we do not consider interactions between the quasiparticles. It is worth mentioning that in the case of Integer Quantum Hall effect the similar problem has been considered recently (see [68]). It has been shown that the expression for Hall conductivity throughĜ has the same form as for the non -interacting case but withĜ replaced by the complete interacting Green function. Based on the approach of [68] we expect that the same refers to expression for N of Eq. (2). However, the consideration of this issue remains out of the scope of the present paper.

II. WEYL-WIGNER PHASE SPACE FORMALISM
In this section we briefly review the technique of Wigner transformation applied to quantum mechanics defined in infinite continuous coordinate space. Phase space formalism allows to describe quantum mechanics using c-functions instead of operators. Weyl-Wigner transformation of the matrix elements of operator (that is Weyl symbol of operator) represents such a correspondence.

A. Weyl symbol of operator and Wigner distribution function
We start from definition of an average of operatorÂ with respect to quantum state Ψ Here by x, y or p, q we denote the continuous coordinates or momentum respectively. For simplicity we consider the case of one -dimensional space R 1 . The generalization of our expressions to the case of D -dimensional space R D is straightforward. Let us change the coordinates: Then This gives Weyl symbol of operator A W (x, p) and Wigner distribution W (x, p) are defined as follows (where the momentum space representation may be obtained using the similar way) The Weyl symbol of the product of two operators, called the Moyal product, is defined as follows (this time in the momentum representation) changing variables where An example of the fermionic system is given by that of the Dirac fermions with the action whereψ and ψ are the Dirac spinor fields. Depending on the nature of the given problem they may be understood either as complex -valued spinors or as operators or as the Grassmann -valued fields. Dirac operator is defined here asD Action may be written as where we introduce shorthand notations ψ | and |ψ . Their meaning is ψ | x =ψ(x), and x |ψ = ψ(x).
C. Relation between the Green function and Dirac operator in Weyl-Wigner formalism In continuous case we have which can be rewritten in the "operator" form Applying Weyl-Wigner transformation, we obtain Eq. (20) is the Gronewold equation.

D. Wilson fermions
One of the methods to discretize the Dirac field is Wilson fermions model (see Appendix B) for more details. Action for Wilson fermions in Euclidean space discretized using rectangular lattice has the form where Here lattice sites are referred to as n, m. n +μ means the lattice site situated one lattice spacing ahead of n in direction µ. Indices α and β correspond to Dirac matrices γ k . Inserting Fourier transform of the field (where |M| = (2π) D is volume of momentum space) in the action (B1) we obtain the terms like Here e i is the unit lattice vector in the i -th direction. Then, the action in momentum space becomes where

E. Wilson fermions in the presence of gauge field
In continuous coordinates space, the transition from the Dirac operator in coordinates representation to momentum representation in presence of a gauge field A, is obvious.
In lattice theory it demands some work. In coordinates space, in the presence of gauge field, the Dirac operator takes the form while We restrict ourselves by the case of the U (1) gauge field A and then this parallel transporter is given by Using the Peierls substitution the partition function may be written in the momentum representation In fact, the same refers to the other lattice models defined by operatorsQ different from that of the model of Wilson fermions.

F. Approximate generalization of Weyl-Wigner fromalism from continuous space to lattice
Definitions of phase space formalism for the continuous case (7), (8), (9), (44) are modified for the case of discrete coordinates x n as follows, Where M is the first Brillouin zone and x n are the lattice points. We denote the trace of Weyl symbol as follows: In one -dimensional case we will have Let us consider the Moyal product of the Weyl symbols of two operators in the one -dimensional case. Weyl symbol of the product of two operators, called the Moyal product, is defined as follows (this time in momentum representation) changing variables we come to In case of the near diagonal operators, the only important region is around the origin, hence, we can change the integration area back to the square form and get the same result, for the Moyal product, as for the continuous space Although this expression has been derived for the case of one dimensional lattice, obviously it remains valid for the lattice models in any number of dimensions.

G. Weyl -Wigner transform -general properties
Weyl symbol of operators introduced above establishes correspondence between operators and functions defined on phase space. It satisfies a certain set of properties typical for the constructions of deformational quantization. In fact, we may use the other definitions of Weyl symbol, which posses the same properties in order to explore various nondissipative transport effects. Those basic properties are:

Star product identity
2. First trace identity 3. Second trace identity 4. Weyl symbol of identity operator

Star product
By tr we understand the trace of the operator itself in the original Hilbert space, while Tr is the trace operation defined for the Weyl symbols. In particular, for the approximate Wigner -Weyl calculus defined above it is given by Eq. (37).
Let us mention also the following useful property (See Appendix C.) It is worth mentioning, that the precise Wigner -Weyl calculus on the lattice has been proposed in [? ]. It obeys the above abstract properties precisely unlike the case of the introduced above approximate Wigner -Weyl calculus. This calculus has been designed for the models defined on rectangular lattices, and it allows to derive useful expressions for Hall conductivity to be discussed partially in the present paper as well. (Next we will use the similar constructions to investigate chiral separation effect.) It is important, however, that throughout the present paper we are limited to the case, when sums over lattice points in expressions for the total currents may be substituted by integrals. This is possible, when various inhomogeneities existing in the theory, are sufficiently small at the distances of the order of lattice spacings. Under these conditions the precise constructions of [? ] are not, in fact, necessary, and we are able to use the approximate Wigner -Weyl calculus for the lattice models described above.

A. Partition function variation
In this section we repeat briefly the considerations of [70] that lead to the construction of topological expression for the quantum Hall conductivity of non -homogeneous systems. We will use later this technique to consider the nonhomogeneous CSE.
In Euclidian space-time the partition function of a noninteracting fermionic system is expressed through the inverse bare Green function. We call it further for simplicity the Dirac operator and denote byQ. The partition function is given by Here ψ,ψ are the Grassmann -valued fields, while S is the action where we used Weyl symbols of operators Here by |ψ ψ| we denote operator with Grassmann -valued matrix elements ψ(x)ψ(y). For simplicity of notations we discretize both space coordinates and imaginary time. This is usual for the lattice discretized relativistic field theory and unusual for the lattice models of condensed matter physics. In the latter case we are able to take off the discretization of imaginary time at any step of calculations in order to arrive at the conventional expression Q = iω −Ĥ, whereĤ is one -particle Hamiltonian. Using Peierls substitution [44], in the presence of gauge field (54) takes the form Propagator of fermions is defined aŝ Its expression in momentum space is Variation of partition function may be expressed as follows We obtain that is here we use definitions from section II G. From now on we use continuum limit for the coordinates r n → x. This is possible if variations of fields on the distances of the order of lattice spacings are neglected.
In the presence of gauge field we substitute Variation with respect to the gauge field A → A + δA gives and Electric current is given by

B. Groenewold equation
Dirac operator and Green function obey the following equation Weyl-Wigner transform gives Groenewold equation As a result of the variation with respect to the gauge field A → A + δA hence

C. Topological invariance
Integrating (or summing on the lattice) local current density of (66), we obtain the total integrated current It is worth mentioning that J is not the conventional current I, which is defined as an integral of current density over the cross section of a given sample. Relation between the two may be understood easily for the homogeneous system of rectangular form with length L at finite temperature 1/β. Then J = βLI. The last equality in Eq. (71) is valid if spacial boundary conditions are periodic. For the variation of J we obtain Using identity and as well as and periodic boundary conditions in momentum space we get The cyclic properties of the trace will give Hence, J i is topological invariant in the presence of periodic spacial boundary conditions. Notice, that the above consideration fails in the presence of external electric field, when periodic boundary conditions cannot be imposed. Therefore, the appearance of non -vanishing response of J to external electric field does not contradict with the statement that J is topological invariant for the systems with periodic boundary conditions.

D. Linear response
Let us consider the case when the external gauge field C is present We assume here that Q W has an additional space dependence to that coming from the gauge field. The gauge field itself is divided to the background one B(x) and to that for which we are looking a linear response A(x).
Hence, the Dirac operator may be written as where The propagator may also be presented as a perturbation substituting this variation back into Groenewold equation (68), and leaving only the first order, we get using the identity (C5) we may write, up to the linear terms in we used above the following identity as well as Hence, we may represent δG W as follows

E. Electric current and gradient expansion
In case of sufficiently weak inhomogeneity discussed above (61) becomes or, generally speaking In the presence of gauge field As it was mentioned above the definition of current gives The total integrated current (defined as an integral over space -time of the current density) is given by response to external uniform electric field Using property (48) (it is valid if periodic spacial boundary conditions are imposed), we obtain Here and the first two terms in (100) are Since up to the linear terms in A we have The third term in (100) gives (we keep only the linear terms in A ij ): Writing current of Eq. (96) in terms of A k and A mn , we get We may define where and in case of periodic boundary conditions in momentum space j The j (2) i(mn) (x) is the local electric conductivity tensor since it is a coefficient in front of electromagnetic tensor. The average electric conductivity (we assume A ij = const) is to be obtained from where V (4) is the overall volume of Euclidean space -time while V is the three -dimensional volume. We defined Thus the term in electric current giving response to external field strength is

F. Hall conductance
Using expressions of Appendix A we obtain in the presence of external electric field Using (112) and (A17) we get We can define an averaged Hall conductivity For the 3 + 1 D systems it may be rewritten as follows [70]:

IV. CHIRAL SEPARATION EFFECT
A. Axial current Now we are equipped by all necessary tools to study chiral separation effect in non -homogeneous systems. This is an appearance of axial current in the fermionic systems with finite chemical potential and external magnetic field. The latter is supposed to be uniform, but the non -homogeneity of an arbitrary nature is present in the system even when the external magnetic field is off. It is assumed that the lattice model is similar somehow to the model of Wilson fermions -its "Dirac operator"Q is a 4 × 4 matrix expressed through the gamma matrices. However, the form ofQ not necessarily repeats that of the Wilson Dirac operator. Moreover, the considered operatorsQ are in general not diagonal with respect to momentum. The local axial current density may be defined as Repeating all steps of the previous section we come to the following term containing the linear response to external field strength: Integrating (or summing on the lattice) the local current given in (117), we get Here v is volume of the lattice cell. We have a useful formula v|M| = (2π) D . Dividing by the total 4 -volume we obtain the average axial current

B. (The absence of ) topological invariance for the total axial current
Integrating over all space -time (or summing over the lattice) the local current density of (117), we get the total integrated axial current As for the case of electric current, the total integrated axial current differs from the conventional total current I 5 through the cross -section of the sample. For example, for the case of a uniform rectangular sample in the case when nothing depends on time we have J 5 = βLI 5 , where β = 1/T is inverse temperature (assumed to be large) while L is the length of the sample. Variation of J 5 is given by Using identities and as well as we obtain (for the case of periodic boundary conditions) that under the trace we may substitute As a result we come to If γ 5 commutes (or anti -commutes) with G and Q then The latter condition means that the model possesses precise chiral symmetry. Under this (very restrictive) condition we obtain In the above considerations we also implied that there are no singularities of the Green function. The latter condition means that the fermions are gapped, and the Fermi energy is within the gap. We conclude that for the systems with gapped fermions in the presence of precise chiral symmetry the total integrated axial current J 5 i would be a topological invariant. In practise, however, the corresponding requirements are too restrictive. For example, for lattice Dirac fermions the presence of a gap (mass) excludes chiral symmetry. Therefore, in practise the total axial current cannot be topological invariant unlike the total electric current. (Recall that for the latter we need periodic boundary conditions, which exclude, in particular the presence of external electric field.) Below we will see, that the topological invariance ever appears in the consideration of axial currents of realistic systems in the form of robustness of the response of axial current to external magnetic field and chemical potential.

C. Axial current for gapless fermions at finite temperature
We are going to regularize the theory by finite (but small) temperature in order to deal with gapless fermions.
Matsubara frequencies are p 4 = ω n = One can see that ω n never equals to zero. Therefore, even for the system of massless/gapless fermions the propagator never has poles in momentum space. As a result the expression for the axial current is well -defined Introducing the chemical potential ω n → ω n − iµ we obtain Here |M 3 | is volume of the three -dimensional Brillouin zone. We represent the above expression as has the meaning of the CSE conductivity when external field strength corresponds to a constant magnetic field H: Here we assume antisymmetrization with respect to indices i and j. We will see below that for the wide range of systems −ǫ ijk σ ijk ′ = δ kk ′ σ CSE with a scalar CSE conductivity. We represent expression for the CSE conductivity as The limit of small temperature and CSE conductivity The limit of small temperature T → 0, N t → ∞, π Nt = ǫ → 0 allows to replace the sum by an integral. However, the point ω = 0 is excluded from this integral due to the above mentioned properties of finite temperature theory: Then (133) becomes we obtain In the static case both the Green function G and its inverse Q do not depend on time. Therefore all possible singularities of the above expressions are situated at ω = 0. Our integrals avoid these singularities due to the small but finite values of ǫ. In the absence of inhomogeneity (when the stars may be omitted in the above expressions) at ω = 0 the singularities of expressions standing in the integrals mark positions of Fermi surfaces. The presence of inhomogeneity changes positions of those singularities. However, weak inhomogeneity cannot force those singularities to approach boundary of the Brillouin zone. On the language of effective low energy continuum theory of our lattice model we say that the inhomogeneity cannot force singularities of the Green functions and their products to approach infinity.

E. CSE conductivity as a topological invariant
In Eq. (139) the integrals entering σ cancel each other except those in the small vicinities of the mentioned above singularities. That's why we may restrict integrations in Eq. (140) by the small regions of the Brillouin zone above/below the singularities. Here, in this region we assume the presence of precise chiral symmetry, which means that the effective low energy theory of our lattice model is chiral invariant (if the chiral anomaly is ignored). Recall, that the chiral symmetry cannot be maintained in the whole Brillouin zone of the majority of physical models. This is why the chiral anomaly appears. In the expression of Eq. (140), however, we restrict integrations to the region, where γ 5 commutes/anti -commutes with Q and G. We will see below that as a result the sum of the integrals in Eq. (139) represents a topological invariant, which does not depend on the form of the surface in 4D momentum space surrounding the singularities. We may deform this surface arbitrarily in such a way that it remains surrounding the singularities. This way instead of the two pieces of the infinitely close planes (situated above and below the singularities) we may integrate over the sphere in momentum space (this is illustrated by the Figure).
Thus we rewrite Here the integral is over Σ 3 , which is the 3D hypersurface in 4D momentum space. It consists of the two infinitely close pieces of the planes situated above and below the singularities of expression standing in the integral. Since γ 5 commutes/anti -commutes with G and Q in this region, we may rewrite this expression as and This expression is topological invariant provided that γ 5 commutes or anti -commutes with Q W and G W in the vicinity of Σ 3 . This may be proved easily using methods of Sect. IV B. Therefore, we may deform the surface Σ 3 in such a way that the deformed Σ 3 does not cross the singularities of expression standing inside the integral. As it has been mentioned above, we may deform the surface, for example, in such a way that it will have the form of a sphere (this is illustrated by the figure). When space inhomogeneities are sufficiently weak we are able to omit the star product in the above expression for the CSE conductivity (see Appendix D). This simplifies considerably calculation of σ CSE . Our expression is then reduced to that of [44]: Let us discuss for the definiteness example of the system with Wilson fermions in the presence of weakly varying external electric potential (see Appendix B). We assume that this model is used for the description of the continuous field theory with one massless fermion. This means that parameter m (0) is set to zero. In the absence of electric potential the model has one Fermi point at p = 0. Calculation of Eq. (144) in this case has been given in Appendix E. It gives N = 1. In the system of N fermions of this type we obtain N = N . The presence of electric potential, which is much smaller than the inverse lattice spacing, does not break chiral invariance of effective low energy theory. At the same time, it turns Fermi point into a Fermi surface. Nevertheless, as far as the introduced electric potential is weak enough, the value of N remains equal to that of the trivial homogeneous model. It is still given by Eq. (144), in which Σ 3 embraces the Fermi surface. Suppose that we modify functions g i and m entering Eq. (B4) in such a way that they become depending on a new parameter of the dimension of length l g . Suppose that pl g is of the order of unity around the singularities of expression under the integral of Eq. (143). In this case we already cannot omit star products in Eq. (143). However the CSE conductivity still remains unchanged if modification of functions g i and m is continuous. It is given by Eq. (143), in which Σ 3 embraces all singularities of an expression of the corresponding integral. These singularities already do not repeat the positions of Fermi surfaces but are supposed not to approach infinity (or the inverse lattice spacing).
We come to an interesting conclusion, that the CSE conductivity remains unchanged if we modify the fermionic system under consideration. There are only two requirements to this modification: 1) it is smooth; 2) it does not break chiral invariance of an effective low energy theory. This observation gives us the simple receipt how to calculate in practise the CSE conductivity for the given system. We should find the simple homogeneous system like that of the Wilson fermions (or that of overlap fermions [44]), which is connected to the given one by a continuous deformation. The Fermi surface is allowed to change its form during such a deformation. We should also require that the chiral symmetry is not broken in effective low energy theory neither in the original inhomogeneous system nor during the deformation to the mentioned simple homogeneous system. Finally, value of the CSE conductivity in the original complicated system is equal to its value in the simple one.

V. CONCLUSIONS AND DISCUSSIONS
In the present paper we discuss Chiral Separation Effect in essentially nonhomogeneous systems of chiral fermions. This is an appearance of non -dissipative axial current along the direction of external magnetic field. Our main result is expression for the CSE conductivity, i.e. for the response of axial current to external magnetic field and to chemical potential (both are assumed to be independent of time and space coordinates): Here G W is the Wigner transformed two point Green function, while Q W is Weyl symbol of lattice Dirac operator. Both are taken in the absence of external magnetic field. Both are matrices 4 × 4, which may model systems of Dirac fermions in the lattice regularized QFT or in the condensed matter systems with emergent Dirac/Weyl fermions. We assume that the inhomogeneities of the system under consideration are negligible at the distance of the order of lattice spacing. In this expression the integral is taken along the closed surface Σ 3 that surrounds positions Ξ(x) of all singularities of expression standing inside the integral. In the case of weak inhomogeneities (when the Moyal (star) products may be replaced by the ordinary ones) Ξ(x) at each x coincides with the position of the coordinate dependent Fermi surface. In a more general case the position of Ξ(x) cannot be predicted easily, and it does not, in general case, coincide with the position Ξ ′ (x) of the singularities of G W . One may consider surface in phase space Ξ (or surface Ξ ′ ) as an extension of the notion of Fermi surface to the case of nonhomogeneous system. In turn, Eq. (146) is an extension of the topological invariant N 3 responsible for the topological stability of the Fermi points/Fermi surfaces in relativistic quantum field theory [62] to the case of non-homogeneous systems. In both cases chiral symmetry remains essential for the topological stability of these objects. We accept assumption that along Ξ(x) at each x matrix γ 5 commutes or anti -commutes with Q W . This reflects the requirement that the considered fermions are chiral in the low energy limit -close to the zeros of Q W or poles of G W . Under this condition expression of Eq. (146) is a topoloigical invariant. It is robust to the smooth deformation of the system as long as Σ 3 does not cross Ξ(x) for any x. This property of N allows us to make a strong statement about the CSE conductivity of nonhomogeneous system. It remains non -sensitive to the particular form of lattice Dirac operatorQ. To be explicit, we may start consideration from the system with lattice Wilson fermions in the presence of weak external electric potential. In such a system the value of N can be easily calculated and is equal to the number of the species of Wilson fermions (see Appendix E). In fact, the same value of N (equal to the number of effective continuum Dirac fermions) can also be obtained for the other lattice regularizations (say, for overlap fermions) in the presence of slowly varying external fields. The details of calculations then repeat those of [44]. Next, we consider deformation of the system. It may result, for example, from elastic deformations in case of solid state systems or, say, from external gravitational field or from rotation in case of lattice regularized relativistic QFT. There may be the other reasons, which lead to modification of lattice Dirac operatorQ. We assume that under such deformations the chiral fermions remain chiral. This means that matrix γ 5 remains commuting/anti -commuting with Q W at low energies (i.e. in the vicinity of Ξ(x)). Surface Σ 3 should remain embracing Ξ(x) for all values of x. Besides, it remains in the region of the Brillouin zone, where γ 5 remains commuting/anti -commuting with Q W . The positions of Σ 3 , Ξ and Ξ ′ are assumed to remain in the small region of Brillouin zone. The size of this region is to be much smaller than the size of the Brillouin zone itself. In the effective low energy theory this region becomes the vicinity of zero momentum. On the language of effective low energy theory we require that the positions of singularities of the Green function do not approach infinity. This requirement is especially natural for the lattice regularized non -homogenous relativistic quantum field theory of chiral fermions. Recall that an ordinary Fermi surface of a homogeneous system cannot approach infinity. It is natural to suppose that Ξ being an extension of the notion of Fermi surface to the case of non -homogeneous systems, also cannot approach infinity.
Thus we come to an interesting conclusion. For any system of chiral fermions the axial current in the presence of constant external magnetic field is proportional to magnetic field. The coefficient of proportionality is equal to const + σ H µ, where µ is chemical potential. Coefficient σ H is universal. Irrespective of the particular form of the system it is equal to N It is worth mentioning that in the above consideration we ignored completely effect of interactions between the fermions. We expect, however, that the topological expression of Eq. (146) remains valid in the presence of interactions if we replace the noninteracting two -point Green functionĜ by the complete interacting Green function. (The same refers to its inverseQ.) This expectation is based on the recent consideration of the similar question for the quantum Hall effect of systems with interactions (see [68]). Notice, that radiative corrections to the CSE in QED calculated in [45] contain singularities in the limit of vanishing electron or photon mass. This reflects the well -known problem of infrared/collinear singularities in the systems with massless fermions interacting via an exchange by gauge bosons. For this reason in high energy physics the interacting fermions are typically considered with finite mass. In the case of Weyl and Dirac semimetals, however, the emergent Dirac fermions are true massless. Instead of the exchange by photons we are to consider Coulomb interactions. We do expect that the topological nature of CSE conductivity for chiral fermions survives in the presence of these interactions. Infrared divergencies are to be treated in this case carefully. However, the explicit consideration of this question remains out of the scope of the present paper.
We expect wide applications of results obtained here both in condensed matter theory and in relativistic high energy physics. In the latter case the chiral separation effect of chiral fermions is specific for the quark gluon plasma. The quark gluon plasma appears, in particular, during the heavy ion collisions. The fireballs existing just after a collision are in the presence of strong external magnetic field. Axial current of the CSE results after the decay of the fireball in an asymmetry of the outgoing particles. This asymmetry has been observed in experiment giving an evidence of the CSE. At the same time, all previous theoretical descriptions of this phenomenon dealt with the simplified homogeneous systems. This description is, of course, far from reality. In practise the system of chiral quarks inside the fireball is highly inhomogeneous. From the very beginning it was not clear which kind of CSE effect we may observe in this case. Our present answer is very simple -the complicated structure of the system does not affect at all the value of CSE conductivity σ H . It remains equal to its value from the homogeneous model. In order to prove this statement, we need to assume, though, that the external magnetic field is homogeneous. This requirement is not realistic as well. But overall, the pattern of the CSE that follows from the above study is much more close to reality than the one that follows from the consideration of naive homogeneous systems.
In condensed matter physics the analogue of the chiral separation effect emerges, first of all, in an effective description of fermionic superfluid 3 He − A. In this system the chiral fermions appear in vicinities of the Fermi points. The superfluid component of liquid may exist in the non -homogeneous state (say, forming various vortices) thus giving rise to the non -homogeneous effective theory of fermionic quasiparticles. We expect that indirectly the results obtained above may be extended to this effective theory. The direct extension is not possible here because atoms of 3 He are neutral, and therefore, only the emergent U (1) gauge field appears in this case. This emergent gauge field is axial rather than vector. And this fact complicates an analogy.
The direct observation of chiral separation effect may be performed for the Dirac and Weyl semimetals, where emergent relativistic chiral Dirac fermions appear in the vicinity of the Fermi points. Here we may apply constant external magnetic field and consider the axial current of the CSE that appears due to the finite chemical potential. Again, the realistic systems are not homogeneous because of the presence of impurities, because of dislocations and disclinations of the crystals, and because of elastic deformations. This is an ideal setup for the observation of the CSE of inhomogeneous systems, where, as we expect, our findings may become an important ingredient of its description. It is worth mentioning that there is a certain technical difficulty related to the experimental detection of the axial current in Dirac/Weyl semimetals. It is not as simple as the detection of electric current or even spin current. However, we hope that the future development of experimental technique will resolve this difficulty. Therefore, and In Euclidean space with Hence and Thus we come to the following relation between Euclidean field strength and real electric field of Minkowski space

Appendix B: Wilson fermions
Action for the Wilson fermions (defined on rectangular lattice) is Here indices m, n enumerate lattice points while α, β are spinor indices. By n +μ we denote shift of the lattice point by one link in the µ -th direction. Let us calculate the Fourier transform (1 − cos(p ν )) (B5) The two-point function The inverse matrix is defined by Then in four -dimensional space Inserting (B8) and (B3) into (B7), we get We are interested to find the inverse of Q in this case in the same way as in (B12). It can be shown that the contribution due to the commutation relations [p, A(x)] is the one lattice spacing shift to each harmonic. So, as long as the space dependence of the field is much larger than lattice spacing, we may neglect this contribution. Hence, the Green's function in the presence of gauge field under the conditions described above (we use (B13)) takes the form G αβ (p − A(i∂ p )) ≈ q −iγ q g q (p − A(i∂ p )) + m(p − A(i∂ p )) αβ k g 2 k (p − A(i∂ p )) + m 2 (p − A(i∂ p )) (B16) in the operator notations (and matrix form) we obtain G(p − A(i∂ p )) ≈ q −iγ q g q (p − A(i∂ p )) + m(p − A(i∂ p )) αβ k g 2 k (p − A(i∂ p )) + m 2 (p − A(i∂ p )) (B17) Appendix C: Translation operator and the Moyal product properties Let us consider translation operator e a∂x (f (x)g(x)) = f (x + a)g(x + a) = e a∂x f (x) (e a∂x g(x)) (C2) One can prove this identity as follows. The left hand side has the form The right hand side is Hence e a∂x (f (x)g(x)) = e a∂x f (x) (e a∂x g(x)) (C5) We may represent Moyal product of Weyl symbols of three operators as follows The poles are space dependent and for A 4 = iφ = 0 they correspond to the single points In the neighborhood of these points In the case of nonzero A 4 (x) = iφ(x) → 0 instead of the singularities concentrated at a point in momentum space for any given x we have singularities concentrated along the closed surfaces in momentum space. The form of these surfaces depends on x. More explicitly, we have spheres with the center at p = A(x) and radius |φ 0 (x)|. The Dirac operator becomes For the Green function we haveĜ (p − A(x)) = q −iγ q ξ q + 1 2 ξ 2 q k ξ 2 k + 1 4 j ξ 2 j 2 (E14) Topological invariant responsible for the CSE effect is given by Here surface Σ 3 surrounds all singularities of the Green functions at all values of x. The key point for the calculation of N is that we may deform the system smoothly removing the fields A, φ at all. This will bring us to a homogeneous system with A = φ = 0. Then we chose Σ 3 of the form of the 3 sphere that surrounds point p = 0. Then and At the same timeĜ (ξ) = − q iγ q ξ q + ξ 2 /2 ξ 2 + ξ 4 /4 ≈ − q iγ q ξ q ξ 2 (E18) and ∂ pi G = ∂ ξi G = −i γ i − 2ξ i (γξ)/ξ 2 ξ 2 (E19) as a result integral over x is irrelevant, and we have an integral over the surface of sphere (dσ i is a vector orthogonal to the surface of the sphere, its absolute value is equal to the area element): dσ i ξ 4 tr γ 5 γ a ξ a γ j γ k γ l = 4ǫ ijkl 48π 2 Σ3 dσ i ξ a ξ 4 ǫ ajkl = 24δ ia 48π 2 Σ3 dσ i ξ a ξ 4 = 1 (E20)