Exact solutions and critical chaos in dilaton gravity with a boundary

We consider (1+1)-dimensional dilaton gravity with a reflecting dynamical boundary. The boundary cuts off the region of strong coupling and makes our model causally similar to the spherically-symmetric sector of multidimensional gravity. We demonstrate that this model is exactly solvable at the classical level and possesses an on-shell SL(2,R) symmetry. After introducing general classical solution of the model, we study a large subset of soliton solutions. The latter describe reflection of matter waves off the boundary at low energies and formation of black holes at energies above critical. They can be related to the eigenstates of the auxiliary integrable system, the Gaudin spin chain. We argue that despite being exactly solvable, the model in the critical regime, i.e. at the verge of black hole formation, displays dynamical instabilities specific to chaotic systems. We believe that this model will be useful for studying black holes and gravitational scattering.


Introduction
The models of two-dimensional dilaton gravity were popular for decades [1][2][3]. Some of them describe spherically-symmetric sectors of multidimensional gravities with dilaton fields φ related to the sizes of the extra spheres 1 . Some others are exactly solvable at the semiclassical [4,5] or quantum [3] levels which makes them valuable for studying black holes and gravitational scattering [6][7][8].  These models become particularly important in the context of information paradox [9,10] confronting an apparent loss of quantum coherence during black hole evaporation with the principles of quantum theory. Since unitarity of quantum gravity is strongly supported by the AdS/CFT correspondence [11,12], modern AMPS argument [13,14] suggests dramatic violation of the equivalence principle ("firewalls") in the vicinity of old black hole horizons, see [15,16] for earlier works. This feature, if exists, may leave "echoes" in the gravitational wave signal [17,18] to be detected by LIGO [19,20], cf. [21,22]. From the theoretical viewpoint, further progress can be achieved by understanding unitary evolution of black holes outside of the explicit AdS/CFT framework. This brings us to the arena of two-dimensional models which may, in addition, clarify relation of black holes to quantum chaos [23][24][25][26][27][28], cf. [29].
Unfortunately, solvable models of two-dimensional dilaton gravity essentially differ from their multidimensional cousins. Consider e.g. the celebrated Callan-Giddings-Harvey-Strominger (CGHS) model [4], see [1,2] for reviews. Its two-dimensional Minkowski vacuum in Fig. 1a, unlike the multidimensional vacua, has disconnected sets of "left" and "right" infinities J ± L and J ± R , and transitions between those are expected [30] to be important for the information loss problem. Besides, the CGHS model is strongly coupled [31] near the "left" infinities which puts its semiclassical results on shaky ground. It was recently suggested [32] that due to the above peculiarities evaporation of the CGHS black holes leads to remnants rather than firewalls.
We consider the modified CGHS model proposed 2 in [33,34], see also [31,[39][40][41][42]. The region of strong coupling in this model is cut off by the reflective dynamical boundary placed at a fixed value φ = φ 0 of the dilaton field, see Fig. 1b. Param-eter e 2φ 0 1 plays the role of a small coupling constant. We explicitly obtain reparametrization-invariant action of the model by restricting CGHS action to the space-time region φ < φ 0 and adding appropriate boundary terms. Note that the original CGHS model is formally restored in the limit φ 0 → +∞ which shifts the regulating boundary in Fig. 1b all the way the left. We do not consider this limit avoiding potential problems with strong coupling, cf. [43][44][45].
As an additional bonus, the above model with a boundary is causally similar to spherically-symmetric multidimensional gravity, cf. Fig. 1b. The price to pay, however, is nonlinear equation of motion for the boundary which, if non-integrable, may damage major attractive property of the CGHS model -its solvability. Note that the previous studies of this or similar models were relying on numerical [40][41][42]46] or shock-wave [33,34,39] solutions.
In this paper we demonstrate that the CGHS model with a boundary is exactly solvable at the classical level. We obtain general solution of the classical field equations and construct an infinite number of particular soliton solutions. The latter describe reflection of matter waves off the boundary at low energies and formation of black holes at energies above some critical values, see Figs. 2a and 2c. Each solution is characterized by N integers or half-integers s 1 , . . . , s N and the same number of real parameters. The parameters of the solitons satisfy inequalities ensuring positivity of energy.  We establish one-to-one correspondence between the above solitons and the eigenstates of the auxiliary integrable system -the rational Gaudin model [47][48][49]. This allows us to classify these solutions and study their properties.
We find that equation of motion for the boundary is invariant under conformal transformations v → w(v), where v is the light-cone coordinate, w(v) is an arbitrary function. These transformations relate physically distinct solutions, and one should not confuse them, say, with the residual reparametrization symmetry in [39,40]. In particular, the transformations from the global SL(2, R) subgroup change massless matter field(s) f of the model as the standard zero-weight fields. They also map the solitons into solitons. The transformations with nonzero Schwarzian derivative act non-linearly on f , and we do not consider them in detail.
Finally, we study dynamics of the model in the critical regime, i.e. at the verge of black hole formation, cf. Fig. 2b. We demonstrate that in this limit scattering of matter waves off the boundary displays instabilities specific to chaotic systems: the final state of the process becomes extremely sensitive to the initial Cauchy data. This feature is in tune with the near-horizon chaos suggested in [25]. We argue that it impedes global integrability of the model, i.e. prevents one from choosing a complete set of smooth conserved quantities in the entire phase space.
In Sec. 2 we introduce dilaton gravity with a boundary and study its properties. We construct exact solutions in Sec. 3. Critical chaos is considered in Sec. 4. In Sec. 5 we discuss possible applications of our results.
where 3 the integrand in the first line is the CGHS Lagrangian [4] describing interaction of the metric g µν and dilaton φ with massless scalar f ; the dimensionful parameter λ sets the energy scale of the model. In Eq. (2.1) we modified the CGHS action by restricting integration to the submanifold φ < φ 0 and adding the boundary terms 4 at φ = φ 0 . We introduced the proper time of the boundary τ , its extrinsic curvature K = g µν ∇ µ n ν , and unit outer normal n µ ∝ ∇ µ φ.
In fact, the choice of the boundary action in Eq. (2.1) is limited. First, the Gibbons-Hawking term with extrinsic curvature ensures consistency of the gravitational action. Without this term the boundary conditions following from Eq. (2.1) would be incompatible with the Dirichlet condition φ = φ 0 , see [51] and cf. Appendix A.1. Second, we assume no direct interaction of the matter field f with the boundary. Then the only natural generalization of our model would include an arbitrary constant in the last term of Eq. (2.1). However, this parameter needs to be fine-tuned in order to retain Minkowski solution (see below). Thus, the action (2.1) describing interaction of the boundary with the gravitational sector of the CGHS model is fixed [33].
The quantity e 2φ 0 is a coupling constant controlling loop expansion in the model (2.1). Indeed, change of variablesφ = φ − φ 0 ,f = e φ 0 f brings this parameter in front of the classical action, S =S/e 2φ 0 . Thus, e 2φ 0 plays the role of a Planck constant implying that the model is classical at e 2φ 0 1. It is clear that the bulk equations in the model (2.1) are the same as in the original CGHS model [4] [1,2]. However, extremizing the action with respect to the boundary values of g µν and f , one also obtains the boundary conditions satisfies the CGHS equations, cf. Appendix A.1. In this case the boundary φ = φ 0 is static, x boundary = −φ 0 /λ, and the first of Eqs. (2.2) is automatically satisfied. Note that the Minkowski vacuum (2.3) is a solution in our model due to exact matching between the bulk and boundary terms with λ in the action (2.1).

