(Non)-Integrability of Geodesics in D-brane Backgrounds

Motivated by the search for new backgrounds with integrable string theories, we classify the D-brane geometries leading to integrable geodesics. Our analysis demonstrates that the Hamilton-Jacobi equation for massless geodesics can only separate in elliptic or spherical coordinates, and all known integrable backgrounds are covered by this separation. In particular, we identify the standard parameterization of AdS_p X S^q with elliptic coordinates on a flat base. We also find new geometries admitting separation of the Hamilton-Jacobi equation in the elliptic coordinates. Since separability of this equation is a necessary condition for integrability of strings, our analysis gives severe restrictions on the potential candidates for integrable string theories.


Introduction and summary
Over the last two decades, AdS/CFT correspondence [1] has led to great advances in our understanding of gauge theories and string theory on Ramond-Ramond backgrounds. A special role in this progress has been played by integrability, a surprising property of field theories, which allows one to compute spectrum, correlation functions, and scattering amplitudes [2] using an infinite set of conserved charges [3]. Originally integrable structures were discovered in the N = 4 super-Yang-Mills theory, but they have been extended to other systems 1 , and it is important to classify field theories admitting integrability. A promising approach to such classification, which is based on analyzing behavior of strings on a dual background, has led to ruling out integrability for the superconformal theory on a quiver [8] and for a certain deformation of N = 4 SYM [9] 2 . In this paper we will analyze integrability of strings on a large class of Ramond-Ramond backgrounds, rule out integrability for a wide range of field theories, and identify the potential candidates for integrable models.
To put our results in perspective, let us briefly review the status of integrability in N = 4 SYM (or in string theory on AdS 5 × S 5 ). In the planar limit, the field theory can be solved by the Bethe ansatz [11], and the spectrum of strings on the gravity side can be found by solving the Landau-Lifshitz model [12]. The agreement between these two exact solutions provides a highly nontrivial check of the AdS/CFT correspondence. The methods of [11,12] are applicable only to the light states, whose conformal dimension obeys the relation While the techniques of [11] are not applicable when inequality (1.1) is violated, the integrability might still persist in this case, at least for some sectors of the theory. Violation of (1.1) implies that excitations of AdS 5 × S 5 might contain D-branes in addition to the fundamental strings, and generic excitations of this type are very complicated. Fortunately, some states violating the condition (1.1) still have very simple behavior: these are the BPS states with 3 ∆ = J. On the gravity side of the correspondence, the BPS states are represented by supergravity modes or by D-branes, depending on the value of J. To have interesting dynamics, one can introduce some fundamental strings in addition to these BPS branes and to replace (1.1) by As we already mentioned, the planar techniques of [11] are not applicable to the states (1.2) which violate (1.1), and in this paper we will use alternative methods to study integrability of such states. The most useful version of the AdS/CFT duality involves field theory on R × S 3 , then the bulk configurations satisfying (1.2) are represented by fundamental strings in the presence of giant gravitons [13]. The interactions between these objects can be very complicated [14], but additional simplifications occur for semiclassical configurations of giant gravitons with J ∼ N 2 , which can be viewed as classical geometries [15]. In this regime of parameters, integrability of the sector (1.2) reduces to integrability of strings on the bubbling geometries constructed in [15]. If the CFT is formulated on R 3,1 , the counterparts of the bubbling geometries are given by brane configurations describing the Coulomb branch of N = 4 SYM [16]. In the latter case, one can introduce an additional deformation which connects AdS 5 × S 5 asymptotics to flat space and see whether integrability persists for such configurations.
It turns out that the answer to the last question is no, and this result discovered in [9] was the main motivation for our investigation. As demonstrated in [9], addition of one to the harmonic function describing a single stack of D3 branes destroys integrability of the closed strings on a new asymptotically-flat background. Since continuation to the flat asymptotics destroys the dual field theory, this procedure appears to be more drastic than a transition to the Coulomb branch, which corresponds to a normalizable excitations, so the latter might have a chance to remain integrable. In this paper we focus on geometries dual to the Coulomb branch of N = 4 SYM (either on R 3,1 or on R × S 3 ) and on similar geometries involving other D branes. A different class of theories, which involves putting D branes on singular manifolds, was explored in [8], where it was demonstrated that strings are not integrable on the conifold. From the point of view of field theory, this result pertains to the vacuum of N = 1 SYM with a quiver gauge group, which is complementary to our analysis of excited states in N = 4 SYM.
To identify the backgrounds leading to integrable string theories, one has to analyze the equations of motion for the sigma model and to determine whether they admit an infinite set of conserved quantities. Instead of solving this complicated problem, we will focus on necessary conditions for integrability and demonstrate that strings are not integrable on a large class of backgrounds created by D-branes. Integrability on a given background should persist for string of arbitrary size, and in the limit of point-like strings it leads to integrability of null geodesics 4 , which implies that the motion of a particle is characterized by 10 conserved quantities, matching the number of the degrees of freedom x i . Massless geodesics can be found by solving the Hamilton-Jacobi (HJ) equation, where S is the action of a particle, and g M N is the background metric. The system is called integrable if the HJ equation separates [17], i.e., if there exists a new set of coordinates Y M , such that S(Y 0 , . . . Y 9 ) = 9 I=0 S I (Y I ). (1.4) This also implies that the HJ equation has ten independent integrals of motion. Non-trivial examples of geometries leading to integrable geodesics include Kerr-Neumann black hole [18] and its generalizations to Kerr-NUT-AdS spacetimes in higher dimensions [19]. To rule out integrability of geodesics on a particular background, it is sufficient to demonstrate that separation (1.4) cannot be accomplished in any set of coordinates.
In this paper we will analyze the motion of massless particles in the geometries produced by stacks of parallel Dp branes and identify the distributions of branes which lead to integrable HJ equation (1.3). Specifically, we will focus on supersymmetric configurations of Dp-branes with flat worldvolume 5 , and assume that Ramond-Ramond (p + 1)-form sourced by the branes is the only nontrivial flux in the geometry. This implies that metric g ij has the form where the first term represents (p + 1)-dimensional Minkowski space parallel the branes, and H is a harmonic function on the (9 − p)-dimensional base space. We will further assume that the base space is flat. For a single stack of Dp-branes, the HJ equation separates in spherical coordinates, and this well-known case is reviewed in section 2. This section also includes another example, separation in elliptic coordinates, which plays an important role in the subsequent discussion.
Section 3 describes our main procedure, which is subsequently used to study geodesics on a variety of backgrounds. In subsection 3.1 we demonstrate that the HJ equation (1.3) does not separate unless the metric on the base space has the form ds 2 base = dr 2 1 + dr 2 2 + r 2 1 dΩ 2 and H depends only on r 1 and r 2 . The most general harmonic function H leading to integrable HJ equation is derived in section 3.2, and equations (3.50)-(3.52) summarize the main result of this paper for branes with flat worldvolume. The brane configurations giving rise to geometries (3.50)-(3.52) are analyzed in section 3.3. The results of section 3 imply that (Y 0 , . . . , Y 10 ) leading to separation (1.4) must reduce to the elliptic coordinates discussed in section 2. Section 4 discusses physical properties of the geometries leading to integrable geodesics. We demonstrate that separability persists for the wave equation beyond the eikonal approximation, a property that have been observed earlier for various black holes [18,20]. In section 4.2 we show that separability of the wave equation is associated with a hidden symmetry of the background, and we construct the conformal Killing tensor associated with this symmetry. In section 4.3 we apply the techniques of [8,9] to demonstrate that most backgrounds with separable wave equation do not lead to integrable string theories.
In section 5 our results are generalized to Dp-branes dissolved in D(p + 4)-branes, the system which plays an important role in understanding the physics of black holes [21]. We find that in asymptotically-flat space there are no integrable solutions apart from the sphericallysymmetric distribution of branes. However, there are several separable configurations in the near-horizon limit of D(p + 4)-branes, and the most general Dp-D(p + 4) configurations leading to separable HJ equation are presented in (5.17), (5.11).
In section 6 we consider another generalization by allowing the branes to rotate, i.e., by breaking the Poincare symmetry on the brane worldvolume. Although the general analysis of rotating branes is beyond the scope of this paper, we consider the special class of rotating solutions which cover all microscopic states of the D1-D5 black hole [22,23]. Such solutions are parameterized by curves in eight-dimensional space, and our analysis demonstrates that HJ equation is not separable unless this curve is a simple circle. For such configuration the separable coordinates have been found before [24], and we will demonstrate that these coordinates reduce to a special case of the general elliptic coordinates discussed in section 3.2.
In section 7 we consider a different class of rotating solutions, which describes all half-BPS states of IIB supergravity supported by the five-form field strength [15]. We demonstrate that there are only three bubbling solutions leading to separable HJ equation for the geodesics: AdS 5 ×S 5 , pp-wave and a geometry dual to a single M2 brane. In section 8 we discuss the equation for geodesics on bubbling geometries in M theory, and we demonstrate that the elliptic coordinates emerging from the separation of variables in the geometries of [15] coincide with standard parameterization of AdS 5 ×S 5 , AdS 7 ×S 4 , and AdS 4 ×S 7 .