Solution in the bulk and reflection laws
The CGHS equations in the bulk are exactly solvable [1,2] in the light-cone frame (u, v), where ds 2 = −e 2ρ(u, v) dudv . In the frame (2.4) the matter field satisfies ∂ u ∂ v f = 0 and therefore splits into a sum of incoming and outgoing parts, The respective energy fluxes are This specifies the Cauchy problem in our model: one prepares f in or T vv at the past null infinity and calculates f out or T uu at J + , see Fig. 1b. The solution for the scale factor ρ and dilaton field φ is where We fixed the integration constants in these expressions by requiring, first, that the space-time is flat in the infinite past, i.e. no white hole preexists the scattering process. Second, we chose the coordinates in such a way that the quadrant u ∈ (−∞; 0), v ∈ (0; +∞) covers all space-time accessible to the distant observer. In particular, the limits u → −∞ at v > 0 and v → +∞ at u < 0 lead to J − and J + , respectively, see Fig. 3. Now, consider the boundary φ = φ 0 described by the function u = U (v) in the "Kruskal" coordinates. Substituting the bulk solution (2.5), (2.7) into the boundary conditions (2.2), one obtains equation for U (v) and reflection law for the matter field f , see Appendix A.2 for the derivation of these equations and proof that they are compatible with the definition φ(U (v), v) = φ 0 of the boundary. Note that the second of Eqs. (2.9) relates the incoming and outgoing waves by conformal transformation v → U (v). The first equation implies that the boundary is always time-like, dU/dv > 0. When rewritten in the appropriate terms, it coincides 5 with the boundary equation obtained in [33,34,39] using energy conservation. One easily finds solution in the empty space using Eqs. (2.9) and (2.7) with T vv = T uu = 0, where the integration constant in the first expression was chosen to make U (v) smooth and invertible in the interval 0 < v < +∞. Solution (2.10) is the linear dilaton vacuum 6 : coordinate transformation brings it to the standard form (2.3). In what follows we impose flat asymptotics (2.10) in the infinite past v → 0, u → −∞.
Note that the space-time (2.7) is always flat far away from the boundary, i.e. at large |u| and v. Below we transform to the asymptotic Minkowski coordinates (t, x) using Eq. (2.11).
We have got a receipt for solving the Cauchy problem in the CGHS model with a boundary. In this case the initial Cauchy data are represented by the incoming wave f in (v) or its energy flux T vv (v). One solves Eqs. (2.9) with the initial condition (2.10) at v → 0 and finds U (v), f out (u). The scale factor of the metric, dilaton and matter fields are then given by Eqs. (2.7) and (2.5).