Examples: spherical and elliptic coordinates
We begin with discussing two known examples of brane configurations which lead to integrable equations for the geodesics. First we recall the situation for a single stack of Dp branes. In this case the metric has the form and dΩ 2 8−p is the metric on a (8 − p)-dimensional sphere: 3) The Hamilton-Jacobi equation for a particle propagating in the geometry ( Solution (2.5) has (p + 1) integrals of motion p µ , 8 − p independent integrals coming from S L (this is ensured by the isometries of the sphere, the explicit form of S L is given in appendix A.1), and one integration constant coming from the differential equation (2.6). This implies that action (2.5) can be written in terms of 10 conserved quantities, so the geometry (2.1) leads to integrable geodesics. As demonstrated in [9], integrability does not persist for strings, unless a = 0 in (2.2). Spherical coordinates (2.1) will play an important role in our construction since any localized distribution of D branes leads to a harmonic function which approaches (2.2) at infinity. Thus, any set of separable coordinates must reduce to (2.1) far away from the branes. Our second example deals with two stacks of Dp branes separated by a distance 2d (see figure 1). The metric produced by this configuration is given by Introducing coordinates y i on the (7 − p)-dimensional sphere and separating variables in the action, This equation describes the motion of a particle in a two-center potential, and it is well-known that for p = 5 it can be separated further by introduction of the elliptic coordinates [25]: In appendix A.2 we review the separation procedure that leads to (2.10) and write down the explicit form of R (see (A.16)). For future reference we quote the asymptotic relation between elliptic and spherical coordinates: This completes our review of spherical and elliptic coordinates, and in the remaining part of this paper we will investigate whether the separation of variables persists for more general geometries produced by Dp-branes.

Geodesics in D-brane backgrounds
We now turn to the main topic of this paper: the analysis of geodesics in the geometry produced by D-branes: We will further assume that the base space is flat: The massless geodesics in the background (3.1) are integrable if and only if the HJ equation (1.3) separates in some coordinates (Y 0 , . . . , Y 9 ), as in (1.4). In section 3.1 we will use this separation and the Laplace equation to demonstrate that function H can only depend on two of the ten coordinates (Y 0 , . . . , Y 9 ). The further analysis presented in section 3.2 demonstrates that (Y 0 , . . . , Y 9 ) must reduce to a slight generalization of the elliptic coordinates presented in section 2. Although the Laplace equation (3.3) is satisfied away from the branes, any nontrivial function H must have sources, and the brane configurations leading to a separable HJ equation are analyzed in section 3.3.

Reduction to two dimensions
and H becomes a function of (r 1 , . . . , r k ). Moreover, since all rotational symmetries have been isolated, we conclude that 6 If the branes are localized in a compact region, then at sufficiently large values of and asymptotic behavior of function H is given by Here a is a parameter, which is equal to zero for the near-horizon geometries and which can be set to one for asymptotically-flat solutions.
Our goal is to classify the backgrounds (3.1) that lead to separable HJ equations for geodesics, and we begin with separating variables associated with symmetries. Poincare invariance of (3.1) and rotational invariance of the base metric (3.4) allow us to write the action appearing in the HJ equation (3.4) as Here p µ is the momentum of the particle in p + 1 directions longitudinal to the branes, and S is the part of the action that depends on coordinates y 1 , . . . , y d i of the sphere Ω d i . The label L i represents the angular momentum of a particle along this sphere, and it is defined by relation An explicit construction of S A special case of this equation with k = 1 separates in spherical coordinates (see section 2), and we will now prove the equation (3.10) does not separate if k > 2. The separation of (3.10) for k = 2 will be discussed in section 3.2.
Separability of equation (3.10) should persist for all values of angular momenta, so we begin with setting all L j to zero and p 1 = . . . = p p = 0. The resulting equation (3.10) can be viewed as a HJ equation on an effective (k + 1)-dimensional space A general theory of separable HJ equations on curved backgrounds has a long history (see [26]), and a complete classification is presented in [27,28]. In particular, this theory distinguishes between ignorable directions (which correspond to Killing vectors) and non-ignorable ones. Clearly, the time direction in (3.11) is ignorable, but (3.11) does not have additional Killing vectors which commute with ∂ t . Indeed, any such vector would be a Killing vector of the kdimensional flat space, so it must be a combination of translations and rotations in r i . However, the asymptotic behavior of function H (3.7) breaks translational symmetry, and our assumption (3.5) destroys the rotational Killing vectors, so if the HJ equation for the metric (3.11) separates in some coordinates (t, x 1 , . . . x k ), only one of them (specifically, t) can correspond to an ignorable direction. Moreover, the discrete symmetry t → −t of (3.11) guarantees that t does not mix with (x 1 , . . . x k ) in the metric, and such orthogonality leads to simplifications in the general analysis of [27,28]. Specifically, according to theorem 6 of [28] separation of variables in (3.11) implies that 7 (1) There exist k independent conformal Killing tensors A (a) with components (A (2) Each of the one-forms dx l = ω l i dr i is a simultaneous eigenform of all A ij (a) with eigenvalues ρ (a,l) . This implies that a projector P (l)i j onto dx l satisfies equation (3) The metric in coordinates dx l is diagonal, so projectors P (l)i j commute with A ij (a) and g ij , and projectors with different values of l project onto orthogonal subspaces. Since the number of projectors is equal to the number of coordinates, we arrive at the decomposition (4) The components of A (a) along the Killing direction satisfy an overdefined system of differential equations: Notice that in [28] the theorem is formulated in terms of coordinates x i , so it does not use the projectors. For our purposes the covariant formulation given above is more convenient, in particular, to rule our the separation of variables, we will have to work in the original coordinates r i and demonstrate that the required Killing tensors A (a) do not exist. Notice that equation (3.14) can be rewritten in the form which does not refer to projectors (and thus to coordinates x i ): multiplying this equation by g ij and using (3.13), we find This relation is equivalent to (tti) component of the equation for the Killing tensor A (a) .
The theorem quoted above implies that A (a) ij are conformal Killing tensors on the k-dimensional base of (3.11), and, for the flat base, all such tensors can be written as quadratic combinations of k(k + 3)(k + 4)(k + 5)/12 Killing vectors [29] 8 : Equations (3.12) and (3.15) give severe restrictions on coefficients b (a) m,n , but fortunately the consequences of (3.12) have been analyzed elsewhere. Indeed, equation (3.12) does not involve g tt , so it remains the same for H = 1, when (3.11) gives the flat space, and the corresponding HJ equation gives an eikonal approximation for the standard wave equation. It is well-known that in 3 + 1 dimensions (k = 3) the latter can only be separated in ellipsoidal coordinates and their special cases [26], and generalization of this result to k > 3 is presented in [27,28]. This leads to the conclusion that the HJ equation in the metric (3.11) with k > 2 can only separate in ellipsoidal coordinates or in the degenerate form thereof. Before ruling out this possibility, we briefly comment on the peculiarities of the two-dimensional base. In this case the conformal group becomes infinite-dimensional, so the base space admits an infinite number of the conformal Killing tensors. This situation will be analyzed in section 3.2.
To summarize, we concluded that for k > 2, the HJ equation can only be separable in some special case of ellipsoidal coordinates (x 1 , . . . , x k ), which are defined by [30,31] Here (a 1 , . . . , a k ) is the set of positive constants, which specify the ranges of variables x i : Rewriting the metric (3.11) in terms of x i and substituting the result into (3.10), we find the HJ equation in ellipsoidal coordinates (see appendix B for detail): Here h i is defined by Function H appearing in (3.19) must satisfy equation (3.6) away from the sources, and appendix B we demonstrate that for such functions equation (3.19) never separates in ellipsoidal coordinates. This shows that the HJ equation can only be integrable for k = 1 (the situation considered in section 2) and for k = 2, which will be analyzed in the next subsection.

Separation of variables and elliptic coordinates
In the last subsection we have demonstrated that the HJ equation (3.10) is not integrable unless the flat base has the form (3.4) with k ≤ 2 and H is a function of r 1 and r 2 only. In section 2 we have already discussed k = 1 and this subsection is dedicated to the analysis of k = 2. To simplify some formulas, we slightly deviate from the earlier notation and write the metric (3.1) with the base (3.4) for k = 2 as The connection to coordinates of (3.4) is obvious: with d 1 = m, d 2 = n. The arguments presented in the last subsection ensure that H appearing in (3.21) can only depend on r and θ, i.e., the distribution of Dp branes that sources this harmonic function is invariant under SO(m + 1) × SO(n + 1) rotations. Notice that The Poincare and SO(m + 1) × SO(n + 1) symmetries of (3.21) lead to a partial separation of the Hamilton-Jacobi equation (1.3) for geodesics (this equation is a counterpart of (3.8)): where S (m) We recall that H(r, θ) is a harmonic function describing the distribution of D branes, so away from the sources it satisfies the Laplace equation on the base of the ten-dimensional metric (3.21): 1 r m+n+1 ∂ r (r m+n+1 ∂ r H) + 1 r 2 sin n θ cos m θ ∂ θ (sin n θ cos m θ∂ θ H) = 0. (3.27) Let us assume that the massless HJ equation (3.26) separates in coordinates (x 1 , x 2 ). In particular, this implies that the metric in (x 1 , x 2 ) must have a form [32,33] where g 1 (x 1 , x 2 ) and g 2 (x 1 , x 2 ) satisfy the Stäckel conditions: We wrote the Stäckel conditions (3.29) for an arbitrary number of coordinates to relate to our discussion in section 3.1, but in the present case (k = 2) equation (3.29) gives only two relations (i = 1, j = 2, l = 1, 2): In particular, we find that so by adjusting function A in (3.28) we can set We can further redefine variables, x 1 →x 1 (x 1 ), x 2 →x 2 (x 2 ) to set g 1 = g 2 = 0, at least locally. 9 Introducing x =x 1 , y =x 2 , we rewrite (3.28) as dr 2 + r 2 dθ 2 = A(x, y) dx 2 + dy 2 . (3.33) To summarize, we have demonstrated that in two dimensions the Stäckel conditions (3.29) imply that coordinates x i separating the HJ equation are essentially the same as the conformally-Cartesian coordinates (x, y) in (3.33) 10 . In higher dimensions, conditions (3.29) are less stringent than the requirement for coordinates x i to be conformally-Cartesian: for example, conditions (3.29) are satisfied by spherical coordinates that have g 1 = 0, g 2 = 2 ln x 1 , g 2 = 2 ln(x 1 sin x 2 ), (3.34) but there is no change of coordinates of the formx i (x i ) that allows one to write implies that A must be equal to constant, so if conformally-Cartesian coordinatesx i separate the HJ equation, thenx i must be Cartesian 11 . Returning to k = 2, we will now find restrictions on A(x, y) and H(r, θ). First we definẽ R(x, y) byR (x, y) ≡ R(r, θ). (3.37) Then equation (3.33) can be used to rewrite (3.26) in terms ofR, and we find the necessary and sufficient conditions for the Hamilton-Jacobi equation (3.26) to be separable: Here M is an effective mass in (p + 1) dimensions defined by The construction of the most general harmonic function H(r, θ) that admits the separation of variables (3.38)-(3.39) will be performed in three steps: 1. Determine the restrictions on function A(x, y) imposed by equation (3.38).
2. Use equation (3.39) to find H(r, θ) corresponding to a given A(x, y). To implement the first step, it is convenient to introduce complex variables Here l is a free parameter which has dimension of length. Rewriting equation (3.38) in terms of complex coordinates, we conclude that 12 ∂z ∂w ∂z ∂w = 0, (3.44) so z(w,w) is either holomorphic or anti-holomorphic. Without loss of generality we assume that then equation (3.43) gives an expression for A(x, y): The second step amounts to rewriting equation (3.39) as Implementation of the third step amounts to finding expressions for h(w) and U 1 (x), U 2 (y) which are consistent with Laplace equation (3.27) for function H. Physically interesting configurations correspond to branes distributed in a compact spacial region, so at large values of r function H behaves as Here Q is the total brane charge, and a is a parameter, which can be set to zero for asymptotically flat space, and which is equal to zero for the near-horizon geometry of branes (cf. equation (2.2)). Keeping only the two leading terms in (3.48), we conclude that the Hamilton-Jacobi equation (3.26) separates in coordinates (r, θ), and the solution is similar to the first example discussed in section 2. The subleading corrections in (3.48) obstruct the separation in spherical coordinates, but for some harmonic functions H the Hamilton-Jacobi equation (3.26) separates in a different coordinate system, as in the second example discussed in section 2. Notice that at large values of r the new coordinate system must approach spherical coordinates since (3.48) depends only on r, and for elliptic coordinates such asymptotic reduction was given by equation (2.13).
In appendix C we construct the most general expressions for h(w), U 1 (x), U 2 (y) by starting with asymptotic relations (3.48) and 49) and writing expansions in powers of l/r. Notice that these asymptotics lead to the unique value of l for a given configuration of branes (for example, l = d/2 in (3.55)). Requiring that function (3.47) satisfies the Laplace equation (3.27), the boundary condition (3.48), and remains regular at sufficiently large r, we find three possible expressions for h and U 1 , U 2 : 13 The derivation of these constraints is presented in appendix C. For six-and seven-branes 14 (i.e., for m + n < 2), harmonic function is slightly more general: III : The terms proportional to P are ruled out by the boundary condition (3.48) for p < 6, but they are allowed for p = 6, 7.
Notice that case III can be obtained from case II by interchanging the spheres S m and S n , so, without loss of generality, we can focus on solutions I and II. Case I corresponds to spherical coordinates, and case II corresponds to elliptic coordinates discussed in the second example of section 2: as demonstrated in appendix C, expression (3.51) for h(w), is equivalent to Thus (x, y) are equivalent to the elliptic coordinates (ξ, η) defined by (2.12). For future reference, we rewrite the metric (3.21) in terms of x and y (see (C.12)):

Properties of the brane sources
In section 3.2 we have classified the geometries which lead to separable Hamilton-Jacobi equations for geodesics. The Laplace equation (3.27) played an important role in our construction, but this equation is only satisfied away from the sources. In this subsection we will analyze the solutions (3.50)-(3.52) to find the sources of the Poisson equation and to identify the corresponding distribution of branes. As already discussed in section 2, the spherically symmetric distribution (3.50) corresponds to a single stack of Dp branes. Since solution (3.52) can be obtained from (3.51) by interchanging S m and S n , it is sufficient to discuss only (3.51). Function H defined by (3.51) becomes singular at for all values of (m, n), at 15 and at x = 0 if n > 1. For m = 0, we have an additional locus (3.58), and repeating the steps above, we find a counterpart of (3.61): (3.62) Equations (3.61) and (3.62) describe two point-like sources, and we have already encountered these points in the original elliptic coordinates discussed in section 2. In the remaining part of this section we will focus on m > 0. For n = 0, 1, the m-dimensional sphere described by (3.61) is the only singularity of the harmonic function, and for n > 1 there is an additional locus given by (3.59). To formulate (3.59) in terms of r and θ, we first use the expression (3.51) for h(w), to rewrite (3.59) as The last relation can be rewritten as e w + e 2w − 4 = 2 ⇒ re iθ + r 2 e 2iθ − 4l 2 = 2l, (3.64) so there must exist an angle ψ, such that Solving the last equation for r, we find The right-hand side of the last equation must be real, this implies that (3.59) is equivalent to n > 1 : Substituting this relation into (3.55) and recalling that d = 2l, we conclude that ψ = y.
To give a geometric interpretation of (3.67), we recall that, for m > 0, y ranges from zero to π/2, so (3.67) describes a line connecting r = 0 with r = 2l. The geometry (3.21) has an m-dimensional sphere attached to every of this line, so the singular locus (3.67) has a topology of (m + 1)-dimensional disk. For m = 0, y, and θ range from zero to π, so (3.67) represents a line connecting two singular points (3.61), (3.62). The pictorial representation of singular loci in (r, θ) plane is given in figure 2.
We will now combine (3.61) and (3.67) to analyze the brane distribution for m > 0.
In this case, the m-dimensional sphere described by (3.61) is the only singularity of the harmonic function, and in the vicinity of this singularity we find To give a geometric interpretation of this expression, we consider the base metric (3.21) in the vicinity of singularity (3.61): It is convenient to introduce polar coordinates (R, Φ) by Recalling that θ varies from − π 2 to π 2 when n = 0, we conclude that Φ ∈ [0, 2π) (or Φ ∈ [−π, π), see below), as long as R remains small. Then metric (3.69) takes the standard form To write (x, y) in terms of (R, Φ), we use the real form (3.55) of the holomorphic map (3.51). In particular, for small x and y we find Extraction of the square root from the last expression should be done carefully: positivity of the harmonic function (3.68) (or (3.51)) requires that x > 0. Thus we can write as long as Φ ∈ [−π, π). The range Φ ∈ [0, 2π) is equivalent from the point of view of (3.70), but it leads to a more complicated counterpart of (3.73), and it will not be explored further. Substitution of (3.72) and (3.73) into equation (3.68) leads to a simple expression for the harmonic function in the vicinity of the sources: We recall that R = 0 corresponds to the m-dimensional sphere (3.61) with radius 2l. Notice that this expression never becomes negative since Φ ∈ [−π, π).
Here again, the m-dimensional sphere described by (3.61) is the only singularity of the harmonic function, and in the vicinity of this singularity we find As before, we used (3.72) to rewrite the harmonic function in terms of coordinates (3.70), and in this case Formula (3.75) gives a standard harmonic function in three dimensions transverse to S m , so this configuration corresponds to D-branes uniformly distributed over S m .
In this case the singularity consists of a line (3.67) connecting r = 0 and (3.61) with sphere S m fibered over it. In the vicinity of θ = 0, the base of the metric (3.21) becomes and the locus (3.67) is an m + 1-dimensional ball with metric The singularity corresponds to x = 0 (recall (3.57) and (3.59)), and the leading contribution to the harmonic function (3.51) for small x is For small x, metric (3.56) becomes Away from y = 0, we can neglect x 2 in comparison with sin 2 y, then function (3.78) describes the Coulomb potential produced by D-branes uniformly distributed over S m . Rewriting (3.79) as we find the charge density ρ: As expected, the charge density vanishes on the boundary (3.61) of the ball, where y = 0.
To summarize, we found that for n ≤ 1 the brane sources are localizes on the m-dimensional sphere, they produce a Coulomb potential (3.75) for n = 1 and a potential (3.74) with a fractional power of the radial coordinate R for n = 0. For n > 1, the branes are located on a line connecting r = 0 and (3.61) with sphere S m fibered over it ((R, Ω m ) subspace of (3.80)). These sources produce a Coulomb potential in (n + 1) transverse directions with charge density (3.81).

Beyond geodesics: wave equation, Killing tensors, and strings
The main goal of this paper is identification of backgrounds which can potentially lead to integrable string theories. As discussed in the introduction, integrability can be ruled out by looking at relatively simple equations for the light modes of strings (massless particles), and in the last section we demonstrated that classical equations of motion for such particles are integrable only for the harmonic functions given by (3.50)-(3.52). It is natural to ask whether separability of the HJ equation (1.3) in the elliptic coordinates (3.55) persists at the quantum level and whether it is related to some hidden symmetry of the system. In section 4.1 we analyze integrability of the wave equation, a quantum counterpart of (1.3), and in section 4.2 we identify the Killing tensor responsible for the separation. Finally in section 4.3 we investigate the question whether integrability of geodesics persists for finite size strings.

Separability of the wave equation
In this subsection we will analyze the wave equation which governs dynamics of a minimallycoupled massless scalar: Most Dp-branes generate a nontrivial dilaton, so the last equation would look differently in the string and in the Einstein frames 16 , and here we will focus on the most interesting case of the Einstein frame: Since the HJ equation (1.3) arises in the eikonal approximation of (4.1), the arguments presented in sections 3 imply that (4.1) is not integrable unless the metric ds 2 base has the form (3.56) and H is given by (3.51) 17 , although these conditions are not sufficient for integrability of (4.1) 18 .
Let us now demonstrate that the wave equation (4.1) does not separate in the string metric unless p = 3. The string-frame counterpart of (4.6) can be obtained by formally setting e 2Φ = 1 in that equation: Clearly, the multiplicative separation (4.7) does not work for this equation unless p = 3. Since coordinates (x, y) are uniquely fixed by the discussion of the HJ equation (which is an eikonal limit of (4.11)) presented in section 3, to rule out the separation, it is sufficient to show that a substitution of into (4.11) does not lead to separate equations for F and G for any fixed function P (x, y). To demonstrate this, we perform such substitution and rewrite the result as Since H and P are fixed functions, the first two lines of (4.13) remain the same for all F and G, so equation (4.13) does not separate unless its third line is only a function of x and the forth line is only a function of y. This implies that (4.14) Functions P 1 and P 2 can be absorbed into F and G (recall (4.12)), so we set P 1 = P 2 = 1. Direct substitution into (4.13) shows that the third line of that equation obstructs separation unless p = 3.
To summarize, we have demonstrated that while the HJ equation is integrable in the elliptic coordinates (3.55), its quantum version may or may not be separable depending on the frame. In particular, the equation for the minimally-coupled massless scalar separates in the Einstein, but not in the string frame.

Killing tensor
Separation the Hamilton-Jacobi and Klein-Gordon equations implies an existence of nontrivial conserved charges which are associated with symmetries of the background. In the simplest case, such symmetries are encoded in the Killing vectors, which correspond to invariance of the metric under reparametrization where V M (x) satisfies the equation This symmetry guarantees a conservation of the charge and momenta p µ appearing in (3.8) were examples of such charges. Although not every separation of variables can be associated with Killing vector (separation between x and y coordinates found in sections 3 and 4.1 is our prime example), the general theory developed in [18,27,28,34] guarantees that any such separation is related to a symmetry of the background, which is encoded by a (conformal) Killing tensor. In this subsection we will discuss the conformal Killing tensor associated with separation (3.51) and (4.8)-(4.9). The conformal Killing tensor of rank two satisfies equation which is solved in appendix D. Here we will deduce the same solution by using the separation of variables found in section 3. A conformal Killing tensor K M N always implies that is an integral of motion of the massless HJ equation 20 , so K M N can be extracted from the know separation. Going to the eikonal approximation in (4.6) with harmonic function (4.10), we find For separable solutions, both sides of this equation must be constant, and identifying this constant with −I in (4.19), we find In appendix D this expression is derived in a geometrical way by solving the equation (4.18) for the Killing tensor.

Non-integrability of strings
In section 3 we have classified all supersymmetric configurations of Dp-branes that lead to integrable equations for null geodesics. Specifically, we demonstrated that a metric (3.1)-(3.2) leads to a separable HJ equation (1.3) if and only if it has the form with the harmonic function H from (3.50)- (3.52). In this subsection we will investigate whether integrability of geodesics extends to strings with finite size. Our discussion will follow the logic presented in [9], and to compare our results with ones from that paper we rewrite the metric in terms of a new function f = H 1/4 so the metric (4.22) becomes: To demonstrate integrability of strings on a particular background, one has to find an infinite set of integrals of motions, and this ambitious problem has only been solved for very few geometries [35,5]. However, to rule out integrability of sigma model on a given background, it is sufficient to start with a particular solution and look at linear perturbations around it. If such linearized problem has no integrals of motion, one concludes that the original system is not integrable. This approach has been used in [8,9] to rule out integrability of strings on the conifold and on the asymptotically-flat geometry produced by a single stack of Dp-branes. The analysis presented in this subsection is complimentary to [9]: we still focus on the near-horizon limit (where strings are known to be integrable for a single stack), but allow a nontrivial distribution of sources. The equation for linear perturbations around a given solution of a dynamical system is known as Normal Variational Equation (NVE) [36], and to determine whether NVE is integrable, one can use the Kovacic algorithm 21 [37]. Thus to demonstrate that the string theory on a particular background is not integrable one needs to perform the following steps: 4. algebrize NVE (rewrite equations as differential equations with rational coefficients) and transform them to normal form to make NVE be suitable for using the Kovacic algorithm 5. apply the Kovacic algorithm to the obtained NVE, if it fails the system is non-integrable.
Now we apply this method to check integrability of strings in the background (4.23). We begin with looking at the Polyakov action supplemented by the Virasoro constraints For a specific string ansatz on 2-sphere, the system has an effective Lagrangian density Equations of motion for cyclic variables t, φ lead to two integrals of motion (E and ν), and combining this with Virasoro constraint (4.26) we finḋ These equations can be derived from the effective Lagrangian Expanding around a particular solution, we find the following NVE for η ≡ δθ (see appendix E.2 for detail): To proceed we need to choose a particular configuration corresponding to the specific function f . In section 3 we have demonstrated that equations for geodesics are integrable only if function H is given by (3.50)-(3.53), and here we consider (3.51) ignoring the P -term in (3.53) 22 : Here we used the map (3.55) between the coordinates (x, y) and (ρ + , ρ − ). After carrying out all calculations one obtains the following NVE 22 We will discuss NVE associated with this term in appendix E.3.
Application of the Kovacic algorithm to (4.34) shows that the system is not integrable unless d = 0, m + n = 4 (this corresponds to AdS 5 ×S 5 ). Detailed description of the method used in this section and complete calculations are presented in appendices E.1-E.3.
To summarize, we have demonstrated that integrability of geodesics discussed in section 3 does not persist for classical strings, and AdS 5 ×S 5 is the only static background produced by a single type of D-branes placed on a flat base that leads to an integrable string theory 23 .

Geodesics in static Dp-D(p + 4) backgrounds
In the sections 2-3 we analyzed geodesics in the backgrounds produced by a single type of D branes. However, some of the most successful applications of string theory to black hole physics [21] and to study of strongly coupled gauge theories [38] involve intersecting branes, and in this section our analysis will be extended to a particular class of brane intersections. Specifically, we will extend the results of sections 3 to 1/4-BPS configurations involving Dp and D(p + 4) branes. The geometries produced by such "branes inside branes" continue to play an important role in understanding the physics of black holes, and a progress in understanding of the infall problem and Hawking radiation requires a detailed analysis of geodesics and waves on the backgrounds produced by Dp-D(p + 4) systems. In this section we will continue to explore static configurations, and a large class of stationary solutions produced by D1-D5 branes will be analyzed in the next section.
Let us consider massless geodesics in the geometry produced by Dp and D(p + 4) branes: The first and the second terms describe the spaces parallel/transverse to the entire Dp-D(p + 4) system, and the four-dimensional torus represented by dz 2 4 is wrapped by D(p + 4) branes. We assume that Dp branes are smeared over the torus 24 . Metric (5.1) contains two harmonic functions, H 1 and H 2 , which are sourced by Dp and D(p + 4) branes. Away from the sources, these functions satisfy the Laplace equation on the (5 − p)-dimensional flat base with metric ds 2 base . Let us assume that geometry (5.1) leads to a separable Hamilton-Jacobi equation. Then arguments presented in section 3.2 imply that H 2 can only depend on two coordinates, (r 1 , r 2 ), where metric ds 2 base has the form (3.4). To see this, we separate the Killing directions in the action and rewrite the HJ equation (1.3) as 23 Analysis presented in this subsection does not rule out integrability on backgrounds containing NS-NS fluxes in addition to D-branes or on geometries produced by several types of branes, such as Dp-D(p + 4) system discussed in the next section. 24 Some localized solutions are also known [39], but we will not discuss them here.
Here we define If a given distribution of branes corresponds to integrable geodesics, then equation (5.9) should be separable for all allowed values of M and N . Since equation (5.9) is analytic in these parameters, separability must persist even in the unphysical region where M = 0 and N is arbitrary 25 . In this region, equation (5.3) reduces to (3.26) with H → H 2 , p µ p µ → N 2 , then arguments presented in section 3.1 reduce the problem to H 2 (r, θ) with base space Let us now demonstrate that H 1 can only be a function of r and θ. Indeed, separation for M = 0 implies that where y i are coordinates on S m , andỹ j are coordinates on S n . Substituting (5.5), (5.6) into (5.3) and differentiating the result with respect to y k , we find Rewriting this relation as we conclude that H 1 develops unphysical singularities at θ = 0, π 2 for arbitrarily large r unless (∂H 1 /∂y k ) = 0. Similar argument demonstrates that (∂H 1 /∂ỹ k ) = 0, so H 1 can only depend on (r, θ).
To summarize, separability of the HJ equation (5.3) requires the functions H 1 and H 2 to depend only on (r, θ), then the action has the form (3.24), (3.25). This results in the HJ equation for (x, y) and H 2 . As already discussed in section 3.2, solutions (3.51) and (3.52) are related by interchange of two spheres, so without loss of generality, we will focus on (3.50) and (3.51).
Spherical coordinates (3.50) lead to separable equation (5.9) if an only if H 1 and H 2 do not depend on θ, and since these functions have to be harmonic we find Here a = 1 for asymptotically-flat space, and a = 0 for the near-horizon solution. The corresponding metric (5.1) gives the geometry produced by a single stack of Dp-D(p + 4) branes [40]. Separation in the elliptic coordinates (3.51) leads to and now we will determine the corresponding function H 1 . Equation (5.9) separates in coordinates (5.12) if and only if Since H 2 is already given by (5.11), the last relation implies that 26 First we set a = 0 in (5.11), then equation (5.14) becomes To determineṼ 1 (x) andṼ 2 (y), we recall that, away from the sources, function H 1 must satisfy the Laplace equation (C.13): sin n y cos m y ∂H 1 ∂y = 0 (5.16) Substituting (5.15) into (5.16) and performing straightforward algebraic manipulations, we find the most general solutions for H 1 : 26 Notice that M = 0 due to equation (5.9).
(m, n) = (0, 3) Here Π[n; φ|m] is the incomplete elliptic integral. So far we have assumed that a = 0 in H 2 . The case a = 1,Q = 0 corresponds to Dp branes only, so it is covered by discussion in section 3. Solutions with nonzero a andQ can be analyzed by looking at formal perturbation theory in a, and it turns out that H 1 must be constant for such solutions.

Geodesics in D1-D5 microstates
In the last three sections we have analyzed geodesics in a variety of static backgrounds produced by D-branes. In general, supersymmetric geometries are guaranteed to have a time-like (or light-like) Killing vector, so they must be stationary, but not necessarily static. In particular, an interest in stationary geometries produced by the D1-D5 branes has been generated by the fuzzball proposal for resolving the black hole information paradox [22,41]. According to this picture, microscopic states accounting for the entropy of a black hole have nontrivial structure that extents to the location of the naïve horizon, and the black hole geometry emerges as a course graining over such structures. Although the vast majority of fuzzballs is expected to be quantum, some fraction of microscopic states should be describable by classical geometries, and study of this subset has led to important insights into qualitative properties of generic microstates [42].
The fuzzball program has been particularly successful in identifying the microscopic states corresponding to D1-D5 black hole, where all microstate geometries have been constructed in [22,23]. Moreover, a strong support for the fuzzball picture came from analyzing the properties of these metrics [43], and success of this study was based on separability of the wave equation on a special classes of metrics. In this paper we have been focusing on separability of the HJ equation as necessary condition for integrability of strings, but such separability also implies separability of the wave equation. This provides an additional motivation for studying the HJ equation for microscopic states in the D1-D5 system.
In this section we will mostly focus on the HJ equation for particles propagating on metrics constructed in [22], and extension to the geometries for the remaining microstates of D1-D5 black holes [23] will be discussed in the end. The solutions of [22], generalize the static metric (5.1) with p = 1 by allowing the branes to vibrate on the four dimensional base, which is transverse to D1 and D5, and geometries of [23] account for fluctuations on the torus. While the metric (6.1), supplemented by the appropriate matter fields given in [22], always gives a supersymmetric solution of supergravity away from the sources, the bound states of D1 and D5 branes, which are responsible for the entropy of a black hole, have additional relations between H 1 , H 2 and A. Such bound states are uniquely specified by a closed contour F i (v) in four non-compact directions, and the harmonic functions are given by [22] 27 Remarkably, the resulting metric (6.1) is completely smooth and horizon-free in spite of an apparent coordinate singularity at the location of the contour [23]. To avoid unnecessary complications, we will focus on a special case |Ḟ| = 1 (which leads to H 1 = H 2 ≡ H), although our results hold for arbitraryḞ.
Applying the arguments presented in section 3.1 to metric (6.1), we conclude that this geometry must preserve U (1) × U (1) symmetry of the base space, i.e., the profile F i (v) must be invariant under shifts of φ and ψ in 28 This implies that the singular curve, can only contain concentric circles with radii r = (R 1 , R 2 , . . . , R n ) in the θ = π 2 plane and concentric circles with radii r = (R 1 ,R 2 , . . . ,R l ) in the θ = 0 plane. First we focus on circles in the θ = π 2 plane and demonstrate that separability of the HJ equation implies that n = 1. Then we will show that n = 1 also implies that l = 0.

Circles in the
For a single circular contour (r = R 1 , θ = π 2 ), the integrals (6.2) have been evaluated in [24], and superposition of these results gives the harmonic functions for several circles: Here Q i is a five-brane charge of a circle with radius R i , and s i is a sign that specifies the direction for going around this circle.
Let us assume the the HJ equation (1.3) separates in some coordinates. Metric (6.1) has eight Killing directions (t, u, φ, ψ and the torus), they can be separated in the action: Our assumption of integrability amounts to further separation ofS(r, θ) in some coordinates (x, y). In particular, for J ψ = J φ = p u = 0, the HJ equation (1.3) in the metric (6.1) can be written as This equation looks very similar to (3.26), but in practice it is easier to analyze: we don't need to impose the Laplace equation (as we did for (3.26)), since the explicit forms of H and A are known (see (6.4)). The discussion of section 3.2 implies that (6.6) should be viewed as a relation between h(w) (which is defined by (3.42) and (3.45)) and (R 1 , . . . , R n ).
In particular, substitution of the perturbative expansion (C.2) for h(w) leads to an infinite set of constraints on (R 1 , . . . , R n ). To write these constraints, we introduce a convenient notation: Then separation in (k + 2)-rd order of perturbation theory gives a constraint Already the first nontrivial relation (k = 2) implies that R 1 = R 2 = . . . = R n , so it is impossible to have more than one circle (see below). We conclude that the HJ equation does not separate on the background (6.1), (6.4) unless n = 1. The remaining constraints (6.8) for k > 2 are automatically satisfied for this case.
We will now prove that equation (6.8) with k = 2 implies that R 1 = . . . = R n , and the readers who are not interested argument can go directly to part 2. Let us order the radii by R 1 ≥ R 2 ≥ . . . ≥ R n and define a function Then the derivative is positive unless R 1 = . . . = R n and s 1 = . . . = s n , so function G reaches its minimal value when all radii are equal, and it is this value, that gives (6.8) for k = 2. We conclude the equation G = 0 implies that R 1 = . . . = R n .

Circles in orthogonal planes.
Having established that separation requires to have at most one circle in θ = π 2 plane, we conclude that the same must be true about θ = 0 plane, but in principle it is possible to have one circle in each of the two planes. In this case we find and for J ψ = J φ = p y = 0 the HJ equation becomes Fifth order of perturbation theory gives a relation Q 1 Q 2 = 0, which implies that there is no separation in the geometry produced by two orthogonal circles.
The perturbative procedure implemented in part 1 also gives the expression for h(w) in terms of D k for the configurations satisfying (6.8): We have already encountered this holomorphic function in (3.51)-(3.52) (depending on the sign of D 1 ), and demonstrated that it corresponds to elliptic coordinates (3.55) or (C.28). For completeness we present the expression for H 2 that clearly demonstrates the separation of variables in (6.6): So far we have demonstrated that the HJ equation with J ψ = J φ = p u = 0 separates only for microstate whose harmonic functions are given by (6.4) with k = 1, and separation takes place in the elliptic coordinates (6.14), (3.55). A direct check demonstrates that this separation persists for all values of momenta, when the relevant HJ equation is and 29 To summarize, we have demonstrated that the HJ equation (1.3) separates for the stationary D1-D5 geometry (6.1)-(6.2) with |Ḟ| = 1 if and only if the string profile is circular, i.e., the harmonic functions are given by (6.17).
Notice that integrability of (6.16) follows from separability of the wave equation in the background (6.1), (6.17), which has been discovered long time ago [20,43]. Let us clarify the relation between variables used in these papers and the elliptic coordinates (3.55), (6.14).
To separate the wave equation in the metric (6.1) with harmonic functions (6.17), one can use the coordinates (r ′ , θ ′ ) which appear naturally if the D1-D5 solution is viewed as an extremal limit of a black hole [44,20,43]: The relation between (r, θ) and (r ′ , θ ′ ) was found in [24] 30 Looking at the (r ′ , θ ′ ) sector of the metric: we arrive at natural "conformally-Cartesian" coordinates (x, y): Substituting this into (6.19) we find the expressions for (r, θ) in terms of (x, y): r = a sinh 2 x + sin 2 y, cos θ = sinh x cos y sinh 2 x + sin 2 y , (6.22) Using the definitions (3.42), we conclude that w ≡ ln r + iθ = 1 2 ln sinh 2 x + sin 2 y + i arccos sinh x cos y sinh 2 x + sin 2 y = ln(sinh z). (6.23) is a holomorphic function of z = x + iy, as expected from the general analysis presented in section 3.2. Inverting the last expression, we recover the relation (3.54). z = ln 1 2 e w + e 2w + 1 . Thus we conclude that coordinates (6.21) used in [44,43] are essentially the elliptic coordinates up to a minor redefinition of r ′ . We will come back to this feature in section 8.

Geodesics in bubbling geometries.
Given integrability of sigma model on AdS 5 × S 5 , it is natural to look at deformations of this background which might preserve integrable structures. In particular, reference [9] demonstrated that deformation of AdS 5 × S 5 to asymptotically-flat geometry by adding one to the harmonic function destroys integrability of sigma model, although the Hamilton-Jacobi equation for geodesics remains separable. The extension from AdS 5 × S 5 to flat geometry is only possible if one choses flat metric on the worldvolume of D3 branes 31 , and in section 3 we analyzed several classes of geometries produced by flat Dp branes. From the point of view of AdS/CFT correspondence, it is equally interesting to look at field theories on R × S 3 , which are dual to geometries produced by spherical D3 branes. The most symmetric geometry of this type is a direct product of global AdS 5 and a five dimensional sphere, but less symmetric examples are also known [15,45]. In this section we will apply the techniques developed in section 3 to identify the most general geometries of [15] with separable geodesics. We begin with recalling the metrics of 1/2-BPS geometries constructed in [15] 32 : The solutions are parameterized by one function z(x 1 , x 2 , Y ) that satisfies the Laplace equation, and obeys the boundary conditions 3) The regions with z = 1 2 form droplets in (x 1 , x 2 ) plane, and any configuration of droplets leads to the unique regular geometry. Solutions (7.1) are dual to half-BPS states in N = 4 super-Yang-Mills theory, and droplets in (x 1 , x 2 ) plane correspond to eigenvalues of the matrix model describing such states [46] (see [15] for detail).
A generic distribution of droplets in (x 1 , x 2 ) plane leads to solution (7.1), which has a nontrivial dependence upon three coordinates (x 1 , x 2 , Y ). Repeating the arguments presented in section 3.1, one can show that geodesics can only be integrable if at least one of these coordinates corresponds to a Killing direction. Such configurations can be obtained by performing a dimensional reduction of (7.1) along one of the directions in (x 1 , x 2 ) plane. Only two such reductions are possible 33 : This corresponds to field theory living on R 3,1 , which is dual to the Poincare patch of AdS5 32 These metrics are supported by the five-form field strength, and expression for F5 can be found in [15]. Notice that our notation in (F.1) slightly differs from one in [15]: we replaced y of [15] by Y to avoid the confusion with coordinate y introduced in section 3.2. 33 There are also counterparts of (7.4) with ∂rz = 0, which correspond to wedges in the (x1, x2) plane ( figure  4). However, such configurations lead to singular geometries, see [47] for further discussion. The counterpart of (7.5) with ∂2z = 0 is related to (7.5) by rotation.
The complete solution of the Laplace equation (7.2) and expression for V φ for this case were found in [15]: Summation in (7.6) is performed over n circles with radii R i , and following conventions of [15] we will take R 1 to be the radius of the largest circle. For example, a disk corresponds to one circle, a ring to two circles, and so on.
The HJ equation for the solutions specified by (7.6) is analyzed in the appendix F.1, where it is demonstrated that integrability leads to an infinite set of relations between radii R i . Specifically, the expressions defined by must satisfy the relations (F.11): As demonstrated in appendix F.1, this requirement implies that n < 2 in (7.6), so variables separate only for flat space (n = 0) and for AdS 5 ×S 5 (n = 1) ( figure 5(a)). Moreover, construction presented in the appendix F.1 gives the unique set of separable coordinates (F.22) for AdS 5 ×S 5 and their relation with standard parameterization of this manifold will be discussed in section 8.
(b) Geometries with pp-wave asymptotics We will now discuss the geometries with translational U (1) symmetry (7.5), which correspond to parallel strips in the (x 1 , x 2 ) plane (see figure 3(b,c)). It is convenient to distinguish two possibilities: z can either approach different values x 2 → ±∞ (as in figure 3(b) or approach the same value on both sides (as in figure 3(c)). Here we will focus on the first option, which corresponds to geometries with plane wave asymptotics, and the second case will be discussed in part (c).
Pp-wave can be obtained as a limit of AdS 5 × S 5 geometry by taking the five-form flux to infinity [48]. This limit has a clear representation in terms of boundary conditions in (x 1 , x 2 ) plane: taking the radius of a disk ( figure 5(a)) to infinity, we recover a halffilled plane corresponding to the pp-wave (see figure 5(c)). Taking a similar limit for a system of concentric circles ( figure 3(a)), we find excitations of pp-wave geometry by a system of parallel strips (see figure 3(b)). Since strings are integrable on the pp-wave geometry [48], it is natural to ask whether such integrability persists for the deformations represented in figure 3(b). We will now rule out integrability on the deformed backgrounds by demonstrating that even equations for massless geodesics are not integrable.
Solutions of the Laplace equation (7.2) corresponding to the boundary conditions depicted in figure 3(b) were found in [15], and their explicit form is given by (F.24). Such solutions are parameterized by the strip boundaries, and the black strip number i is located at d 2i−1 < x 2 < d 2i . In appendix F.2 we use the techniques developed in section 3 to demonstrate that the HJ equation can only be separable when n = 0 in (F.24), i.e., when the solution represents an unperturbed pp-wave: Interestingly, some special solutions also separate for the pp-wave with an additional strip (see figure 5(d)). Specifically, the geodesics which do not move on the spheres and along x 2 direction (i.e., geogesics with p = 0, L 1 = L 2 = 0 in (F.26)) separate in coordinates (x, y) defined by (F.34): x + iy = w + ln (c) Geometries dual to SYM on a circle Finally we consider configuration depicted in figure 3(c). As discussed in [15], these configurations are dual to Yang-Mills theory on S 3 × S 1 × R, and since we are only keeping zero modes on the sphere, the solutions (7.1) correspond to BPS states in two-dimensional gauge theory on a circle.
The solution of the Laplace equations corresponding to figure 3(c) is given by (F.35), and in appendix F.3, it is shown that only a single strip ( figure 5(b)) leads to a separable HJ equation. For completeness we present the expression for the natural coordinates: x + iy = 2 ln 1 2 e w − (d 2 /l) + e w/2 (7.12) w = ln r l + iθ.
To summarize, we have demonstrated that from the infinite family of the 1/2-BPS geometries constructed in [15], only AdS 5 ×S 5 , pp-wave, and a single M2 brane give rise to integrable geodesics ( figure 5). This implies that the short strings can only be integrable on these backgrounds. Notice, however, that our results do not extend to equations of motion for D3 branes (which are expected to be integrable, as least for 1/2-BPS objects) since the HJ equation (1.3) did not take into account coupling to the RR field. We also found that the separable coordinates (x, y) given by (7.9) are the same elliptic coordinates (or their limits) as the one encountered in section 3.2, and in the next section we will discuss the relation between (x, y) and the standard parameterization of AdS 5 ×S 5 .

Elliptic coordinates and standard parameterization of AdS p ×S q
In the last section we have demonstrated that the HJ equation for geodesics on 1/2-BPS geometries of [15] separates only for AdS 5 ×S 5 and for its pp-wave limit. Moreover, this separation happens in the elliptic coordinates (7.9). On the other hand, equations for supergravity fields on AdS 5 ×S 5 are usually analyzed in the standard parameterization 34 : ds 2 = L 2 − cosh 2 ρdt 2 + dρ 2 + sinh 2 ρdΩ 2 3 + dχ 2 + cos 2 χdφ 2 + sin 2 χdΩ 2 (8.1) and the detailed study of [49] uses the explicit SO(4, 2) × SO(6) symmetry of this metric to separate the resulting equations and to find the mass spectrum of supergravity modes on AdS 5 ×S 5 . This suggests a close relation between the elliptic coordinates and (8.1), which will be clarified in this subsection. We will also discuss the elliptic coordinates for AdS 7 ×S 4 , AdS 4 ×S 7 , and AdS 3 ×S 3 .

1/2-BPS geometries in M-theory.
Let us now turn to the LLM geometries in M-theory [15]: 34 We denoted an azimuthal direction on S 5 by χ to avoid confusion with coordinate θ introduced in (7.4).

Metric (8.5) is parameterized by one function D satisfying the Toda equation,
on a three-dimensional base, and some known boundary conditions in the Y = 0 plane. Although (8.5) is much more complicated than the Laplace equation (7.2), we expect that repetition of the arguments presented in section 7 ensures that function D can only depend on two rather than three variables, and in this case (8.5) can be rewritten as a Laplace equation via a nonlocal change of variables [50]. Specifically, for the rotationally-invariant case, it is convenient to rewrite the metric on the base as 35 FunctionD = D + 2 ln R defined above satisfied the same Toda equation (8.5) as D, and in terms of (X, Y ) this equation can be rewritten as The Toda equation with translational invariance along x 1 can also be written as (8.7)-(8.8) after a replacement (x 1 , x 2 , D) → (φ, X,D). A nonlocal change of coordinates [50,15], Unfortunately, the boundary conditions for V are rather complicated [15] (see [51] for a detailed discussion), and simple expressions for D and V φ similar to (7.6) are not known. Nevertheless, the results presented in sections 6 and 7 strongly suggest that the HJ equation on the geometries (8.5) would only separate on the most symmetric backgrounds: AdS 7 ×S 4 , AdS 4 ×S 7 , and the pp wave. Let us discuss the relation between the standard parameterizations of these backgrounds and the elliptic coordinates. Since elliptic coordinates are only defined in flat space, we begin with rewriting the (X, Y ) sector of (8.7) in a conformally-flat form. It turns out that this is accomplished by going to coordinates (ξ, η) defined by (8.9), and metric (8.7) becomes After introducing the standard polar parameterization (r, θ), ζ = r sin θ, η = r cos θ, (8.12) 35 We introduced an obvious notation: x1 + ix2 = Re iφ , X = ln R.
one can define the elliptic coordinates by (3.55). We will now compare these coordinates with the standard parameterization of AdS 7 ×S 4 and AdS 4 ×S 7 .
Standard parameterization of AdS 3 ×S 3 . As our final example of elliptic coordinates, we consider AdS 3 ×S 3 in global parameterization, which can be obtained by taking the near horizon limit (H → Q/f ) in (6.18) [52]: Rewriting the metric in terms of the elliptic coordinates defined by (6.21): we conclude that these coordinates give the standard parameterization of AdS 3 ×S 3 .
is obtained by taking z = 1 2 x 2 in (7.1) and setting Y = r 1 r 2 , x 2 = 1 2 (r 2 1 − r 2 2 ). (8.27) This leads to a very simple relation for the polar coordinates defined by (7.5) We conclude that in the pp-wave limit, the elliptic coordinate degenerate into the radii of the three-spheres.
The pp-wave limit of AdS p ×S q in M theory, L → ∞, fixed r 2 = Lx, r 5 = L π 2 − y . (8.33) As in the case of the type IIB pp-waves, we conclude that the elliptic coordinates degenerate into the radii of the spheres.
The pp-wave limit of the AdS 3 ×S 3 geometry (8.24) is obtained by writing Since we are dealing with translational rather than rotational symmetry,D = D and X = x2 in (8.9) and sending Q to infinity, while keeping all hatted variables fixed. This results in the metric Once again, elliptic coordinates degenerate to the radii of the one-spheres.
To summarize, we have demonstrated that in all examples of AdS p ×S q , where the HJ equation separates between the sphere and AdS in global coordinates, such separation emerges as a particular case of integrability in elliptic coordinates, and standard parameterization of AdS p ×S q coincides with elliptic coordinates on the relevant flat base. In the pp-wave limits of AdS p ×S q , the elliptic coordinates reduce to the radii of the appropriate spheres.

Discussion
Integrability of geodesics and Klein-Gordon equation has led to numerous insights into physics of black holes. While the black hole solutions are few and far between, the large classes of supersymmetric geometries are known, and in this article we have classified such solutions with integrable geodesics. This integrability is demonstrated to imply that the HJ equation must separate in the elliptic coordinates. For branes with flat worldvolumes, such separation, that extends the known result for the spherical coordinates, can only occur for special distributions of sources, which are analyzed in section 3. For the curved supersymmetric branes, the elliptic coordinates can only be introduced in the most symmetric cases, and as demonstrated in section 7, all these situations reduce to AdS p ×S q or their pp-wave limits.
Our results rule out integrability of N=4 SYM beyond the large N limit. Specifically, we proved that the excitations of strings around heavy supersymmetric states (∆ ∼ N 2 ) are not integrable. While this is consistent with general expectations, it is somewhat surprising that none of the 1/2-BPS geometries give rise to integrable sectors. It would be interesting to extend this result to states with fewer supersymmetries.
Our results also have unfortunate consequences for the technical progress in the fuzzball program. While a large number of geometries corresponding to microscopic states of black holes have been constructed in the last decade, the detailed calculations of the absorption/emission rates have only been performed for the simplest cases. Such calculations are based on solving the Klein-Gordon equation, and as we demonstrated in section 6, this equation, as well as the Hamilton-Jacobi equation for geodesics, cannot separate beyond the known cases. Our results do not imply that a study of geodesics on a particular background is hopeless. A lot of useful information can be extracted by performing numerical integration of the equations of motion and by studying some special configurations rather than generic geodesics.

A Examples of separable Hamilton-Jacobi equations
In this appendix we provide some technical details pertaining to derivation of the results presented in sections 2 and 3. In particular, we write down the explicit expressions for S L on the sphere (A.7), the complete integral for R in the elliptic coordinates (A.16) and make the connection between the holomorphic function introduced in (3.45) and the standard elliptic coordinates (2.12).

A.1 Motion on a sphere
While discussing separation of variables in the Hamilton-Jacobi equation, on several occasions we have encountered equation on a k-dimensional sphere S k (see, for example, (2.7), (3.25)). In this appendix we will write the complete integral of (3.25) using induction. Writing the metric on S k as and splitting S Solving the first equation 38 , and applying induction, we arrive at the complete integral of (A.1) that depends on k parameters L j : This explicit solution should be substituted into (2.5). 38 Although the integral in (A.6) can be performed, the result is not very illuminating.

A.2 Two-center potential and elliptic coordinates
In section 2 we reviewed separation of variables in elliptic coordinates, and here we present some details of that construction. The motion of a particle in the geometry produced by two stacks of Dp branes is governed by the Hamilton-Jacobi equation (2.11) where H is given by (2.9) To rewrite (A.8) in terms the elliptic coordinates (ξ, η) defined by (2.12), we notice that (A.8) can be viewed as a Hamilton-Jacobi equation corresponding to an effective Lagrangian for r and θ: Rewriting the last expression in terms of ξ, η, and going back to the Hamilton-Jacobi equation for S(ξ, η), we find This equation separates in variables (ξ, η) is and only if Rewriting the harmonic function (A.9) in terms of elliptic coordinates, we find that the left hand side of the last expression, separates only for p = 5. In this case equation (A.12) becomes and its complete integral is Here λ is a separation constant and M 2 = −p µ p µ .
To embed the elliptic coordinates in the general framework presented in section 3.2, we have to find the holomorphic function h(w) defined by (3.45). Comparing the kinetic terms in (A.12) with equation (3.38), we find the relation which leads to expressions for (ξ, η) in terms of (x, y): The relation between coordinates (r, θ) and (x, y) looks rather complicated (see (2.12)), but it can be simplified by making use of the complex variables (3.42). First we rewrite equations (A.18) as an expression for z in terms of ξ and η: Next we recall the definition (3.42) of the complex variable w and use (A. 19) to write ρ ± in terms of it: To simplify the expressions for various ingredients appearing in (A.20) it is convenient to define holomorphic functions W ± : Then we find Substitution of these results into (A.20) leads to the desired relation between z and w: Finally, the asymptotic behavior (3.49) determines l in terms of d: l = d/2. As expected from the discussion in section 3.2, z turns out to be a holomorphic function: Finally, by comparing (A.15) with (3.38), (3.39), we extract the expressions for the potentials U 1 (x) and U 2 (y):

B Ellipsoidal coordinates
As demonstrated in section 3.1, if the HJ equation (3.10) separates for k > 2, such separation must occur in ellipsoidal coordinates, including degenerate cases. In this appendix we will demonstrate that a combination of the Laplace equation (3.6) and the requirement (3.5) rules out such separation. To avoid unnecessary complications, we will first give the detailed discussion of the k = 3 case, and in section B.2 we will comment on minor changes which emerge from generalization to k > 3.

B.1 Ellipsoidal coordinates for k = 3
The ellipsoidal coordinates have been introduced by Jacobi [30], and there are several equivalent definitions. We will follow the notation of [53].
Ellipsoidal coordinates, (x 1 , x 2 , x 3 ) are defined as three solutions of a cubic equation for x: where a > b > c > 0 and r 1 , r 2 , r 3 correspond to our r k from (3.4). The roots are arranged in the following order: Coordinates (r 1 , r 2 , r 3 ) can be expressed through (x 1 , x 2 , x 3 ) by 39 In terms of the elliptic coordinates, the metric of the flat three dimensional space becomes Expressions for h 2 , h 3 , R 2 , R 3 are obtained by making a cyclic permutation of indices. To simplify the formulas appearing below, it is convenient to introduce d ij ≡ x i − x j . In ellipsoidal coordinates equation (3.10) becomes Before imposing the Laplace equation (3.6) we will demonstrate that separation requires that (3.4). Indeed, the separation implies that Then multiplying equation (B.4) by d 12 d 13 d 23 and applying ∂ 2 1 ∂ 2 2 to the result, we find a relation which does not involve S: Let us assume that d 1 > 0 in (3.4), then the last relation must hold for all values of L 1 , thus This condition can only be satisfied if a = b, this degenerate case corresponds to oblate spheroidal coordinates: We will now demonstrate that separation of the HJ equation in spheroidal coordinates implies that ∂ φ H = 0, i.e., violation of (3.5) for i = 1, j = 2. This will falsify our assumption d 1 > 0, and similar arguments will show that separability of the HJ equation requires d 2 = d 3 = 0.
To simplify notation and to connect to our discussion in section 3.2 it is convenient to use an alternative form of the oblate spheroidal coordinates (B.8): This gives expressions for the radii: and the metric becomes Notice that these are precisely the elliptic coordinates found in section 3.2. The HJ equation (3.10) in coordinates (x, y) has the form and it does not separate unless Recalling that at large values of x function H does not depend on φ (see (3.7)), we conclude that Φ ′ (φ) = 0, then ∂ φ H = 0 everywhere. As already mentioned, this relation violates the condition (3.5), so our assumption d 1 > 0 was false. Repeating this arguments for the remaining two spheres, we conclude that d 1 = d 2 = d 3 = 0 in (3.4) and (3.6).
We will now combine the integrability condition (B.6), with Laplace equation (3.6) to rule out separability of the HJ equation in ellipsoidal coordinates, Since we have already established that d 1 = d 2 = d 3 = 0, equation (3.6) reduces to the Laplace equation on a flat three-dimensional space formed by (r 1 , r 2 , r 3 ), and it can be rewritten in terms of ellipsoidal coordinates using (B.3): In the remaining part of this subsection we will demonstrate that equations (B.14), (B.16), (3.5), (3.7) are inconsistent with separability of (B.15). Let us assume that the HJ equation (B.15) separates. Although the ellipsoidal coordinates are restricted by (B.2) the formal separation (B.5) must persist beyond this range. In particular, it is convenient to look at the limit x 1 → x 2 , where x 2 is kept as a free parameter. The Laplace equation (B.16) guarantees that H remains finite in this limit, then equation (B.15) reduces to This implies that with the same function F . Repeating this argument for x 1 = x 3 , we conclude that Rewriting (B.15) in terms of F , we conclude that H must be symmetric under interchange of its arguments 40 , and combination of this symmetry with integrability conditions (B.14) leads to severe restrictions on the form of H. Taking (i, j) = (1, 2) in (B.14), and solving the resulting equation, we find where G 1 , G 2 , G 01 , G 02 are some undetermined functions of their arguments. Applying (B.14) with different values of (i, j), we find further restrictions on the form of H: Symmetry of H under interchange of any pair of coordinates implies that G 0i = 0, matrix G ij is anti-symmetric, and all G ij can be reduced to a single function 41 G: Notice that constant and linear functions G lead to H = 0, and quadratic G gives constant H. Substitution of (B.22) into (B.16) leads to a complicated equation for function G: Here we used a shorthand notation: . Equation (B.23) should work around x 2 = x 3 as long as x 1 is sufficiently large 42 , 40 Of course, such formal interchange takes us outside of the physical range (B.2). 41 Specifically, G12 = −G21 = G23 = −G32 = G31 = −G13 so it is convenient to introduce G ≡ G12. 42 All sources of the harmonic function are localized in some finite region of space, which cannot protrude to large x1, which is analogous to radial coordinate. so it can be expanded in powers of d 23 . Taking the leading piece proportional to d 23 , multiplying by d 3 12 , differentiating the result four times with respect to x 1 , we find a closed-form equation for F (x 1 ) ≡ G ′′′ 1 : To analyze this equation, it is convenient to set c = 0 by shifting x i , a 2 , b 2 by −c 2 (see definition (B.3)), and set a = 1 by rescaling x i . The resulting equation (B.24) has a general solution where The low limit of integration, q, will be defined in a moment. Function H should remain finite and smooth as long as x 1 is sufficiently large, this implies that functions G(x 2 ) and G(x 3 ) must be finite for all 0 ≥ x 2 ≥ −b 2 ≥ x 3 ≥ −1 (recall the region (B.2) and our normalization c = 0, a = 1), so function G(x) must be well-defined for all 0 > x > −1. This gives restrictions on C 2 and C 3 .
We will now demonstrate that function G(x) cannot remain finite for x < −b and x > −b unless C 3 = 0. Let us choose q = −b − ε and assume that function G(p) is finite. Recalling the definition of F (F (x) = G ′′′ (x)), we find As x changes from −b − ε to −b + δ, the last three terms as well as contributions proportional to C 1 and C 2 remain finite, so we focus on the term containing C 3 : Here we introduced two functions, which remain finite and non-zero in some vicinity of x = −b. Changing the order of integration in (B.28), we find The integral in the right-hand side does not make sense at x > −b unless which is clearly not the case. We conclude that G 3 (x) (and thus G(x)) is ill-defined at x > −b unless C 3 = 0. 43 Once C 3 is set to zero, a similar analysis in a vicinity of x = √ 1 − b 2 − 1 leads to C 2 = 0, and integrating the resulting F (x), we find Putting back a and c, we can rewrite the last expression as This harmonic function diverges as x 2 and x 3 approach −b 2 unless C 1 = 0.
To summarize, we have demonstrated that separation of the HJ equation in ellipsoidal coordinates with k = 3 implies that the transverse space is three-dimensional, and the harmonic function is given by (B.22) with Such harmonic function, H = C 3 , does not satisfy the required boundary condition (3.7), so ellipsoidal coordinates for k = 3 do not no separate the HJ equation for geodesics in D-brane backgrounds. In the next subsection we will outline the extension of this result to k > 3.

B.2 Ellipsoidal coordinates for k > 3
After giving a detailed description of ellipsoidal coordinates for k = 3, we will briefly comment on the extension of the results obtained in the last subsection to k > 3. The ellipsoidal coordinates in higher dimensions were introduced in [30] 44 : Here r i are the k radii introduced in (3.4), and x i are ellipsoidal coordinates corresponding to k root of an algebraic equation The ranges of ellipsoidal coordinates are analogous to (B.2) in the three-dimensional case: In terms of the ellipsoidal coordinates, the radial part of the metric (3.4), To simplify some formulas appearing below, it is convenient to define Assuming that this equation separates in ellipsoidal coordinates, S(x 1 , . . . , x k ) = S 1 (x 1 ) + . . . + S k (x k ), (B.43) 44 We use a slightly modified notation to connect with discussion in the last subsection.
The last relation is false, and the easiest way to see this is to notice that the leading contribution in the vicinity of x 1 = −a 2 1 comes when all ∂ 1 derivatives hit (x 1 + a 2 1 ) 2 , giving n-order pole, however, the remaining ∂ 2 derivatives do not kill D. As before the degenerate case (e.g., a 2 = a 1 ) requires a separate consideration, and it can be eliminated using the arguments that followed equation (B.8).
Next, we demonstrate that functions S and H must have a formal symmetry under interchange of their arguments. Indeed, taking the limit x 1 → x 2 in (B.42), we find a counterpart of (B.17): This and other similar limits imply that then equation (B.42) ensures the symmetry of H. 45 Modifying the arguments that led to (B.22), we conclude that Indeed, integrability condition (B.44) for (i, j) = (1, 2) gives Symmetry of H implies that C l andC l have the same functional form. Repeating this argument for other pairs (i, j), we find where P l is a polynomial of l-th degree, which is anti-symmetric under interchange of any pair of its arguments. Any such polynomial is proportional tõ which already has degree (k − 1)(k − 2), so P (k−1)(k−2) =P and using the relatioñ we find This proves (B.48).
Since we have already established that d 1 = d 2 = . . . = d k = 0, it becomes easy to write the Laplace equation (3.6) in ellipsoidal coordinates (cf. (B.16)): Substitution of (B.48) into (B.53) leads to a counterpart of equation (B.23): Repeating the analysis which led from equation (B.23) to (B.33) for k = 3, we find the most general solution of (B.54) that exists for all x ∈ (−a 2 j , −a 2 j+1 ): Requiring H to remain finite at sufficiently large x 1 , we conclude that most integration constants in the last relation must vanish, and H must be a constant, just as in the k = 3 case.

C Equation for geodesics in D-brane backgrounds
In this appendix we implement the program outlined in section 3.2. Specifically, we introduce the perturbative expansions (C.2) for functions h(w), U 1 (x) and U 2 (y), substitute the results into (3.47), and find the restrictions imposed by the Laplace equation (3.27).

C.1 Particles without angular momentum
First we analyze the harmonic function (3.27) for L 1 = L 2 = 0, and angular momenta will be added in the next subsection. To stress the fact that we are dealing with this special case, the potentials will be denoted asŨ 1 (x),Ũ 2 (y). Far away from the sources, the harmonic function is given by (3.48), and Laplace equation separates in spherical coordinates, so we find the leading asymptotics 46 : x + iy ≈ ln r l + iθ, H ≈ Q r 7−p . (C.1) In this appendix we will construct the most general function H and coordinates (x, y) which satisfy four conditions: and corresponding potentials are given by 47 We now start with coordinates defined by (C.3) and find the most general separable solution of the Laplace equation (3.27). First we observe that starting with an arbitrary solution (C.3) and adjusting l, we can set a 2 to one of three values (−1, 0, 1). Indeed, solution (C.3), definitions (3.42), and boundary conditions (C.1) remain invariant under the transformation l → e λ l, w → w − λ, x → x − λ, a 2 → a 2 e −2λ (C.5) 46 Due to linearity of the Laplace equation (3.27), a constant term can always be added to H, so the boundary condition (C.1) can be easily extended to the asymptotically-flat case H ≈ a + Q r 7−p . 47 For p ≥ 5 we find additional solutions forŨ1 andŨ2, which will be discussed below.
for any real λ. Since a 2 is a real number 48 , it can be set to zero or to ±1 by choosing an We recall that in this section we are focusing on the special case L 1 = L 2 = 0.
(b) We will now find the most general functionŨ (r, θ) =Ũ 1 (x) +Ũ 2 (y) which satisfies the Laplace equation (3.27) in the metric ds 2 9−p = dr 2 + r 2 dθ 2 + r 2 sin 2 θdΩ 2 n + r 2 cos 2 θdΩ 2 m . (C.11) Using a relation we can rewrite (C.11) in terms of x and y: ds 2 9−p = d 2 (cosh 2 x − cos 2 y)(dx 2 + dy 2 ) + sinh 2 x sin 2 ydΩ 2 n + cosh 2 x cos 2 ydΩ 2 m (C. 12) and to rewrite the Laplace equation (3.27) in these coordinates: sinh n x cosh m x ∂H ∂x + 1 sin n y cos m y ∂ ∂y sin n y cos m y ∂H ∂y = 0 (C.13) We are looking for a separable solution of equation (C.13) which has the form (C.10): Since the first line of this equation is the sum of x-and y-dependent terms, the second line must have the same structure, so it can be rewritten as
This change of variables can be obtained from case 2 by making replacements 50 Repeating the steps which led to (C.25), we arrive at while the counterpart of (C.26) becomes The corresponding harmonic function is given by (C.34).

C.2 General case
As shown in the last subsection, the harmonic function (3.47), solves the Laplace equation (3.27) as long as L 1 = L 2 = 0 and U 1 (x), U 2 (y) are given by (C. 25) or (C.26). We will now find the potentials U 1 (x) and U 2 (y) for non-zero values of L 1 and L 2 .
Recalling the change of coordinates (C.9), we can write the first two terms in (C.31) as Then the harmonic function (C.31) can be rewritten as whereŨ 1 andŨ 2 are given by (C.25) or (C.26). Reparametrization (C.28) can be obtained from (C.9) by making the replacement (C.27), then instead of (C.33) we find 51 In this case, functionsŨ 1 andŨ 2 are given by (C.29) or (C.30).

D Equations for Killing tensors
In section 4.2 we discussed the Killing tensors associated with separation of the HJ equation in elliptic coordinates. The detailed calculations are presented in this appendix.

D.1 Killing tensors for the metric produced by D-branes
In this appendix we will solve the equations for Killing tensors in the geometry (4.2)-(4.3) and derive the expression (4.21). To simplify some formulas appearing below, we modify the parameterization of (4.2)-(4.3): We begin with analyzing (xxx) and (yyy) components of this relation: Substituting this into (xxy) component and recalling that g xx = g yy , we find Using the explicit expressions for the Christoffel symbols, equation (D.5) can be rewritten as Combining this with a similar relation coming from (xyy) component of (D. 3), we conclude that so K xy and N ≡ K xx − K yy are dual harmonic functions. We also quote the expressions for W i , which can be obtained by evaluating the covariant derivatives in (D.4): (D.10) Next we look at the Killing equation (D.3) with three legs on R 1,p : Using the expression for the relevant Christoffel symbol, we can rewrite (D.11) as an equation for a conformal Killing tensor on R 1,p : Tensors satisfying this equation have already been encountered in section 3.1, and the general solution of (D.13) has been found in [29]: Here K µ are responsible for separation on R 1,p (which is already taken into account in the ansatz (3.8)), and to study the separation in (x, y) coordinates 52 , it is sufficient to focus on the first term in (D.14), which is invariant under ISO(p, 1) transformations. Repeating this arguments for equation (D.3) with three legs on the spheres, we conclude that the Killing tensor responsible for separation of x and y is invariant under SO(m + 1) × SO(n + 1) × ISO(p, 1), in particular, this implies that Here h ab and hȧ˙b are metrics on S m and S n . Ansatz (D.16) ensures that all components of equation (D.3) which do not contain legs along x or y are satisfied.
Next we look at (x, µ, ν) and (y, µ, ν) components: Taking into account the expressions for the Christoffel symbols, we can rewrite equation (D.17) as Similar analysis of equations along sphere directions gives: In the special case K xx = K xy = 0 (which corresponds to the Killing tensor (4.21)), we find Integrability conditions require that

D.2 Review of Killing tensors for flat space
In the last subsection we have encountered Killing equations on R 1,p space and on a sphere, and we chose the simplest solution, 25) In this section we review the most general solution of (D.24). Killing tensors on symmetric spaces were studied in [29], where it was shown that the most general solution of (D.24) on a sphere can be written as µ are conformal Killing vectors satisfying equation Although an explicit solution of equation (D.27) was not given in [29], it can be easily derived, and in this section we will present such derivation for flat space since this case played an important role in the construction presented in section 3.1 53 . Out result can be easily extended to the conformal Killing tensors on a sphere. Consider the equation for the Conformal Killing Vector (CKV) on flat space: From equations with i = j we conclude that Then applying ∂ 1 ∂ 2 to the (1, 2) component of equation (D.28), we find a restriction on f : In two dimensions this exhausts all equations for the CKV, thus the most general solution is parameterized by one holomorphic function F :  33) we conclude that Repetition of this argument for all i = j gives Substituting this into (D.31), we find a relation between D 1 and D 2 : Similar arguments show that To summarize, we have demonstrated that f is a linear function, and (D.29) leads to the following expressions for K 1 , K 2 : Substitution of these relations into (D.30) gives an equation for h 1 and h 2 : 39) and the solution reads Repeating this argument for other pairs, we find the most general expression for K i : To summarize, we have demonstrated that any conformal Killing vector in flat space has form (D.41) with constant parameters A i , B, C ij . The arguments of [29] imply that the general Killing tensor on such space can be written as combination of such Killing vectors (D.26) with additional constants A, B ab .