Simple equation for the boundary
One notices that Eq. (2.9) for U (v) is, in fact, a Riccati equation. The standard substitution brings it to the form of a Schrödinger equation for the new unknown ψ(v), Note that ψ(v) is defined up to a multiplicative constant. Now, one can solve for ψ(v) given the initial data T vv (v). After that the entire solution is determined by Eq. (2.12) and expressions from the previous Section. For example, the outgoing energy flux equals , (2.14) 5 It does not conform, however, with the boundary conditions introduced at one-loop level in [43][44][45]: in the classical model the latter conditions imply that the boundary is space-like. 6 Recall that we excluded solutions with eternal black holes in Eq.
We obtained Eq. (2.14) by substituting the reflection law (2.9) into the definition (2.6) of the flux and then expressing the derivative of U (v) from the first of Eqs. (2.9) and Eq. (2.12). Importantly, Eq. (2.13) is well-known in mathematical physics. Similar equation appears in Liouville theory at classical and semiclassical levels [52]. Besides, the eigenstates of the Gaudin model [47] can be related to the solutions of Eq. (2.13) with monodromies ±1 and rational T vv (v) [48]. In what follows we exploit these similarities for studying exact solutions in dilaton gravity.
The function ψ(v) in Eq. (2.12) has simple geometric meaning. First, the value of ψ is related to the proper time τ along the boundary, where we used Eqs. (2.4), (2.9), (2.12) and introduced the arbitrary constant ψ 0 related to the origin of τ . Function τ (v) is illustrated in Fig. 3. Second, recall that v is the exponent of the flat light-cone coordinate (t + u) far away from the boundary, Eq. (2.11). Thus, ψ(v) maps the affine coordinate at J − to τ . Equation (2.13) relates this coordinate-independent function to the asymptotic Cauchy data T vv (v).  Consider general properties of classical solutions in the model with a boundary. Expression (2.15) implies that ψ(v) vanishes in the infinite past, Indeed, behavior ψ → c 0 v as v → 0 corresponds to the linear dilaton vacuum (2.10) in the beginning of the process. To simplify the next argument, we set 7 c 0 = 1. We consider well-localized T vv (v) and therefore linear asymptotics of the solution to Eq. (2.13). If T vv is small, one has C ≈ 1. The respective "lowenergy" solutions describe reflection of matter waves off the time-like boundary, see Figs. 4a,b. As T vv grows, the function ψ(v) becomes more concave and C decreases because ∂ 2 v ψ ∝ −T vv < 0. For some large fine-tuned T vv (v) one obtains critical solutions with C = 0. In this case the boundary is null in the asymptotic future because its proper time τ (v) in Eq. (2.15) remains finite as v → +∞. The respective "critical" solution in Figs. 4 is at the brink of black hole formation: we will see that the asymptotically null boundary sits precisely at the horizon of would-be black hole.
At sufficiently high energies we get C < 0 and therefore ψ(v) has a maximum (point A in Fig. 4a). The boundary is null at this point: where Eqs. (2.9), (2.12), and (2.7) were solved to the leading order in Besides, one discovers that the condition φ = φ 0 defines two intersecting curves 3 near A, and only one of those is the time-like boundary considered so far. The second curve is space-like, it is shown by the dashed line in Fig. 4b. The boundary conditions (2.9) are not met at this line. We obtained the analog of the black hole singularity in the model with a boundary. Indeed, our model is formulated at φ < φ 0 i.e. in the space-time region to the right of both solid and dashed graphs in Fig. 4b. The space-like "edge" φ = φ 0 swallows all matter at u > 0 limiting the region accessible to the outside observer to u < 0. The line u = 0 is a horizon.
Except for the point A itself, the solution is smooth at the space-like "singularity" φ = φ 0 . This fact was not appreciated in the previous studies. The mass of the formed black hole, by energy conservation, is related to the value of the dilaton field at the future horizon, where we subtracted the final matter energy from the initial one in the first equality (cf. Eq. (2.11)), integrated by parts and used Eqs. (2.8) in the second equality, and then expressed the result in terms of φ, Eq. (2.7). Since φ < φ 0 , this implies that all black hole masses are larger than see detailed discussion in [5,42]. Black holes with M bh = M cr have boundary sitting precisely at the horizon. They are formed in the critical solutions. The solutions in Fig. 4b, when replotted in the finite-range coordinates (ū,v) = (arctan u, arctan v), look like Penrose diagrams, see Fig. 2. From now on, we will exploitū andv for visualizing the solutions. We will also mark the (smooth) spacelike "singularities" φ = φ 0 by zigzag lines, see the one in Fig. 2c.