E Non-integrability of strings on D-brane backgrounds
In section 3 we have classified Dp-brane backgrounds leading to integrable geodesics, and in section 4.3 we outlined the procedure for studying dynamics of strings on such geometries. This appendix contains technical details supporting the analysis of section 4.3. Appendix E.1 reviews one of the approaches to integrable systems, and section E.2 applies this general construction to specific configurations of strings on Dp-brane backgrounds. In particular, in appendix E.3 we focus on the geometries leading to integrable geodesics 54 and rule out integrability of strings for the majority of such backgrounds.

E.1 Review of analytical non-integrability
A dynamical system is called integrable if the number of integrals of motion is equal to the number of degrees of freedom. To demonstrate that a system is not integrable, it is sufficient to look at linear perturbations around a particular solution and to demonstrate that the resulting linear equations do not have a sufficient number of integrals of motion. The equation for linear perturbations is known as the Normal Variational Equation (NVE), and the Kovacic algorithm 54 As discussed in the Introduction, integrability of geodesics is a necessary condition for integrability of strings.
[37] gives a powerful analytical tool for studying this equation [36]. Let us outline the procedure developed in [36]. Consider a Hamiltonian system parameterized by canonical variables as (x i , p i ), (i ranges from 1 to N ). To rule out integrability using NVE, one has to perform the following steps.

Start with a Hamiltonian and write down the equations of motion
2. Write the full variational equations of (E.2): 3. Choose a particular solution x i =x i , p i =p i , i = 1, ..., N and one pair of the canonical coordinates, (x k , p k ). The Normal Variational Equation (NVE) is a subsystem of (E.3) with i = k: 4. Rewrite the system (E.4) as a second-order differential equation for δx k : δẍ k + q(t)δẋ k + r(t)δx k = 0, (E.5) or in the standard notation η + q(t)η + r(t)η = 0, δx k = η. (E.6) Equation (E.6) is known as NVE [36].
5. Make the NVE suitable for using the Kovacic algorithm, namely, algebrize the NVE (i.e. rewrite it as differential equation with rational coefficients) by using a change of variables t → x(t): The Kovacic algorithm requires the coefficients in this equation to be rational functions of x. The next step is to convert (E.7) to the normal form by redefining function η: . Apply the Kovacic algorithm to (E.8), and if it fails then the system is not integrable.
A more detailed explanation of these steps, containing several examples, can be found in [9].

E.2 Application to strings in Dp-brane backgrounds
Let us now apply the NVE method described in Appendix E.1 to equations for strings in the Dp-brane background given by (3.21): To compare with [9] we introduce a notation f = H 1/4 . The Polyakov action for the string, leads to the equations of motion and Virasoro constraints, To demonstrate that system is not integrable, it is sufficient to look only at a specific sector and show that the number of conserved quantities does not match the number of variables. Specifically, we consider the following ansatz: All other coordinates are considered to be constants. Substitution (E.13) into (E.10) leads to the Lagrangian L = −f −2ṫ2 + f 2ṙ2 + f 2 r 2 (− sin 2 θφ ′ 2 +θ 2 ). (E.14) Solving the equations of motion for cyclic variables t, φ t = Ef 2 , , n + m ≥ 2 (E. 25) where we used the mapping between (x, y) and (ρ + , ρ − ), namely The NVEs for P = 0 are the same in both cases (E.24), (E.25): Application of the Kovacic algorithm shows non-integrability of the system unless d = 0 and m + n = 4. The integrable case corresponds to strings on AdS 5 × S 5 , then equation (E.27), coincides with (3.26) from [9]. Non-vanishing P -term in (E.24) for n = 0, m = 0 and n = 0, m = 1 gives rise to divergences at θ = π/2 (f (r, π/2) = 0), so the NVE around solution (E.20) is not well defined. The last remaining case of (E.24) is n = 1, m = 0. SettingQ = 0 gives the following NVE Using the Kovacic algorithm we see that this term corresponds to an integrable system. Unfortunately the corresponding geometry is unphysical since it has singularities at arbitrarily large values of r. To see this we recall that metric (E.9) becomes singular when f = 0, and function f vanishes at ρ + − ρ − = 2d, which corresponds to θ = 0 and arbitrary value of r (recall (2.9)). To summarize, in this appendix we have demonstrated that supergravity backgrounds with integrable geodesics do not lead to integrable string theories, with the exception of AdS 5 × S 5 .