On-shell conformal symmetry
We find that the boundary equation which change ψ(v) as an h = −1/2 primary field and T vv (v) as an energy-momentum tensor with large negative central charge Transformations from the SL(2, R) subgroup of (2.19), (2.20), have vanishing Schwarzian derivative and therefore change f in the standard way . Besides trivial translations of v they include v-dilatations due to shifts of the asymptotic coordinate t + x in Eq. (2.11) and inversion v → 1/v related to PT-reflection t + x → −(t + x). These transformations constitute the global symmetry group of our model. As a side remark, let us argue that (2.19), (2.20) is a symmetry of the gravitational degrees of freedom but not of the matter sector. To this end we introduce the field χ(u) = e −λτ (u) /ψ 0 which is T -symmetric with respect to ψ(v) and therefore satisfies where Eq. (2.13) was used in the left equality; similar expression for T vv , see Eq. (2.14). All these equations and boundary conditions can be summarized in the flat-space action In this setup (2.19), (2.20) is an apparent conformal symmetry of Φ far away from the boundary, whereas the symmetry of the matter sector is hidden in the reflection laws.

Integrable sector 3.1 General solution
One can use Eq. (2.13) to express the entire solution in terms of one arbitrary function. Indeed, introducing we find, Then U , T uu , φ, and f are given by Eqs. (3.1), (2.14), (2.7), and (2.6). We obtained general classical solution in the model with a boundary. By itself, this solution is of little practical use because the function ψ(v) has a zero at v = 0 and, possibly, another one at v =ṽ 1 > 0, see Fig. 4a. In general, the incoming flux T vv (v) in Eq. (3.2) is singular at these points. Indeed, Eq. (3.1) gives where R(v) is regular at v ≥ 0. As a consequence, T vv (v) has first-order poles at v = 0 andṽ 1 . Requiring zero residuals at these poles, we obtain two constraints Choosing multiparametric R(v) and solving the constraints, one finds an arbitrary number of smooth solutions. The physical ones satisfy In what follows we will concentrate on a large class of soliton solutions with power-law singularities. We will argue that some of them satisfy Eq. (3.3).