F Geodesics in bubbling geometries.
In this appendix we will present the analysis of geodesics in the 1/2-BPS geometries constructed in [15]. Our conclusions are summarized in section 7.
The BPS geometries constructed in [15] are supported by the five-form field strength, and the metric is given by where function z(x 1 , x 2 , Y ), which satisfies the Laplace equation, obeys the boundary conditions As discussed in section 7, only three classes of configurations can potentially lead to integrable geodesics, and we will discuss these classes in three separate subsections.

F.1 Geometries with AdS 5 ×S 5 asymptotics
The boundary conditions depicted in figure 3(a) lead to geometries (F.1) which are invariant under rotations in (x 1 , x 2 ) plane, and such solutions are conveniently formulated in terms of coordinates introduced in (7.4): The complete solution of the Laplace equation (F.2) and expression for V φ for this case were found in [15]: Summation in (F.5) is performed over n circles with radii R i , and following conventions of [15], we take R 1 to be the radius of the largest circle. For example, a disk corresponds to one circle, a ring to two circles, and so on.
To write the HJ equation (1.3) in the metric (F.4), we notice that coordinates (t, φ), as well as two spheres, separate in a trivial way, this gives S = −Et + Jφ + S We begin with analyzing equations for geodesics with L 1 = L 2 = J = 0: which is similar to (3.26). To evaluate the last term in (F.8), we need expressions (F.5) as well as relation Applying the techniques developed in section 3 to (F.8), we conclude that this equation is separable if and only if there exists a holomorphic function g(w) = x + iy, such that Complex variable w is defined in terms of (r, θ) by (3.42). To find g(w), we employ the same perturbative technique that was used to derive (3.50)-(3.52): starting with a counterpart of (C.2), g(w) = w + k>0 a k e −kw , U 1 (x) = e −2x 1 + k>0 b k e −kx , (F. 10) and solving (F.9) order-by-order in e −x , we find the necessary condition for the existence of a solution: (D 2 ) k−1 D 2(k+1) = (D 4 ) k . (F.11) Here we defined combinations of R j appearing in (F.5): Relations (F.11) are clearly satisfied for n = 0 (flat space) and for n = 1 (AdS 5 ×S 5 ), and we will now demonstrate that they fail for n ≥ 2.
If n ≥ 2 in (F.12), we can define ε ≡ R 2 /R 1 < 1 and (F.14) Rewriting an infinite set of relations (F.11) as Then inequality (F.14) implies that The left-hand side of this inequality does not depend on k, and the right-hand side goes to zero as k goes to infinity (n is an arbitrary but fixed number), then we conclude that Substituting this into (F.15), we express allD k throughD 2 and R 1 : then inequality (F.14) ensures thatD 2 = 0: On the other hand, the definition (F.13) implies thatD 2 is strictly negative 55 , as long as n ≥ 2. We conclude that relation (F.11) holds only for n = 0, 1. To show this for all values of n we used an infinite set of relations (F.11), however, for any given n is it sufficient to use (F.11) with k = 0, . . . , n. Note that the equations with k = 0, 1 are trivial. Going back to the solution with n = 1, we can extract the relevant holomorphic function and potential U 1 (x) (see (F.9)) g(w) = ln 1 2 e w + e 2w − (R 1 /l) 2 , U 1 (x) = e 2x E 2 R 2 1 l 2 (e 2x + (R 1 /2l) 2 ) 2 , U 2 (y) = 0. (F.21) We recall the w is defined by (3.42), and this relation has a free dimensionful parameter l, which can be chosen in a convenient way. Setting l = R 1 /2, we recover the standard elliptic coordinates (3.51), which can also be rewritten as (3.55): The relation between (x, y) and standard coordinates on AdS 5 ×S 5 will be discussed in section 8.

F.2 Geometries with pp-wave asymptotics
We will now discuss the geometries with translational U (1) symmetry (7.5), which correspond to parallel strips in the (x 1 , x 2 ) plane (see figure 3(b,c)). It is convenient to distinguish two possibilities: z can either approach different values x 2 → ±∞ (as in figure 3(b)) or approach the same value on both sides (as in figure 3(c)). In this subsection we will focus on the first option, which corresponds to geometries with plane wave asymptotics, and the second case will be discussed in subsection F.3. Pp-wave can be obtained as a limit of AdS 5 × S 5 geometry by taking the five-form flux to infinity [48]. This limit has a clear representation in terms of boundary conditions in (x 1 , x 2 ) plane: taking the radius of a disk to infinity, we recover a half-filled plane corresponding to the pp-wave (see figure 5(c)). Taking a similar limit for a system of concentric circles, we find excitations of pp-wave geometry by a system of parallel strips (see figure 3(b)). Since strings are integrable on the pp-wave geometry [48], it is natural to ask whether such integrability persists for the deformations represented in figure 3(b). In this subsection we will rule out integrability on the deformed backgrounds by demonstrating that even equations for massless geodesics are not integrable.
Solutions of the Laplace equation (F.2) corresponding to the boundary conditions depicted in figure 3(b) are given by [15] (F.24) Here n is the number of black strips (strip number i is located at d 2i−1 < x 2 < d 2i ), and x 2 = d 0 is the boundary of the half-filled plane. We will assume that the set of d j is ordered: Although one can shift x 2 to set d 0 = 0, we will keep this value free and use the shift symmetry later to simplify some equations.
Repeating the steps performed in the last subsection, we write the counterpart of (F.6) S = −Et + px 1 + S L 1 (y) +S where coordinates (r, θ) are defined by (7.5). The HJ equation is given by the counterpart of (F.7): As before, we begin with analyzing this equation for L 1 = L 2 = p = 0: Results of section 3 ensure that equation (F.28) is separable if and only if there exists a holomorphic function g(w), such that where complex variable w is defined by (3.42) and x + iy = g(w). Solving equation (F.29) in perturbation theory, we find relations for d i , which are much more complicated that (F.11).
Introducing the convenient notation D p = 2n j=1 (−1) j (d j ) p as in the previous subsection we find two options for the relations between D i :