Soliton solutions with power-law singularities
Let us follow the Painlevé test [54] and guess the form of T vv (v) which guarantees that the general solution ψ(v) of Eq. (2.13) has power-law singularities ψ ∼ where the expansion of T vv starts from (v − v 0 ) −2 due to Eq. (2.13). Substituting Eqs. (3.4) into Eq. (2.13), we obtain an infinite algebraic system for ψ k−s , implying k 0 = 2s + 1. One concludes that s is integer or half-integer. Note that the two equations from the system (3.5) which do not determine the coefficients of ψ, constrain {T k }. For example for s = 1/2 one gets, where we expressed all ψ k−1/2 via {T k } and ψ −1/2 . For larger s, one obtains T −2 = s(s + 1) and higher-order equations listed in Table 1.
We arrived at the practical method for obtaining the soliton solutions in our model. One specifies N singularities of ψ(v): selects their integer or half-integer powers s n and complex positions v n . The function T vv (v) has second-order poles at v = v n , see Eq. (3.4). This analytic structure gives expressions, where we required T vv → 0 as v → +∞ and introduced a polynomial in the nominator of ψ(v) with M zeroesṽ m and a normalization constant C. Next, one solves equations in Table 1 at each singularity and determines T n −1 . After that ψ(v) is obtained by substituting Eqs. Example. Consider the soliton with two s = 1/2 singularities 10 . Solving the finite-energy conditions (3.8), one obtains It is straightforward to check that T vv (v) with these parameters satisfies Eqs. (3.6) at v = v 1 and v = v 2 . To make the solution real at v ∈ R, we take v 1 = a + ib and v 2 = a − ib. Then Eqs. (3.7) give, where ψ(v) was obtained by substituting Eqs. (3.7) into Eq. (2.13). One observes that the matter flux (3.10) peaks near v ∼ a, its total energy E in = 3 2 M cr 1 + a b arcctg(−a/b) is controlled by the ratio a/b , where M cr = 2λe −2φ 0 is the minimal black hole mass. Since ψ → −av as v → +∞, the solution (3.10) describes reflection of matter waves off the boundary and formation of black holes at a < 0 and a > 0, respectively, see Fig. 4a. This fact is clearly seen in Fig. 5 showing the boundary u = U (v) at different a in the finite-range coordinates (ū,v). In Fig. 5c we also plotted the space-like "singularity" φ = φ 0 and horizon u = 0 (zigzag red and solid black lines, respectively). Note that the critical solution in Fig. 5b corresponds to E in = 3 2 M cr . The simplest exact solution in Eq. (3.10) describes the incoming matter flux with a single peak. Solutions with multiple peaks can be obtained by adding singularities at v = a n ±ib n , see Fig. 6. Unfortunately, it is hard to find these solutions explicitly at large N . Besides, it is not clear whether they will satisfy the positivity condition (3.3). We will clarify these issues in the subsequent Sections. v a 1 + ib 1 Figure 6. Singularities of solitons in the complex v-plane.

Simplifying the coefficient equations
Instead of solving the equations in Table 1, one can extract T vv (v) from the general solution. Namely, substituting the solitonic ψ(v) into the first of Eqs. (3.2), we find, Then the second of Eqs.
which are, in fact, equivalent to the ones in Table 1. Indeed, after solving Eqs.

SL(2, C) symmetry
The global SL(2, C) transformations (2.21) are invertible and therefore preserve the singularity structure of the solitons. One obtains, This symmetry relates solitons with different parameters. Real solutions at v ≥ 0 transform under SL(2, R). The transformation (2.21) sends the point v = −δ/γ to infinity. If the original solution was regular at this point, its image receives asymptoticsψ → Cw + D and T vv → O(w −4 ) as w → +∞. In Eq. (2.17) we obtained the same asymptotics from physical considerations. Solutions with other asymptotics, i.e. those violating the finite-energy conditions (3.8) or Eq. (3.9), have singularities "sitting" at infinity.

Relation to the Gaudin model
In this Section we establish one-to-one correspondence between the solitons (3.7) and eigenstates of the auxiliary integrable system, the Gaudin model [47][48][49]. This will allow us to count the number of solitons and explain some of their properties. The Gaudin model [47] describes a chain of N three-dimensional spinŝ s n = {ŝ 1 n ,ŝ 2 n ,ŝ 3 n } with the standard commutation relations [ŝ α n ,ŝ β l ] = iδ nl αβγŝγ n . The model is equipped with N commuting Hamiltonianŝ where v n are complex parameters and (ŝ n ,ŝ l ) ≡ αŝ α nŝ α l is the scalar product. The eigenstates |Ψ of the model simultaneously diagonalize all Hamiltonians, T n |Ψ = T n |Ψ , where T n are complex eigenvalues.

Now, the eigenvectors satisfyT (v)|Ψ = T (v)|Ψ .
A complete set of eigenvectors and eigenvalues in the model (3.16) is provided by the algebraic Bethe Ansatz [47][48][49]. We review this method in Appendix B and list its main results below.
One fixes the representations (ŝ n ) 2 = s n (s n +1) of all spins, where s n are integers or half-integers. The simplest eigenstate |0 of the Gaudin model has all spins down, for all n , (3.18) whereŝ − n ≡ŝ 1 n − iŝ 2 n are the lowering operators. The other eigenstates are obtained by acting on |0 with rising operatorsŝ (3.19) at certain pointsṽ m which satisfy the Bethe equations, The eigenvalue ofT (v) corresponding to the state (3.19) is To sum up, one solves Eqs.
depends on the sumŝ 1 +ŝ 2 . The corresponding solutions have singularities at v = v 2 of powers |s 1 − s 2 |, |s 1 − s 2 | + 1, . . . , (s 1 + s 2 ) in accordance with the irreducible representations ofŝ 1 +ŝ 2 . For instance, consider coalescence of two s 1,2 = 1/2 singularities as v 1 → v 2 . The second-order equations (3.6) at these singularities have four solutions corresponding to four eigenstates of two s = 1/2 spins. In the limit v 1 → v 2 the spins sum up and we obtain 14 one s = 0 (non-singular) solution and three solutions with s = 1 singularity. Finally, one can obtain more general solutions with infinite number of singularities using the thermodynamic Bethe Ansatz for the Gaudin model [55].

Positivity condition
Physical solutions have real ψ(v) at real v. Thus, their singularities v n and zerosṽ m are either real or organized in complex conjugate pairs like in Fig. 6. Besides, the singularities v n may not be placed at the physical part v ≥ 0 of the real axis.
The remaining nontrivial condition is T vv (v) ≥ 0 at v ≥ 0, Eq. (3.3). This inequality is not satisfied automatically. For example, our solutions with two singularities (3.15) have negative and positive T vv (v) at v 1,2 < 0 and v 1,2 = a ± ib, respectively. In fact, any solution with all singularities placed at v < 0 is unphysical. In this case the operatorŝ(v) at real v is Hermitean, and thereforeT (v) in Eq. (3.17) has positive-definite eigenvalues T (v) ∝ −T vv (v).
In the opposite case when all singularities are organized in complex conjugate pairs v 2k−1 , v 2k = a k ± ib k with s 2k−1 = s 2k , one expects to find at least one physical solution. Indeed, consider the state |Ψ 1 (not an eigenstate) of the Gaudin model satisfying (ŝ 2k−1 +ŝ 2k )|Ψ 1 = 0 for all k. Explicit calculation shows that Ψ 1 |T (v)|Ψ 1 < 0 at real v. On the other hand, the variational principle implies that for any N real points w n there exists an eigenstate |Ψ minimizing all Ψ|T (w n )|Ψ . The respective eigenvalue T (v) is negative at all v = w n suggesting that T vv (v) ∝ −T (v) is positive at the entire real axis.
Let us explicitly select the above physical solution at b k → 0. In this case T vv (v) falls into a collection of peaks at v ∼ a k near the singularities v 2k−1 , v 2k .
At |v − a k | b k and yet, far away from other singularities, the operator (3.17) takes the formT (v) ≈ (ŝ 2k−1 +ŝ 2k ) 2 /(v − a k ) 2 . Its eigenvalue T (v) ∝ −T vv (v) is positive-definite unless the eigenstate satisfies (ŝ 2k−1 +ŝ 2k )|Ψ = 0. Thus, in the limit b k → 0 the physical eigenstate coincides with the state |Ψ 1 introduced above. The respective energy flux T vv (v) is the sum of two-spin terms (3.15), One expects that this solution remains physical at finite b k .
Example. In general case the positivity condition bounds parameters of the solutions. Consider e.g. the soliton with three s = 1 singularities at v 1,2 = a ± ib, v 3 < 0, see Fig. 8a. Solving Eqs. (3.8), (3.9), one obtains, . (3.23) The second (negative) term in this expression represents contribution of the singularity v 3 < 0. It can be compensated by the first term if the singularities v 1 and v 2 are close enough to v 3 . Namely, the function , see the gray region in Fig. 8b. The solutions with these parameters involve one peak of the incoming flux, just like the solutions (3.15).  [46], cf. [57,58]. At energies somewhat below critical the boundary has long almost null part ("plateau"), see Fig. 9a. The energy flux reflected from this part is strongly amplified by the Lorentz factor of the boundary and forms a high and narrow peak in T uu (u), see Fig. 9b. We will argue below that in the critical limit the peak tends to a δfunction (shock-wave) with energy equal to the minimal black hole mass M cr . In the overcritical solutions the shock-wave is swallowed by the black hole. Besides, we will see in the next Section that the structure of the peak is highly sensitive to the initial data. This feature impedes global integrability of the model.  Let us find the boundary U (v) in the "plateau" region where v is large and T vv (v) is small. In this case Eq. (2.13) can be solved perturbatively by representing ψ = 1 + ψ (1) + ψ (2) + . . . , where ψ (k) ∝ (T vv ) k . Using ψ ≈ 1 in the r.h.s. of Eq. (2.13), we obtain, where the function g(v) is introduced in Eq. (2.8) and g ∞ is its value at v → +∞. Note that the linear asymptotics Cv 1 of the solution appears at first order of expansion in Eq. (4.1) because in the near-critical regime ∂ v ψ ≈ C is small at large v.
In what follows we will regard C as a parameter of the expansion. Using ψ ≈ 1 + ψ (1) in the r.h.s. of Eq. (2.13), we get The higher-order corrections ψ (n) are obtained in similar way. Now, we compute the reflected energy flux T uu (u) and the boundary function U (v) using Eqs. (2.14) and (2.12), We kept one and two orders of the expansion in Eqs. (4.2) and (4.3), respectively. Note that the leading (first) term in U (v) is constant; this behavior corresponds to the "plateau" in Fig. 9a. At the same time, the reflected flux (4.2) has a peak at large v corresponding to ∂ v g ∼ Ce −2φ 0 . This peak is narrow in terms of slowly-changing u = U (v) in agreement with Fig. 9b. Using the soliton asymptotics T vv ∝ v −4 and ∂ v g ∝ v −3 , one finds that the peak in Eq. (4.2) occurs at v ∝ C −1/3 , and its width ∆v is of the same order. The respective value of U (v) is approximately given by the first term in Eq. (4.3), while the peak width ∆U ∝ C 2/3 U is controlled by the second-order terms. In the critical limit C → 0 the peak of T uu (u) is infinitely high and narrow.
Calculating the total energy within the shock-wave at C → 0, we obtain, where Eqs. (4.2), (4.3) were used. The value of E peak coincides with the minimal black hole mass M cr implying that the peak of T uu (u) tends to a δ-function in the critical limit.