F.3 Geometries dual to SYM on a circle
Finally we consider configuration depicted in figure 3(c). As discussed in [15], these configurations are dual to Yang-Mills theory on S 3 × S 1 × R, and since we are only keeping zero modes on the sphere, the solutions (F.1) correspond to BPS states in two-dimensional gauge theory on a circle. Following the discussion presented in the last subsection, we arrive at equations (F.27) and (F.28), however, instead of (F.24) we find where n is the number of black strips (strip number i is located at d 2i−1 < x 2 < d 2i ), and we assume that the set of d j is ordered: so by choosing an appropriate α, we can set D 2 = 0. This choice leads to great simplifications in (F.38), in particular, we find that D 4 = 0. We will now show that D 2p = 0 for all values of p.
First we observe that the structure of our perturbative expansion guarantees that restriction at order p − 1 has the form where P is some polynomial of p arguments. Next we notice that if any configuration of strips leads to a separable Hamilton-Jacobi equation, the configuration which is obtained by the reflection about x 2 axis must have the same property. The reflection corresponds to the transformation so if equations (F.41) are solved by the set {D k }, they should also be solved by the set {D ′ k }. We can now use induction to demonstrate that any solution of (F.41) with D 2 = 0 must have D 2k = 0 for all k. (F.43) The statement is trivial for k = 1, and we assume that it holds for all k < K.  so equations (F.45), (F.51) are very similar to equations (F.12), (F.11), and we can used the same logic 59 to conclude that n = 1.
To separate equation (F.28) for one strip, it is convenient to set d 1 = 0 by shifting the origin of x 2 . Then we find g(w) = 2 ln 1 2 e w − (d 2 /l) + e w/2 , U 1 (x) = 8d 2 le x E 2 (d 2 + 4le x ) 2 , U 2 (y) = 0. (F.53) Now we go back to the original HJ equation (F.27) containing non-vanishing momenta L 1 , L 2 , p and demonstrate that momenta do not spoil separability. Here we consider the separable case of one strip.
Based on the logic used through the entire paper separability of each momentum p, L 1 , L 2 requires r 2 |g ′ (w)| 2 2pV 1 E + p 2 + Expressing the left-hand side of (F.54) in terms of (x, y) and using the holomorphic function g(w) from (F.53), we find Thus addition of momenta does not spoil separability.
All results obtained in this appendix can be summarized by listing all bubbling geometries leading to separable geodesics (here we use the language introduced in [15] to describe the solutions):