Shock-wave instability
Since our model is equipped with the general solution, one may think that it is integrable, i.e. has a complete set of conserved quantities {I k } smoothly foliating the phase space. In the in-sector these quantities are arbitrary functionals I k [f in ] of conserved f in (v), cf. [59]. Then, I k can be computed at arbitrary space-like line: to this end one evolves the classical fields from this line to J − , extracts the incoming wave 15  Let us argue, however, that {I k } cannot be smoothly defined in the near-critical regime because the map f in → f out in this case is essentially singular. To simplify the argument, we consider solutions with the modulated flux at large v, where C is the small parameter of the near-critical expansion.
If ω is small as well, the asymptotics of T vv is almost power-law, like in the ordinary solitons. However, the shock-wave part of the reflected flux represents squeezed and amplified tail of T vv at v ∼ C −1/3 , see Fig. 9. It should be essentially modulated. For simplicity, let us characterize the outgoing wave packet with a single quantity where we used the flat coordinates (2.11) in the definition of I 3 , then separated the shock-wave part ∆I 3 of the integral at t − x ≡ − log(−λu)/λ log C from the (C, ω)-independent contribution at smaller t − x. In the second line we substituted the shock-wave profile (4.2), (4.3) and extended the integration range to v ≥ 0. Now, one substitutes the asymptotics (4.4) into Eq. (4.5) and finds that ∆I 3 (C, ω) is quasiperiodic. Indeed, change of the integration variable v → ve 2πn/ω with integer n gives relation 16 ∆I 3 (e 6πn/ω C, ω) = e −2πn/ω ∆I 3 (C, ω). Thus, ∆I 3 We see that ∆I 3 has an essential singularity at ω = C = 0. Indeed, taking the limit C → 0 along the paths ω log C = const, one obtains ∆I 3 → −∞, 0, or +∞, see Fig. 10. Thus, any value of ∆I 3 can be obtained by adjusting the limiting path.
The above property ascertains dynamical chaos in the critical limit of our model. Indeed, infinitesimally small changes (4.4) of the initial data at small C produce outgoing fluxes with essentially different values of I 3 . This prevents one from characterizing the critical evolution with a set of smooth conserved quantities I k . Indeed, all functionals I k [f in ], being smooth in the in-sector, are not sensitive to ω at small values of latter. Thus, they fail to describe essentially different out-states f out (u) at different ω. From a more general perspective, one can introduce the integrals which are smooth either in the in-sector or in the out-sector, but not in both.

Discussion
In this paper we considered two-dimensional CGHS model with a regulating dynamical boundary [33,34]. This model is weakly coupled and causally similar to the spherically-symmetric gravity in many dimensions. We demonstrated that classical field equations in this model are exactly solvable. We constructed their general solution and studied in detail a large subset of soliton solutions with transparent properties. We illustrated the results with many explicit examples hoping that this model will serve as a practical playground for black hole physics.
In the critical regime i.e. at the verge of black hole formation, our model displays dynamical instabilities specific to chaotic systems. This property is similar to the near-horizon chaos suggested recently in the context of AdS/CFT correspondence [23][24][25][26][27][28]. We argued that it hinders global integrability of the model. We see several applications of our results. First, exact solvability may extend to one-loop semiclassical level if one adds a reflective boundary to the RST model [5]. This approach, if successful, will produce analytic solutions describing black hole formation and evaporation. The singularities of such solutions should be either covered by the boundary or hidden behind the space-like line φ = φ 0 , see Fig. 4b. Then a complete Penrose diagram for the evaporation process may be obtained, cf. [39,41,42,46].
Second, in the alternative approach one directly adds one-loop corrections to the classical equations of our model with a boundary and integrates the resulting system numerically, cf. [60,61]. By the same reasons as above, the respective solutions should completely describe the process of black hole evaporation.
Third and finally, the model of this paper is ideal for applying the semiclassical method of [62,63] which relates calculation of the exponentially suppressed S-matrix elements to certain complex classical solutions. The results of such calculations may be used to test unitarity of the gravitational S-matrix [63]. this project. We are grateful toÉcole Polytechnique Fédérale de Lausanne for hospitality during our visits. This work was supported by the grant RSCF 14-22-00161.

A Field equations and boundary conditions A.1 Derivation
Field equations in the bulk are obtained by varying the action (2.1) with respect to g µν , φ, and f , and ignoring the boundary terms, 3) The first line here relates the energy-momentum tensors of φ and f , −T µν . The second line implies, in addition, that the rescaled metric e −2φ g µν is flat.
To find the boundary conditions at the line φ = φ 0 , we keep the boundary terms in the variation of the action. For a start, let us consider variations preserving the coordinate position of the boundary φ = φ 0 . We take δφ = 0 along this line and fix the direction of its outer normal, δn µ ∝ n µ . The integration domains in Eq.

A.2 Solution in the conformal gauge
Let us review the general solution [4] of the bulk equations (A.1)-(A.3), see [1,2] for details.
In the light-cone frame (2.4) Eq. (A.3) takes the form ∂ u ∂ v f = 0, its solution is given by Eq. (2.5). Combining Eq. (A.2) with the trace of Eq. (A.1) and substituting R = 8e −2ρ ∂ u ∂ v ρ, we obtain, where the residual coordinate freedom 17 was fixed in the last equation. After that Eqs. (A.1), namely, In this expression M − , u 0 , and v 0 are integration constants; functions g(v) and h(u) were introduced in Eq. (2.8). We fix u 0 = v 0 = 0 by shifting u and v. After that M − represents the mass of white hole in the infinite past [1,2]. Indeed, the past time infinity i − in Fig. 1b  It is worth noting that the patch u ∈ (−∞, 0) and v ∈ (0, +∞) covers all spacetime accessible to the outside observer. Indeed, we already mentioned that the time infinities i − and i + are reached in the limits u → −∞ and v → +∞ at finite values of the dilaton field φ. By Eq. (2.7), the product uv remains finite in these limits implying v → +0 as u → −∞ (i − ) and u → −0 as v → +∞ (i + ), see Fig. 3.
We proceed by deriving equation of motion for the boundary u = U (v) satisfying φ(U (v), v) = φ 0 . Taking the derivative of Eq. (2.7) along this line, we find, where U ≡ dU/dv > 0 because the boundary is time-like. The other two equations come from the boundary conditions (2.2). Introducing the unit outer normal