Floquet Scalar Dynamics in Global AdS

We study periodically driven scalar fields and the resulting geometries with global AdS asymptotics. These solutions describe the strongly coupled dynamics of dual finite-size quantum systems under a periodic driving which we interpret as Floquet condensates. They span a continuous two-parameter space that extends the linearized solutions on AdS. We map the regions of stability in the solution space. In a significant portion of the unstable subspace, two very different endpoints are reached depending upon the sign of the perturbation. Collapse into a black hole occurs for one sign. For the opposite sign instead one attains a regular solution with periodic modulation. We also construct quenches where the driving frequency and amplitude are continuously varied. Quasistatic quenches can interpolate between pure AdS and sourced solutions with time periodic vev. By suitably choosing the quasistatic path one can obtain boson stars dual to Floquet condensates at zero driving field. We characterize the adiabaticity of the quenching processes. Besides, we speculate on the possible connections of this framework with time crystals.


Introduction and main results
The physics of periodically driven many-body systems is a fascinating chapter in the study of out-of-equilibrium dynamics [1,2]. Their behavior can differ substantially from that of their static counterparts. Floquet systems have been investigated at criticality with the tools of conformal field theory (CFT) (see [3] and references therein). The naive expectation that they can only increase indefinitely their energy because of the driving can be avoided in some regions of the driving parameter space. This was confirmed by explicit examples that attain a steady state at finite temperature [3,4]. Even though the actual saturation mechanism is not completely understood, the fact that periodic driving can radically alter the stability of equilibrium points is well known from nonlinear dynamical systems, the Kapitza pendulum being a prominent example [5].
The AdS/CFT correspondence maps equilibration processes onto the dynamical evolution of a dual gravitational system. The holographic dictionary defines the sources of the field theory in terms of asymptotic modes of the bulk fields. It is then possible to study a driven holographic system coupled to an external source. In the cases of interest here, the source corresponds to the leading asymptotic mode of a bulk scalar field in an asymptotically global AdS spacetime.
Periodic driving has been scarcely studied in the context of AdS/CFT. The analysis performed in [6] comes closest in spirit to the present study (see also [7][8][9][10][11], all in the probe approximation). The main difference is that their initial states are already mixed, while here we consider the driving of pure states. Their geometrical ansatz is embedded in the Poincaré patch, with an arbitrarily small planar horizon to cap the IR. They find a monotonous growth of the black hole horizon as a general outcome of the driving. Such growth leads to an infinite temperature final state. The result is, on the one side, compatible with the already mentioned natural expectation and, on the other side, it comes out as a consequence of the ingoing boundary conditions at the horizon. The only remaining freedom may be in the rate of growth of the energy density and entropy.
The authors in [6] find three different phases depending upon the range of values acquired by the dimensionless combination χ = φ b ω b , namely the product of the amplitude and the frequency of the driving at the boundary. Small values of χ correspond to a dissipation dominated regime, where the horizon absorbs all the energy entering the system from the boundary. Here the scalar field synchronizes with the driving but does not change its amplitude. Raising χ they find two other phases in which the scalar field stops being "transparent", and it may absorb part of the energy input and increase its amplitude.
The picture of [6,10] should be compared with the results we find here. The fact that our driving does not necessarily imply decoherence (i.e. a collapse into a black hole) is by no means trivial and it only occurs in certain regions of parameter space. Charting this region, characterizing its features, and studying its stability is the main target in the present paper.
We have examined both complex and real massless scalar fields. While a priori their dynamics is quite different, remarkably, the final results of the real and complex cases exhibit strong similarities, both in the structure of the phase space and in the order of magnitude of specific numerical values like stability thresholds.
Finding the manifold of driven and periodic solutions is the first step. Since such solutions could in principle be all unstable, a thorough inspection of their stability is in order both at the linear and nonlinear level. Linearly stable solutions proved to exist and be robust against rather large fluctuations. For unstable solutions, we follow numerically the evolution in the search for the final endpoint at large times. The picture that emerges is rather interesting: in significant portions of the unstable subspace, the solutions lie at the boundary of two radically different long time behaviors. On one side, fluctuations drive the geometry to a collapse into a black hole, the dual theory loses coherence and thermalizes as expected. On the other side, the geometry is regular forever. It however exhibits a periodic modulation that has the form of a relaxation oscillation 1 where the solution periodically bounces back and stays for a long time close to the initial unstable solution. From the dual field theory point of view, the quantum state remains pure but exhibits an emergent pulsating modulation.
The framework at hand allows us to address the important question of the adiabatic preparation of Floquet condensates [12][13][14]. By slowly varying the amplitude and/or frequency of the driving, one can study if the system keeps up with the source and follows a trajectory on the manifold of periodic solutions. This is naively expected in the limit of slow quenches, as a natural dual counterpart of the quantum adiabatic theorem. However, a careful analysis reveals subtleties whenever the driving frequency approaches that of a normal mode of AdS. Similarly, we study the quenching of the system from AdS through a cyclic protocol starting and ending with a vanishing driving. The possibility that such cycle ends on a periodic solution with vanishing source (but not just AdS) could perhaps relate to the open question of spontaneous breaking of continuous time translations [15,16].
The structure of the paper is as follows. In Section 2 we study the case of the driven complex field. We characterize in detail the space of stationary solutions and focus on their linear and nonlinear stability properties. In Section 3 we analyze the response of the system upon a continuously varying driving. This unravels interesting structures like the presence of critical values of the quench parameters separating classes of qualitatively different time evolutions (ending or not in a black hole, for instance). In Section 4 we show that the unstable solutions are, in fact, attractor solutions in the sense of Type I gravitational phase transitions studied in [17][18][19]. In Section 5 we describe two alternative quenching protocols to construct a boson star starting from AdS. One is quasistatic, while the other involves an unstable solution as an intermediate step. Section 6 is devoted to the late-time evolution of the collapsing black hole solutions and establishes contact with the results of the analysis done in [6]. In particular we are able to pin down the three regimes regimes found in that paper. In Section 7 we comment the analysis for the real scalar field, stressing the fact that it constitutes a nontrivial extension both at the conceptual and technical level. We summarize in Section 8 accounting for and interpreting the results as well as indicating research directions for the future. Some technical material is collected in the appendices.

Periodically driven complex scalar field
Our case study involves the simplest possible setup, namely a complex scalar field in global AdS 4 , with Λ = −3/l 2 for AdS 4 . We will set κ 2 = 8πG = 1, l = 1. The action is invariant under global U (1) transformations φ → e iα φ. Our ansatz for the metric is where x ∈ [0, π/2) is the radial coordinate. We will examine the space of solutions adapted to the following time-periodic ansatz with ρ(x) a real scalar function. The space of solutions contains just three functions f, δ and ρ; they are determined by two parameters, namely, the frequency ω b and a boundary value for ρ.
Introducing the ansatz (2.3) in the Einstein-Klein Gordon equations of motion, gives a set of ω-dependent static equations for ρ(x) which can be solved numerically by standard methods (see (A.11)). Each solution of this system gives rise to a static geometry where the scalar field rotates in the complex plane with angular velocity ω b while keeping constant its modulus ρ(x). The boundary values of ρ(x) are respectively ρ o = ρ(0) at the origin and ρ b at the boundary; ρ b is defined through the asymptotic Taylor expansion with ∆ = 3/2 + 9/4 + m 2 being the conformal dimension of the dual scalar operator O. For non-vanishing ρ b , the dual state is interpreted as being driven by a time-periodic source ρ b e iω b t . We abbreviate these sourced periodic solutions as SPS. In this paper we analyze the case of a massless scalar, m 2 = 0, in detail. We shall also comment on partial results obtained in the case m 2 = −2, where everything seems to follow the same pattern so far. In the massless case, the unsourced solutions with ρ b = 0 are termed boson stars -BS-in (at least part of) the literature [20,21] and we shall adhere in what follows to this name. In making this plot, we have adopted the phase convention that makes ρ o a positive number (see Fig. 3). Hence the sign of ρ b is nothing but the relative sgn(ρ o /ρ b ), which is correlated with the number of nodes of φ(x). Figure 2 will be an essential tool to understand the results of our analysis. The surface represents the complete set of static geometries corresponding to SPS's in the three dimensional space spanned by (ω b , ρ o , ρ b ) (see Fig. 1 for a visualization of these three parameters). The SPS surface cuts at ω b = 0 on a line of configurations where the scalar takes a constant radial profile, in particular ρ o = ρ b . Boson stars lie at the intersection of the SPS surface with the zero source plane ρ b = 0. The boson star curves start from the bottom plane ρ o = 0 at the values ω n = 3 + 2n, given by the spectrum of normal frequencies of the massless scalar in AdS 4 . As the value of ρ o is increased, the nonlinear dressing of the linearized solution shifts the value of the frequency and gives rise to the wiggly curves depicted in Fig. 2. The lower portion of this curve that bends towards smaller values of ω b was already constructed in [21]. The different sectors indicated with I,II,III,... correspond to SPS whose radial profiles ρ(x) have 0, 1, 2... nodes respectively. As the number of nodes increases, the energy density tends to concentrate deeper in the bulk. The reader should understant that this profile is rotating in the transverse plane. The relative sign ρ b and ρ o correlates with the even or odd number of nodes.
In the case of a driving by a relevant scalar field with m 2 = −2, the surface remains qualitatively the same. The unsourced solutions analogous to the boson stars branch instead from the spectrum of linearized fluctuations which is now ω n = ∆ + 2n with ∆ = 1, 2 for the alternative/standard quantization of the scalar field.
Given the external driving force, the existence of regular solutions constitutes by itself a nontrivial result, the natural expectation being that the system would get increasingly excited. The total mass should grow monotonically and, eventually, lead to a collapse, signalling thermalization in the dual field theory side. Before closing this section let us comment why this may be avoided here. The rate at which the total mass changes with time is controlled by the diffeomorphism Ward identity For the case at hand, this identity boils down to the following equatioṅ In the undriven caseφ b (t) = 0, the r.h.s. is identically zero, and the system is isolated. Foṙ φ b (t) = 0, the magnitude and sign of the product on the right hand side is unforeseeable. The product of the two factors can be either positive or negative. This implies that, for a variable source, energy can flow both in and out. The source function φ b (t) is known. In contrast, the 1-point function O(t) is a teleological quantity (in the radial direction), like event horizons are (in time). It can only be extracted from the boundary behaviour of the solution once this has been computed down to the origin and regularity has been imposed. In the particular case of a SPS, the vev oscillates harmonically in phase with the source. This particular case yields an exactly vanishing value for the right hand side. However, as soon as the SPS is perturbed,ṁ will start to fluctuate around zero. The average of this fluctuation will signal whether the mass starts to build up and the system eventually collapses or, else, if the net balance is zero, the perturbed solution stays regular in the future.
An important remark is that there is a chance to find a regular solution only if the quantum system we drive is in a pure state. If there was a horizon, no matter how small, then part of the injected energy would fall behind it, and never reach back to the boundary. Unavoidably, the mass would then grow monotonically and the black hole horizon end up reaching the boundary, i.e. reaching the dual geometry to an infinite temperature state.
In summary, a thorough study of the stability properties of the SPS is compulsory. We will do it first in a linearized approximation and then in a fully nonlinear setup.

Linear stability
In this section we shall establish and chart the regions of stability within the set of SPS's given in Fig. 2 by performing a linearized fluctuation analysis. We find it convenient to move along level curves with constant source amplitude ρ b . These curves are depicted in Fig. 4.
Notice that from the point of view of an observer at the boundary, namely for each pair (ρ b , ω b ), there can be more than one SPS. They have different bulk profiles and, in particular, they reach the origin at different values of ρ o . For instance, in Fig. 4 a duplicity has been highlighted in the case (ρ b = 0.09, ω b = 2) by the red dots on the left. In the case of multiple solutions corresponding to the same boundary data (ρ b , ω b ), it is natural to suspect that one of them is stable since it represents the "ground state" of the sector.
Along the lines of constant ρ b , one finds extremal values of the frequency where two solutions corresponding to the same ω b become degenerate. At these points the spectrum of linearized Turning points of physical quantities are natural locii for the onset of linear instability. The paradigmatic example is the Chandrasekhar mass of a white dwarf. In our analysis we also find some maxima of the mass along the level curves which are accompanied by a squared eigenfrequency transiting from positive to negative values, see for instance Fig. 39. Nevertheless, looking at the extrema of the mass (or ω b ) does not yield exhaustive information about the stability. There are cases where complex eigenfrequencies arise because two real eigenmodes merge. Further information about the mode structure and behavior can be found in Appendix C.
Charting the complete stability region involves a numerical scan of the spectrum of linearized perturbations of SPS's across the space of solutions. The shaded region in Fig. 4 is our best approximation as to where the region of linearly stable solutions extends. Notice that parts of the edge of the stability region are given by the boson star solutions. Here the passage from gray (stable) to white (unstable) occurs precisely because the lowest eigenmode squared turns negative. This observation implies that the spectrum of linearized perturbations around a boson star will always contain a zero mode. In Appendix C we prove that the zero mode generates the boson star line and we provide more information on the building and features of the stability region.
Also notice that the white wedges emerge from integer values of ω b . Such values correspond to special eigenfrequencies of the linearized scalar field problem. Indeed, over global AdS 4 , the scalar wave equation has a general spectrum of regular solutions given by Normal modes come in two families. Solutions with ω b = 3+2k, k = 0, 1, 2... are normalizable and have vanishing source, while solutions with ω b = 2k, k = 1, 2, 3... are non-normalizable and have nonzero source, but vanishing vev. 2 The fact that white wedges descend to the vicinity of these linearized solutions supports the picture of linear instability as a resonance phenomenon. With this remark in mind we would expect also an instability wedge to come down to ω b = 2 but actually we do not find it.

Nonlinear stability
In a nonlinear theory the results of a linear stability analysis are of limited range. A perturbed unstable solution will soon start departing largely from its original, unperturbed 2 This is more easy to see by rewriting the normal modes for AdS d+1 in terms of generalised Legendre polynomials e k (x) = cos ∆ xP state. Therefore, we need to follow it to see what the end result of its evolution can possibly be. In finite dimensional nonlinear systems there are two possibilities: another stable equilibrium point or a limiting cycle. We have built a numerical evolution code which is a minimal adaptation of the ones employed in the study of gravitational collapse in global AdS implementing a standard RK4 finite difference scheme. We have checked that, when initialized at t = 0 with a linearly stable SPS, the output yields consistently the expected periodic geometry over as long as we have let the computer go. Notice that, since the profile never gets particularly spiky, rather low resolutions of 2 10 + 1 points are sufficient. This allows for a considerable increase in computational speed.
The next numerical experiment is to initialize the code with a linearly unstable SPS slightly perturbed with the single unstable normal mode. In other words, if φ 0 (x) is an SPS, and χ 1 (x) its unstable mode with a purely imaginary eigenfrecuency λ 2 1 < 0, then at t = 0 we will insert φ 0 (x) + 1 χ 1 (x) with 1 ∼ O(10 −4 ). In general terms, one could reasonably expect collapse to a thermal phase as the end result of time evolutions starting from unstable initial conditions. Nevertheless, as we shall see now, even in this case, the system provides regular counterexamples.

Re[ϕ(t,0)]
Re[ϕ(t,π/2)] 20 40 60 80 100 t -0.5 0.5 Figure 6: In this plot we appreciate that the driven scalar field keeps in phase with the driving all along its profile. The tip at the origin and at the boundary rotate at the same pace. There is, a priori, no reason why this should be so and is not an artifact or truncation of the simulation. In fact one can find (unstable) solutions where this is not the case.
In order to keep the analysis as systematic as possible, we continue working with the two solutions marked in red on the left hand side of Fig. 4. They both represent SPS's with ω b = 2 and ρ b = 0.09. The lower one lies in the region of linear stability, and we have checked that it also supports fairly "strong kicks", 1 ∼ O(1), leading to oscillations about the initial solution.
In contrast, the upper one contains in its spectrum of linearized normal modes one unstable eigenmode, χ 1 (x), with a purely imaginary frequency λ 2 1 = −0.2861478. This fundamental eigenmode has no nodes and its shape resembles the SPS itself. With 1 ∼ 10 −4 , the time evolution when t → ∞ differs dramatically depending upon the sign chosen for 1 . With positive sign (so that 1 χ 1 (0) has the same sign as φ 0 (0)), the evolution collapses promptly to a black hole geometry. This can be seen in the top line of plots in Fig. 7. The continuous driving reflects itself in the rising of the total mass, while the horizon radius increases monotonically as does the absolute value of the vev itself. The system approaches the expected infinite temperature state. after being perturbed with the first unstable mode φ(0, x) + 1 χ 1 (x). The upper (lower) three plots correspond to 1 = −10 −4 (+10 −4 ). The magnitudes shown are the minimum of the metric function min x∈[0,π/2) (f (t, x)), the energy density m, and the vev | O |. A relaxation oscillation with a pronounced peak and a long lived plateau is apparent on the lower plots.
With the opposite sign, 1 ∼ −10 −4 , the evolution reaches a radically different endpoint. Instead of a collapse, the dynamics stays regular for all times. In Fig. 6 we plot the oscillations of the real part of the scalar field at the origin and at the boundary. Both of them proceed in phase with one another, signalling that the profile of the scalar field stays in a plane while rotating. Moreover, and this is new, the oscillation of the scalar at the origin acquires a violent and sudden modulation in amplitude.
Since we are working with a fully backreacted solution, the modulation of the driven oscillation affects many other physical quantities. Fig. 7 highlights the impact of the modulation on three observables: the minimum of the metric function min x [f (t, x)], the absolute value of the vev O (t) and, finally, the energy density M . Note that min x [f (t, x)] reveals the formation of an apparent horizon whenever it drops towards zero.
The periodicity of the modulation should not obscure the fact that it is highly anharmonic. The oscillations in modulation stay on long lived plateaux followed by sudden pulsations or beats. We will refer to these regular solutions as sourced modulated solutions, SMS. The period T is orders of magnitude away from the one due to the driving frequency (2π/ω b ) and, actually, it is (weakly) dependent on the strength of the perturbation. Hence, in the strict sense, the solution is not a limiting cycle like the one found in nonlinear systems such as the van der Pol oscillator. In such systems the asymptotic dynamics eventually loses memory of the initial state whereas, in our case, it bears resemblance to relaxation oscillations [22]. The closest mechanical analogy for a SMS would be that of an inverted pendulum slightly kicked out of its vertical unstable position, where it will return and stay for a long time until the next sudden turn.
In fact, and very remarkably, these regular solutions also occur with vanishing source. Indeed, unstable boson stars on the black curve slightly above the yellow dots in Fig. 4 also exhibit this modulated dynamics for one sign of the most relevant perturbation, and collapse to a black hole for the other. 3 Such behavior has been observed earlier, both in asymptotically flat [23] and AdS spacetime [24,25]. Also, SMS profiles like the one in Fig. 6 bear some resemblance with the long time modulation of oscillating solutions that appear in confining theories after a global quench that injects energy below the mass gap threshold [26,27]. At first sight, the modulations we find look substantially more anharmonic and strongly peaked. A Fourier analysis should be carried out in order to reach a definite conclusion. Another interesting proposal could relate the SMS's to nonlinearly dressed multi-oscillator solutions [28]. These still have to be constructed in the driven situation, but most likely this is feasible.
The unstable SPS's with low enough amplitude (indicatively of the order of ρ o < 0.3) undergo an exotic evolution. In fact, in these cases the pulsating modulation appears for both signs of the relevant perturbation 1 . With the "potentially lethal" sign (the one that usually would drive the SPS towards a black hole) we now observe a pulsating modulation where the value of ρ o increases instead of decreasing, yet reaching a maximum value and bouncing back again. The mass and the vev also increase and decrease back and the solution stays regular forever. not only perturb the absolute value of the scalar field in the bulk, but also its phase. More precisely, during the pulsation, the two values φ o (t) at the origin and φ b (t) at the boundary start de-phasing in one sense or the other depending upon the sign of 1 , and the scalar profile becomes increasingly helical. After the pulsation has ended they again re-phase and the profile recovers its planar shape.
The study of nonlinear stability is always a long and very much resource dependent task. What we have done so far is to characterize the endpoint of linearly unstable SPS's and we have found two possibilities, a dynamical black hole and a SMS. Concerning this last one, we have indeed checked for their robustness by adding a perturbation to them. Perturbed SMS's develop wiggling plateaux but remain regular all along the simulation, suggesting that their normal mode spectrum contains only real eigenfrequencies. The very late time behaviour of the solution is another delicate point of the nonlinear stability analysis. As far as our codes have run, up to t = O(10 4 ), we have not found the slightest evidence of a nonlinear instability setting in at late times. After the discovery of the AdS instability [29] this is an issue one should be concerned about. However, the essential ingredient found in that context, namely full resonance, is very unlikely to be present in the spectrum of linearized perturbations of a SMS. In the absence of any potential mechanism for destabilization at long times, any runtime is disputable in what concerns any claim for stability.
Finally, let us stress again that all the analysis has been performed in region I of the phase space shown in Fig. 4. It would be interesting to do an accurate scan to see if further exotic evolutions can be obtained along the wedge of unstable solutions. In particular, we have not analyzed the fate of unstable solutions in regions II, III... This remains for a later investigation.

Quenches with periodic driving
The sourced periodic solutions we have found are eternal, extending from t = −∞ to t = +∞. Therefore, it makes sense to try to study if they can be constructed by means of a slowly growing source starting from AdS. In this section we perform this investigation by considering quenches in a generalized sense. The word quench usually refers to a change in some coupling constant, which can be either sudden or slow. In the present context, it will refer to a certain process that interpolates between two different periodically driven Hamiltonians. In short, we shall explore how the system responds when the boundary data that determine the periodic driving of the dual QFT become, themselves, functions of time (ρ b (t), ω b (t)).

Quasistatic quenches
We analyze the system response to very slow changes of the driving parameters, whose variation occur on a typical time scale β ω −1 b . In the context of static quantum systems the adiabatic theorem states that, for sufficiently slow quenches, the ground state of the system follows the change in the Hamiltonian. Even if a non-vanishing transition amplitude to an excited state is generated, there is a well defined way to bring it to arbitrarily small values (by making β large enough). Here we find that a similar phenomenon occurs, albeit the underlying unquenched Hamiltonian is time-periodic. The system follows the slow modulation of the parameters of the periodic driving, moving from one ground state to another.
As we already know, SPS exist on a codimension-one submanifold in (ω b , ρ o , ρ b ) space. In this subsection, we demonstrate that linearly stable SPS's can be reached from the global AdS 4 vacuum by means of a sufficiently slow quench. Conversely, we also show that these linearly stable SPS determine which quench processes cannot be regarded as adiabatic. We employ the following ansatz for the scalar field source i.e., the source starts being zero at t = 0, and reaches its final amplitude after a time span β. Afterward, it oscillates harmonically. During the whole process, the frequency ω b is kept constant. We refer to the regime taking place at t < β as the build-up phase, while the driving phase corresponds to t ≥ β. The quench profile is plotted in Fig. 9, where ρ b (t) ≡ |φ(t, π/2)|.

Quasistatic quench to a SPS
Imagine that the final values of and ω b correspond to a SPS in the region of stability, as given in Fig. 4. For concreteness, let us focus on region I, and consider the lower red dot depicted at ρ b = 0.09 and ω b = 2. In this case, if β is sufficiently large, the state obtained after the build-up process is, to an excellent approximation, the SPS corresponding to the final harmonic driving. Furthermore, during the driving phase the numerical solution also remains stationary (up to the largest times we have simulated and with high accuracy).
Let us choose β = 2500. We have determined numerically that, after the quench (i.e., for t > β), the energy density of the geometry corresponds to m = 0.2121, which agrees with the energy density of the SPS we expect to land on up to the fifth significant digit. A more detailed check is provided in Fig. 10, where we compare the fields ρ and f during the driving phase (solid curves) and their profiles in the corresponding SPS (dashed curves).
Having demonstrated that the linearly stable SPS can be reached from the vacuum by a sufficiently slow quench process, it remains to analyze the adiabaticity of the whole procedure. Specifically, if the system's response to the time-dependent source is perfectly adiabatic, we expect that where, with no loss of generality, Ψ denotes any field of the geometry. Ψ SP S is the value of this field in the SPS at the given instantaneous ρ b and ω b . The equality (3.2) thus entails that the dynamics does not depend explicitly on time, but only implicitly through the instantaneous value of the source. In other words, in the quasistatic large β limit, one can draw the evolution as a path on the surface of SPS's (in this case, the vertical dashed line displayed in Fig. 4).
Note that these observations should also hold for one-point functions such as m and O : if the response of the system to the quasistatic quench process is perfectly adiabatic, it must be the case that and similarly for the scalar vev. From left to right, β = 500, 1000, 1500, 2500. Right: m(t/β) for the quenches of the right figure.
The curves clearly collapse into a universal profile (dotted black). We also plot m SP S (ρ b (t/β), ω b ) in (solid) yellow. The agreement between both curves implies that the system responds adiabatically.
As we illustrate in Fig. 11, these expectations are fulfilled. Let us look first at Fig. 11a. There, we compare the time evolution of m for four quench processes, with β = 500, 1000, 1500, 2500.
Since for the quench profile (3.1) ρ b (t) only depends on the dimensionless ratio t/β, relation (3.3) implies that for adiabatic response the energy density can only be a function of t/β, but not of t and β separately. Therefore, when plotted in terms of t/β, the four instances of the time evolution of m shown in Fig. 11a must collapse to a universal curve. This is indeed what happens, as we illustrate in Fig. 11b (dotted black). Finally, in Fig. 11b we also depict m SP S (ρ b (t/β), ω b ) (solid yellow). The evident agreement between both curves shows that relation (3.3) is satisfied.

Quasistatic quench to a black hole
Even in the β → ∞ limit, there exists a bound to the amplitude or the frequency of the driving that one can reach. Consider a quench at constant ω b from global AdS 4 to a final source amplitude, , such that there is no SPS associated to ω b and . As the source amplitude builds up, and in the quasistatic limit, we expect that the system evolves through a succession of SPS's, at most up to the time t * when the SPS associated to ρ * b ≡ ρ b (t * ) ceases to exist. 4,5 Past this point, the system exits the surface of SPS's and adiabatic evolution cannot proceed further, no matter how slowly the source amplitude increases afterwards.
Let us analyze these questions in a specific example. We consider the quench described by equation (3.1), again at driving frequency ω b = 2, but increasing now the final source amplitude to = 0.1. There does not exist a SPS associated to this particular driving: the highest driving amplitude compatible with the existence of a SPS with ω b = 2 is ρ * b ∈ [0.099200, 0.099225]. According to the quench profile (3.1), this value is reached at a time t * ∈ [0.73482, 0.73581]β ≡ ηβ. Here η is critical value of the scaling variable t/β after which we expect adiabatic evolution to break down.
In Fig. 12a, we plot the time evolution of m for quench profiles with time spans β = 500, 1000, 1500, 2000, 2500, 3000. We clearly see that gravitational collapse takes place at times t < β. In order to understand adiabaticity and its breakdown, in Fig. 12b we plot m(t/β). Two facts are manifest.
The first one is that, for t/β < η (signaled by the vertical dashed red line), all the energy density curves merge into a universal profile. This observation immediately leads to the conclusion that m(t) is set solely by the instantaneous value of ρ b (t/β) and, as a consequence, any nontrivial dynamics that depends explicitly on t is highly suppressed in the quasistatic limit. The system evolves adiabatically, as expected.
The second one is that, for t/β > η, the system undergoes gravitational collapse (as can Left: time evolution of the energy density for the quench processes described in the main text. From left to right, the quench time span corresponds to β = 500, 1000, 1500, 2000, 2500, 3000. Right: time evolution of the energy density, now plotted against the scaling variable t/β. We clearly observe that, prior to the adiabaticity breaking, m(t/β) approaches a limiting curve as β increases. The vertical line marks the time t * /β = η = 0.73482.
be seen in the unbounded increase of its energy density). Given that, for t/β > η, we have , there is no SPS associated to the driving), we conclude that adiabatic evolution breaks immediately after crossing the ρ b (t) = ρ * b threshold.

Non-quasistatic quenches
Let us address now the non-quasistatic regime, namely that of quenches with finite time span β < ∞. From the QFT side the standard lore is, according to the adiabatic theorem, that the system will not keep up with the variation of the source, and will generically transit to an excited state. If this state is highly excited, it is natural to expect a decohering evolution towards an (effective) thermal mixed state.
In order to inspect the nature of this excited state in the dual gravitational theory we consider a one-parameter family of quench profiles (3.1) with varying β ∈ (0, ∞). For definiteness, we focus on the same ω b = 2, = 0.09 case as before. Recall that our findings from the previous section established that, in the β → ∞ limit, the constructed state was the stable ground state SPS at ρ o = 0.33 given by the lowest red dot in Fig. 4. Instead, for β < ∞, the results of our simulations vindicate the existence of a critical value β c separating two radically different regimes.
• For β ≥ β c , the system always remains regular in the driving regime, and settles down to a time-periodic geometry. For finite β, in addition to the harmonic response given by e iω b t , it develops an additional periodic modulation. This can be seen in several one-point functions such as m and | O |, as illustrated in Fig. 13. This additional modulation is tantamount to the fact that our system is not in the ground state anymore, in agreement with expectations. As β is lowered two phenomena can be seen from the plots. On one side, both the amplitude and the periodicity, T , of the modulation grow. At the same time, the injected mass also increases and ends up oscillating around larger mean values. This reflects the fact that the quench has injected more energy and, consequently, the state is more excited. On the other, periodic modulations become less and less harmonic, and start developing plateaux where the system stays for progressively longer times as β → β c . This can be observed in Fig. 13 as we move to plots on the lower right part. Interestingly, the mass curve in the last plot matches very well with the corresponding one in Fig. 7. As we will show in the next section, this is more than a coincidence, and indeed, in the limit β → β c the end result of the non-quasistatic quench is the SMS to which the unstable SPS with the prescribed values of (ω b , ) will decay when perturbed with the appropriate sign of the most relevant fluctuation. To be most transparent, and referring to Fig. 4, the implied unstable SPS is the top red dot sitting at ρ o = 0.56 along the vertical dashed line. In summary, a quench to the same final boundary data ( , ω b ) yields the stable solution (lower red dot) if performed in a quasistatic way β → ∞, and the SMS associated to the unstable SPS (higher red dot) in the critical limit β → β c . However, the dashed line should not cause confusion in the later case. The process, not being quasistatic, has not proceeded through a sequence of SPS's. It cannot be drawn as a path in the solution space. In Fig. 14, we plot the period T of the periodic modulation of the final solutions as well as the average mass as a function of β. Both functions increase monotonically with decreasing β. For β → ∞, m T tends to its value on the stable SPS, while T is compatible with the period of the fundamental normal mode of this SPS.
• For β < β c , the quench results in an energy injection process strong enough to trigger gravitational collapse. After a transitory regime both m and | O | increase without bound, while the value min x [f (t, x)] drops to zero, signaling the approach to an apparent horizon. The duration of this transitory regime grows with smaller |β − β c |. After the collapse, the harmonic driving keeps injecting energy continuously into the system, which increases its mass monotonically. Representative examples of the behavior just described can be found in Fig. 15.
A pertinent question that cannot be definitely answered with our numerical methods is whether the system reaches an infinite mass state in finite or infinite field theory time. 6 A careful analysis of this post-collapse regime should connect with the findings in [6]. We have indeed pinned down different regimes for the growth of the one-point functions and refer the interested reader to Section 6.

The unstable attractor
We have clearly demonstrated that, for non-quasistatic build-up processes, the major property is the existence of the time scale β c , which separates two different phases during the driving regime. These phases are distinguished by the final fate of the solution. The new time scale seems to be related to the existence of an intermediate time attractor, and we mentioned that this attractor was nothing but the unstable SPS associated to ω b and . In the following discussion we provide evidence in favor of this statement.
For concreteness, let us consider the β = 21.222683 < β c case. In Fig. 16a, we plot ρ(t 0 , x), f (t 0 , x), as obtained from the numerical simulation at t 0 = 26.075. In Fig. 16b, these functions are compared with their values on the unstable SPS, ρ SP S (x), f SP S (x), where we plot the differences at t 0 . Clearly, these differences are small.
In order to gain further understanding, let us compute the time evolution of the norms of δρ and δf , defined by 7 We plot these quantities in Fig. 17. It is clearly seen that there is a time window in which the distance between the fields in the driving phase and the unstable SPS is negligible.
In Section 4, we employ the identification of the intermediate attractor and the unstable SPS to argue that the results presented in this section can be naturally understood as a dynamical type I phase transition. There are also apparent differences between our findings and the type II phase transition uncovered by Choptuik for a massless scalar field [30]; actually, our results are akin to the type I phase transitions studied in asymptotically flat space [17][18][19]. In the first of these works, the authors considered a SU (2) non-Abelian gauge field minimally coupled to gravity. There, for certain one-parameter families of initial data, a static soliton -the Bartnik-Mckinnon soliton-8 played the role of an intermediate time unstable attractor. The existence of this attractor entailed that, for supercritical data and as the critical surface was approached, the mass of the black hole formed did not display Choptuik scaling: there was a mass gap, set precisely by the soliton mass.
As summarized in [18], there are two model-independent assumptions that must be satisfied in order to have a first-order phase transition in spherically symmetric gravitational collapse, which we state here for completeness: • Assumption a. There exists a static, regular solution with only one unstable eigenmode. Let Ψ u (x) denote a generic field of this unstable geometry. Having only one unstable eigenmode, the evolution of any small spherically symmetric fluctuation around Ψ u (x) can be decomposed as where λ ∈ R + is the eigenvalue associated to the unstable eigenmode. Modes with n > 1 have Re λ n ≤ 0.
• Assumption b. The final fate of the perturbation δΨ u (t, x) depends solely on the sign of α. For one sign, a black hole forms; for the other, gravitational collapse does not take place. 9 As pointed out in [18], assumption a means that the stable manifold W S of the solution Ψ u has codimension one. Assumption b entails that the stable manifold is a critical surface that divides the phase space into collapsing and non-collapsing initial data.
Consider now a one-parameter family of initial data Ψ 0 (x; p) that crosses the critical surface W S at p = p * . Ψ 0 (x; p * ) flows to Ψ u (x) along W S . By continuity, initial data Ψ 0 (x; p) with sufficiently small |p − p * | have a time development Ψ(t, x; p) that remains sufficiently close to W S so as that, after some time t 1 , the linearization (4.1) holds. The precise numerical values of the coefficients α, α n are set by the initial data: α = α(p), α n = α n (p).
Given that by definition α(p * ) = 0, we must have that and, as a consequence The above linearization is valid up to some time t 2 > t 1 at which the unstable eigenmode becomes dominant and the solution is repelled away from the critical surface. Therefore, at t = t 2 we have that where we have defined ∆t = t 2 − t 1 . This last condition forces ∆t to scale as In Section 3, we have analyzed a one-parameter family of build-up processes with different time spans β, but the same final harmonic driving, characterized by ρ b = 0.09 and ω b = 2. We found out that there exists a critical β c that divides the post-quench response of the system into two different regimes, which correspond to gravitational collapse (β < β c ) or the establishment of persistent and exactly periodic oscillations in the system (β > β c ). Furthermore, we demonstrated that, as β → β c , at intermediate times, the system is attracted to a SPS with just one unstable eigenmode. This observation, together with the fact that the nonlinear evolution of a perturbed, unstable SPS falls generically into two different dynamical regimes depending solely on the sign of the perturbation (as shown in Section 2), implies that assumptions a and b apply naturally to the problem we are considering. As a consequence, we expect the scaling (4.5) to hold, where now p, p * correspond to β, β c . In particular, the lowest eigenfrequency of the unstable SPS, λ, has to be extractable from our numerical simulations.
In Fig. 18 we plot ∆t versus | log(β − β c )| for a series of quenches with different values of β > β c . The linear fit yields a value of λ = 0.542, matching the value of λ we have obtained independently by means of the numerical solution to the eigenmode equations.

Dynamical construction of a boson star
We are starting to get operational control over the phase diagram of periodic solutions. Armed with a quench protocol like (3.1) we are able to target both at SPS's as well as at SMS's by appropriately tuning the parameters ω b , and β. In the cases we have analyzed so far, the end result is some sourced periodic solution, whether modulated or not. An intriguing legitimate question is if we can engineer a protocol such that the end solution is unsourced, yet not AdS but, instead, a boson star.
Boson stars are dual to excited QFT states that, despite being unsourced, exhibit a periodic one-point function. 10 It is more than tempting to try to establish a link between this type of solutions to the so called time crystals. Strictly speaking, these last systems are postulated to be true ground states, yet to break time translation invariance spontaneously [15]. However, 10 The boson star can be envisaged as circular polarized solution for a pair of real scalar fields rather than a complex field. The one-point function of each of these real scalar fields varies harmonically in time. A natural extension of the model would include a gauge field with which one could gauge away the phase rotation by a gauge transformation linear in time. This is a different scenario which also deserves a thorough study. However, the single real scalar field will also exhibit time periodic unsourced solutions. In this case no gauge transformation is available to gauge the periodic motion away. their existence as ground states has been ruled out in [31]. 11 The no-go theorem leaves open the possibility of finding spontaneous time shift symmetry breaking in excited states (see [16] for a review with references).
In parallel to what just described, there is another possible connection to time crystals. Due to interactions, upon a small periodic driving certain systems can respond sub-harmonically; they are dubbed discrete time crystals or Floquet time crystal. Such proposal has been realized experimentally in spin systems in which the response beats at twice the driving periodicity [4,32,33].
We only take these as inspiring considerations trying to highlight similarities and differences. As a major target, we try to devise a quench protocol that connects continuously vacuum-AdS 4 with a boson star solution. The most naive guess immediately faces a seemingly unavoidable obstruction. Imagine, for concreteness, drawing a straight vertical line at ω b = 2.5 that connects AdS with the BS curve (i.e., lying within region I in Fig. 4). Physically, this corresponds to cranking up, from AdS, a periodic driving source with rising module,ρ b (t) > 0, at a constant frequency ω b . If performed quasistatically, β → ∞, the quench protocol of the previous section moves the state along the vertical segment drawn in the figure. As explained before, this can only proceed up to the point where the line exits the gray (stability) region. If once there one continues increasing the source amplitude, unavoidably will exit the surface and face collapse and thermalization. If one starts decreasing quasistatically the source amplitude, the system will proceed backward through the same states until AdS is again recovered. It is impossible to slide back along the vertical line in the white (unstable) sector down to zero ρ b . Bosons stars look unreachable from this side of the phase diagram by means of a quasistatic process. We will provide two ways to get around this caveat. Needless to say, one of the assumptions made in the previous argument will have to be relaxed. For example, if one insists in using a quasistatic quench, one may still reach the BS curve from region II without crossing any instability curve. A remarkable possibility arises by employing, instead, non-quasistatic processes. This takes us out of the surface of SPS's, but an unexpected attractor will bring us back to it.

Quasistatic method
In this section we explicitly construct a boson star, starting from AdS, via a quasistatic quench. As advertised above, the clue to this procedure comes from drawing a quasistatic path on the space of SPS's across region II, always staying within the region of linear stability. It is fortunate that boson stars lie at the boundary of such region. The advantage of this method is that we control the frequency of the final boson star. We must start by extend the quenches considered in the previous sections, and allow now for time dependence both on the modulus and phase of the driving Note that, although ω b (t) has the quench shape of Fig. 9 (interpolating between ω i and ω f ), the actual instantaneous frequency of the source is given rather by the time derivative of its phase, ω ins (t) = d/dt(ω b (t)t) = ω b (t)t + ω b (t), which has the same asymptotic values as ω b (t). This quench process starts from the global AdS 4 vacuum at t = 0; we expect it to land on the stable boson star with ω b = ω f at t = 2β when β is sufficiently large (modulo small corrections that vanish in the β → ∞ limit).
As an example, we consider (5.2), (5.3) with ω i = 3.1, ω f = 2.9, m = 0.01 and β = 1000. See Fig. 19 for a plot of (t) and ω ins (t) in this case. On the other hand, by solving the ODE system that determines the boson star solution with ω f = 2.9, we obtain m BS = 0.0837265, The agreement between them is rather good, in particular the relative error between the values obtained from the numerical simulation and the solution of the ODE system is below 1%. Another useful quantities to monitor the difference between the end product of the quench and the stable boson star are ∆ρ, ∆f -cf. equations (3.6),(3.7)-. We plot them in Fig. 21a. By looking at Fig. 21a, we clearly see that, although small, the discrepancy between the end state of the quasistatic quench process and the target boson star is not negligible. It is natural to wonder about the origin of this difference. Is the boson star oscillating at a frequency different than the target frequency? Does this difference vanish when β → ∞?
Let us address the first issue. In our example we have that, for t > 2β, ρ(t, 0) = 0.230379. The boson star associated to this number has m BS = 0.0843514 and | O | BS = 1.22635. As we see, these quantities agree with the results of the numerical simulation in five significant digits. On the other hand, the frequency of this boson star is ω BS = 2.8992 < ω f = 2.9. The simulated solution and the actual boson star the system equilibrates to are compared in Fig. 21b. In this case, after the quench both ∆ρ and ∆f are compatible with the numerical noise.
The discussion in the previous paragraph shows that the frequency of the boson star obtained after the slow quench does not need to agree with the target frequency. We expect that the difference between both frequencies vanishes in the strict β = ∞ limit, i.e., for a quasistatic quench. Equivalently, we expect that, when we employ the target boson star to perform the comparison, ∆ρ, ∆f → 0 as β → ∞.  This behavior makes plausible to assume that in the strict quasistatic limit we will recover the target boson star.

Non-quasistatic method
As advertised before, we can also reach BS's with a quench that starts from region I but, in turn, we must give up quasistaticity. The crucial observation is that, if instead of using a large value of β, we let β → β c , we reach a SMS associated to an unstable SPS. Once in the SMS, which is stable by definition, we may slowly turn off the source amplitude. If the end result is a stationary solution that corresponds neither to a BH nor to AdS, the only remaining option is a BS. We shall provide numerical evidence that this is indeed the case.
We construct the SMS by the procedure described in section 3.2. We work with the source profile (3.1) with β = β in β c , and wait until a time t out ≥ β in after build up. Once the SMS has formed, we start to turn off the source amplitude, reaching zero at t = t out + β out . The extinction quench follows the time reversed profile with, now, a larger time scale β out As an example, in Fig. 23 we have considered a build-up phase characterized by ω b = 2.5, = 0.09 and a quench-in time β in = 23.811857278, which is the closest to β c that we have managed to approach from above. As shown in Fig. 23, the system first "jumps" on top of the unstable SPS attractor, where is stays for more than 100 seconds, until it falls onto the side of the modulated SMS, as expected. After the first beating, we begin to turn off the source amplitude slowly, with a time scale β out = 1000. The final, sourceless solution has a non-vanishing mass and frequency ω consistent with some BS.
A very important difference between the setup at hand and the quasistatic method described in the previous subsection is that, in the present case, even for large β out , the frequency of the reached BS, ω BS , has no obvious relation with the driving frequency, ω b , that has been used throughout the quench.
A second very important point is that, presumably due to the violence of the build-up phase (i.e. the smallness of β in ), the system does not return to a strict BS state after the source has been turned off: the final geometry can be understood as the linear superposition of a BS and its first nontrivial eigenmode. Our numerical experiments point to the fact that this conclusion holds irrespectively of the value of β out . In fact, a stronger statement seems to be true: the final BS is actually independent of β out and there is a univocal relation between this BS and the SMS generating it; β out only controls the amplitude of a final oscillation about the final BS, which occurs with the frequency of its first nontrivial normal mode. An interesting open question is determining which property of the original SMS sets ω BS .
In order to strengthen these statements above, we have considered a set of non-quasistatic quenches leading to a BS, now with ω b = 2, β in = 22, and quench-out time scales β out = For normal mode fluctuations around this BS, the first nontrivial eigenfrequency is λ 1 = 1.975822. According to our discussion so far, we expect that, after the quench out, the time evolution of the ρ(t, x) drawn from the simulation is given by where χ 1 (x) is the spatial profile of the BS first nontrivial eigenmode, and e(t, x) an error term that takes into account both linear contributions from higher eigenmodes as well as possible nonlinear corrections. The parameters α and δ are fixed by demanding that (5.8) and its first time derivative (both without the error term) hold exactly at the origin at t = t out + β out . 12 We define After the source has been turned off, we expect that ∆(δ 1 ρ) 1 (i.e., we are close to the BS solution) and, furthermore, that ∆(δ 2 ρ) ∆(δ 1 ρ) (i.e., the leading order deviation with respect to the BS solution is controlled by its first nontrivial eigenmode). Both of these expectations hold, as the reader can convince herself by looking at Fig. 24. after the quench-out. The three magnitudes clearly display the hierarchy described in the main text. Right: ∆(δ 1 ρ) (brown) and ∆(δ 2 ρ) (orange) during the whole post-quench-out regime for β out = 5000. We clearly see that the inequality ∆(δ 2 ρ) ∆(δ 1 ρ) holds.

The post-collapse regime
Our numerical techniques allow us to follow the system after gravitational collapse has taken place, at least during some time.
In the examples we have analyzed, we have not found traces of equilibration to a stationary regime: the harmonic driving always led to a monotonic increase of the energy density. This increase is responsible for a steady growth of the apparent horizon, which gets progressively closer to the boundary. 13 The way in which the energy density increases is not universal, and depends on the specific form of the harmonic driving.
Let us discuss some examples in detail. The first one involves a perturbation of the unstable SPS with ρ b = 0.01, ω b = 2.375 by its first, unstable eigenmode. In Fig. 4 this point is to be located to the right of the upper red dot in region I, on the isocurve with ρ b = 0.01. As we have illustrated in Section 2, choosing the right sign of the perturbation leads right away to gravitational collapse. Following the post-collapse evolution of the system, we find the results depicted in Fig. 29: | O | stays approximately constant, 14 while the energy density increases in an approximately linear fashion, as prescribed by the differomorphism Ward identity given in eq. 2.4.
As explained in [6], it is useful to portrait the system trajectory in the (φ b , O ) plane in order to understand its response to the harmonic driving. We plot the real section of this complex curve in Fig. 26a. It is clearly observed that the trajectory remains bounded. In when perturbed by its fundamental, unstable eigenmode with amplitude = 0.001. We start at t = 30. Each period of the source has a different color, whose wavelength increases the later it starts. Right: Relative value of the imaginary part of nonlinear conductivity σ. A nonzero value indicates that the system response is partially in-phase with the driving.
order to understand the phase alignment of the source and the response, it is convenient to introduce the nonlinear conductivity We plot the relative contribution of σ in to the total conductivity in Fig. 26b. It is clearly observed that, while small, it typically has a non-negligible value. The boundedness of the response, together with the fact that it is not in complete phase opposition with respect to the source, leads to identify the post-collapse evolution of the system as belonging to the dynamical crossover tilted regime obtained in [6]. with = 0.09, ω b = 2 and β = 20. We start at t = β. Each period of the source has a different color, whose wavelength increases the later it starts. Right: Relative value of the imaginary part of the nonlinear conductivity σ. A nonzero value indicates that the system response is partially in-phase with the driving.
The second example of post-collapse evolution we would like to discuss involves a SPS with ρ b = 0.09, ω b = 2, perturbed by its fundamental, unstable eigenmode. In Fig. 4 this is precisely the upper red dot in sector I. Since this solution controls the late-time dynamics of the non-quasistatic quench processes we discussed in depth in Section 3, we can instead focus on those. For definiteness, we take β = 20. In Fig. 27, we plot the time evolution of m and | O |, while in Fig. 28a we depict the phase diagram (along the lines of Fig. 26a).
It is clear that this late-time regime is qualitatively different from the previous one: | O | does not remain bounded and, as a consequence, m increases faster than linearly. On the phase portrait, we observe a clear precession of the system's trajectory, which is accompanied by an increase in its amplitude. This precession is due to the fact that the response is not entirely out-of-phase with the source, as illustrated by Fig. 26b. This new late-time phase is best identified with the unbounded amplification regime found in [6], 15 as the trademark of this regime is the unbounded and partially in-phase response to the harmonic source.  We close this section by providing one last example. It corresponds again to a quench profile of the form (3.1), this time with frequency ω b = 20π, amplitude = 10 −4 and time span β = 1. As illustrated in Fig. 29, the energy density grows linearly after the quench, while the absolute value of the scalar vev remains constant. In accordance with these facts, the phase diagram of the system corresponds to a closed, non-precessing trajectory of fixed amplitude (see Fig. 30a). Consistently, the imaginary part of the nonlinear conductivity vanishes, i.e., the system response is in complete phase opposition with respect to the harmonic source (see Fig. 30b). The features discussed so far imply right away that the system finds itself in the linear response regime.

Periodically driven real scalar field
We now turn the attention to the case of a real scalar field. Now, unlike the complex case, the metric of a time-periodic real solution will not be static. This implies a dramatic change in the formalism and the methodology. Naturally, what makes the problem tractable is timeperiodicity itself which, for instance, allows us to replace the linear analysis of normal modes by the equivalent Floquet analysis. Qualitatively, the final results we obtain are surprisingly close to those of the complex case.
Another important difference with respect to the complex case is the fact that the solution, albeit being periodic, is not harmonic. For the complex scalar the time dependence can be factorized as a complex exponential -cf. eqn. (2.3)-, but for the real one the time periodicity requires an infinite spectrum of Fourier modes given by the following ansatz Only the boundary source remains harmonic, φ(t, π/2) = ρ b cos ω b t. As before, ω b = 2π/T is the driving frequency in the boundary gauge δ(π/2) = 0 (the techniques needed to obtain the fully nonlinear solutions are reviewed in appendix D). Although we cannot factorize the radial dependence into a function ρ(x), we will keep the same notation as for the complex case, i.e., ρ o ≡ φ(0, 0) and ρ b ≡ φ(0, π/2). Fig.31 shows an example of the Fourier spectrum of these solutions.
As in the complex case, the space of sourced periodic solutions -SPS-is a surface embedded in the three dimensional space spanned by coordinates (ω b , ρ b , ρ o ). In the ρ o → 0 limit the periodic solutions go to (sourced and unsourced) linear modes of AdS 4 , while the limiting, unsourced cases with ρ b = 0 are sometimes termed nonlinear oscillon solutions in the literature, and were first constructed in [34] (see also [24] for an extensive account of results and methods). Figure 32 is the counterpart of Fig. 2 in the real case. The similarity is remarkable. The quality of the plot is substantially smaller, as oscillatory solutions are much more demanding in terms of computational resources than stationary ones. This implies   Fig. 2, to which the similarity is evident. The range exhibited is considerably smaller in this plot and this is due to the higher technical difficulty in computing each solution individually.

Linear stability
The study of stability is more involved now than in the complex case. There, a set of static equations was perturbed, and the linearized spectrum of normal modes was easy to establish. In the present situation, the linearized fluctuations obey a time-dependent system of equations. The pertinent tool to employ now is provided by the Floquet analysis. In our situation we consider linearized fluctuations for the fields of the form where the subindex p stands for "periodic solution" andx denotes the perturbations. At first order in the perturbations amplitude we obtain a system of the forṁ where L(t) = L(t + T ) is a time periodic linear integro-differential operator (in x) built out of the solution x p . The Floquet theorem establishes the existence of a solution of the form The information about the stability resides in λ i , with i = 1, 2, . . . In the general case, both λ and P will be complex. From equation (7.3), we immediately see that λ * and P * are also solutions. The real solution thus involves the appropriate linear combination of both. Also, though not evident, one can show that solutions come in pairs, with eigenvalues of opposite signs, ±λ. This is ultimately a reflection of the Hamiltonian character of the dynamical system (7.3). Hence the stability analysis collapses to one of two possible cases, depending on whether Re(λ) is zero or not. The first case yields a stable (cycle) solution, and the second one an unstable (saddle point) one.
Computing λ i requires an exact integration of (7.3) over a full period T . The same techniques used to obtain the exact nonlinear SPS's can be applied here; the reader can find further details in appendix D. The numerical analysis gives the region of stability that can be observed in Fig. 33, whose similarity with Fig.4 is manifest. The linear analysis says that this similarity goes beyond the shape, also at boundary of the stability region we find the same algebraic structure, either the lowest eigenvalue λ 1 becoming purely real or two real eigenvalues fusing and developing imaginary components.

Nonlinear stability
As for the complex case, analyzing nonlinear stability requires to study the full numerical evolution of the system departing from the perturbed unstable solution. The main conclusions remain similar to the ones obtained in the case of a complex scalar. The generic end result of an initial unstable SPS is a black hole. However, in the lower part of the white wedge in sector I unstable solutions decay either to a black hole or to a sourced modulated solution, SMS, depending on the sign of the single unstable linear perturbation. As for the complex case, the pulsation is a modulation of the driven oscillation, with sudden beats separated by a substantially larger plateaux.
Like the Boson Stars, nonlinear oscillons become unstable beyond the turning point for the mass [24]. At such high values of the field at the origin, the SMS's are far more fragile than their complex counterparts and usually decay after a couple of modulating beats. This may come from the fact that the whole geometry now shakes. Or maybe from the fact that the harmonic ansatz for the source starts conflicting with the anharmonic response of the field in the interior and calls itself for a generalized periodic ansatz.
Similarly to the complex case, for sufficiently small values of the source, ρ b ∼ 10 −3 , we have also found exotic real SMS's where the bounce occurs in both directions (i.e, towards lower or higher masses) when the initial perturbation is very weak. We plot an example in the second line of Fig. 34.

Type I phase transition in the real case
In the real case, the ingredients necessary for the assumptions a and b to hold are present (see Section 4). Therefore, we expect the same type I phase transition we have uncovered in the complex case. In order to confirm this expectation, we consider a time-dependent source given by the real part of eq. (3.1) (i.e., we replace exp(iω b t) → cos(ω b t)).
In Fig. 35 we plot m and O for the one-parameter family of build-up processes of the harmonic driving ω b = 2.15, = 0.09. We clearly see that there exists a β c that separates two different late-time phases, as for the complex scalar field. Furthermore, as β → β c , the system seems to spend progressively larger time around an intermediate attractor.

Dynamical construction of a nonlinear oscillon
Nonlinear undriven real SPS's (nonlinear oscillons) are the equivalent to BS's in the real case. In this section we will show how nonlinear oscillons can be dynamically constructed, in analogy to the results we presented in section 5.1 for the BS case. Same as there, here we will describe a quasistatic quench connecting the AdS 4 vacuum with a nonlinear oscillon across the linearly stable part of region II. The protocol requires the appropriate tuning of the source parameters, frequency and amplitude, as a function of time The specific profiles for (t) and ω b (t) we shall employ are again given by (5.2)-(5.3), and plotted in Fig. 19.  region corresponds to a highly dense number of fast oscillations compared with the duration of the process. In the first plot this region is not appreciated but as the inset shows it is also present.
As a particular example, we consider a quench protocol given by ω i = 3.1, ω f = 2.9, m = 10 −3 and β = 1500. Figure 36 contains the time evolution of m and O along the process. The main difference with respect to Fig. 20 is that now we are plotting the oscillation of real quantitities. The fast oscillations pile up and, due to their high density, form the blue shaded area visible in the figure. Now, instead of (3.6)-(3.7), we define two new quantities, (∆φ,∆f ), which measure the distance to the nonlinear oscillon (φ NLO , f NLO ) built with the target data, In Fig. 37 we have plotted the time evolution of these two quantities for our particular example. On the left plot we observe that, after the quench, the geometry is very close to the oscillon with frequency ω NLO = 2.9 but, as also happened in the complex case, even more to the one with ω NLO = 2.89992 (see right plot).

Summary and outlook
In this paper we analyzed periodically driven scalar fields on global AdS 4 . This framework allows one to study different aspects of holographic Floquet dynamics, such as dynamical phase diagrams and late-time regimes alternative to thermalization.
We constructed zero-temperature solutions subjected to a constant-amplitude periodic driving of the scalar and dubbed them Sourced Periodic Solutions (SPS). They are dual to so called Floquet condensates. We have characterized the SPS solution space in detail and studied throughly both linear and nonlinear stability properties. SPS extend beyond linearity known perturbative solutions about AdS.
The unstable SPS's feature a rich phenomenology upon time evolution. In particular, there exist horizonless stable late-time solutions which evade gravitational collapse and develop instead a pulsating behavior. We named them Sourced Modulated Solutions (SMS). SMS's are themselves stable solutions where the modulation pulsation impacts in the vev profile and the total mass while, remarkably, its frequency is not imposed by that of the driving.
We addressed various types of quench processes concerning both the amplitude and the frequency of the scalar source. We focused on both slow and fast quenches. The first allow the study of quasistatic processes and showed that they can be used to prepare SPS's, dual to Floquet condensates, starting from the AdS vacuum. The study of non-quasistatic quenches uncovered a nice surprise: by suitably fine-tuning the quench rate, the system stays for an arbitrary long time on an unstable SPS and then decays either to a SMS or a black hole. These SPS's act therefore as attractors in a Type I gravitational phase transition.
Next, we examined the possibility of using such quenches to prepare boson stars starting from AdS. From the gauge/gravity duality, this is motivated by the exciting possibility of preparing experimentally sourceless Floquet condensates in strongly coupled systems. In fact, we find explicit solutions to this problem using both quasistatic and non-quasistatic quench protocols.
Finally we have studied the post-collapse evolution. We have unravelled the three basic behaviours found elsewere [6].
The analysis has been performed both for the case of complex and real scalars. We showed that the phenomenology is qualitatively similar despite the technical treatment is different (and considerably more involved in the real case). We mainly focused on massless scalars, but preliminary results for the massive situation with m 2 = −2 point to the persistence of the same picture.
The complex field setup in the boson star configuration can be regarded as a time-like version of the spatial Q-lattice model [35]. These models have been recently shown to provide a framework where to study spontaneous symmetry breaking of space translations [36][37][38]. This observation gives support to the speculation that relates boson stars to time crystals.
The statement as such could be affected by subtleties which need to be carefully studied. It is however interesting to note that the holographic setup could avoid the known no-go theorems [31,39]. Indeed, the no-go theorems rely on standard arguments about the largevolume thermodynamic limit, while holography in global AdS is related to an alternative thermodynamic limit at finite volume.
The precise origin of the frequency of the pulsating solutions is an interesting open problem to be addressed. It resembles the possible spontaneous generation of a time scale observed for similar late-time solutions in the pumping setup [40], although in the present case the pulsating frequency does depend on the amplitude of the initial perturbation of the unstable harmonic solution. There is a lower bound to this amplitude set by the numerical noise, hence we cannot establish whether in the limit of vanishing amplitude the period of the modulating pulsation diverges or not.
While entering the final stages of this work, paper [11] appeared. Albeit in a totally different context, their results bear intriguing resemblance with ours. In fact their Fig. 3 and our Fig. 2 are very similar. Our boson stars and nonlinear oscillons are, in their language, vector meson Floquet condensates. However, the physics behind looks very different. In their case the instability is related to a Schwinger pair dielectric breakdown. Also their process of building a sourceless Floquet condensate differs from ours. We intend to investigate further these similarities by going to richer models for a scalar field in AdS involving potentials which have a phase transition.
specific action (2.1) corresponds to taking d = 3 and V (|φ|) = −m 2 φ * φ. As the scalar field is complex, the action is invariant under the global U (1) transformations φ → e iα φ. The equations of motion are The isotropic ansatz for the metric of AdS d+1 can be expressed as follows with x ∈ [0, π/2). A useful field redefinition is given by Upon this redefinition, the scalar field equations of motion can be cast in the forṁ From the Einstein equations we obtain, after setting κ 2 = (d − 1)/2 as well as l = 1, In the main text we consider the time-periodic ansatz (2.3) where ω is an angular frequency and ρ(x) is a real function. We are adhering here to the boundary gauge δ(π/2) = 0 that makes t equal to the proper time of an observer at the boundary. For any given ω there are static solutions for ρ(x), f (x) and δ(x) to be obtained from the equations of motion δ + sin x cos x ρ 2 + ω 2 e 2δ f 2 ρ 2 = 0 , that result from the above-mentioned ansatz.

B The role of the pumping solution
In a previous work [40], some of the authors analyzed a massless, real scalar field in AdS 4 subjected to a linearly rising pumping source, φ(t, π/2) = α b t. One of the major findings was that there exist extremely simple pumping solutions, in which the scalar field profile is flat in the radial direction, φ(t, x) = α b t, while the metric functions f (x), δ(x) are nontrivial, but static.
Our objective in this appendix is to comment on the interplay between these pumping solutions and the four-dimensional SPS's we have found in this paper. From Fig. 2 it is apparent that the surface of SPS's extends to values ω b = 0 on a line ρ b = ρ o . Naively these would be related to the mentioned pumping solutions. The correct answer involves a careful scaling limit.

The scaling limit
Consider the SPS equations of motion (A.11) (for d = 3, m = 0) and set By taking the ω b → 0 limit at fixedρ b , these equations reduce to The equations are nothing but the equations of motion for the pumping ansatz (as reported in [40]), provided we identifyρ b = α b and fix ω b = 0. This simple observation establishes that the pumping solution controls the SPS's in the limit of infinite source amplitude and zero source frequency (i.e., on the upper left part of the diagram displayed in Fig. 4).
As demonstrated in [40], there is a critical pumping rate α * b above which no pumping solution exist. Numerically, α * b ≈ 0.785187. This fact, together with the relation between the pumping solutions and the SPS's in the scaling limit, leads to the prediction that, when ω b 1, the upper boundary of the region containing the linearly stable SPS's in the ω b − ρ b plane must approach the curve Note in particular that, if this expectation is correct, ρ b must diverge as ω b → 0. The analogous statement ω b ρ o = α * b also has to hold in this limit, since the difference between ρ b

C Study of the normal modes
This appendix provides details on the spectrum of linearized normal modes of the complex SPS's and BS's presented in the main text.
Linear stability is the first important piece of information one gets from the normal mode spectrum. A point-wise study of several individual solutions in the plane (ω b , φ o ) leads to establish the stability diagram reported in Fig. 4 whose details are highlighted below.
In Fig. 4, the solid lines departing from the φ o = 0 axis at ω b = 3 + 2n with n ∈ N + representing the BS branches that dress nonlinearly the AdS oscillons. Moving along the BS branches by raising φ 0 , one encounters BS's with increasing mass which are stable until the mass reaches a maximum. This points provides a Chandrasekar-like bound beyond which the BS's become linearly unstable. This is seen in Fig. 39 were we plot the evolution of the squared frequencies of the first two normal eigenmodes λ 2 1,2 . Precisely at the Chandrasekar mass, the second eigenmode λ 2 becomes purely imaginary. The second set of plots shows the evolution of the same modes along a line of SPS's with a small source in sector II (profiles with one node), running closely parallel to the BS curve. The situation here changes qualitatively, though continuously. The instability is now triggered by the first two modes fusing at positive values and developing imaginary parts of opposite sign. This now happens before the maximum of the mass is reached. The continuation of this point away from the line of BS deep into sector II draws the upper limit of the stability region displayed on Fig. 4 (in gray) on the segment ω b ∈ (3, 4). Figure 39 shows also that the lowest normal eigenmode in the fluctuation spectrum of BS's is a zero mode (in blue). This eigenmode corresponds to the zero-frequency deformation generating the BS branch itself. To verify this statement, one can compare the profile of the zero-frequency eigenfunction with the increment of the scalar field profile along the branch (see Fig. 40).
The stability diagram displayed in Fig. 4 presents a peculiarity at ω b = 4 and very small  φ o : the solutions are unstable, even though surrounded by stable regions. Such qualitative behavior is repeated for ω b = 2n with n ≥ 2 but not for ω b = 2. A direct look to the modes shows that the instability here is again triggered by the fusion of the first two modes developing opposite imaginary parts (see Fig. 41). It should be noted that the "encounter" of two normal modes does not lead always to a fusion and an instability. For instance, at ω b ∼ 2 for the SPS's with very small source the modes do not fuse and do not spoil linear stability.
When two modes encounter they can either cross or repel. Deciding which is the case could be however numerically demanding, especially at small values of the source φ b . Let us show this by means of an explicit example. Consider two SPS's corresponding to φ b = 0.2 and φ b = 0.1 respectively. In a region where ω b ∼ 1 the second and third modes approach and repel each other. The two modes however get closer as the source is lowered. Specifically, the minimal distance between two modes decreases more than linearly with the source, as illustrated in Fig. 42. The repulsion phenomenon departs from the perturbative expectation about the mode behavior. In fact a perturbative analysis at ω b 1 suggests that the eigenfrequencies λ n behave as λ n = λ * n ± ω b , where λ * n are the eigenfrequencies of AdS. This behavior is maintained also at higher values of ω b as the numerics shows, see Fig. 43. Nevertheless, the linear behavior of the modes has to break down as two modes approach in order to be compatible with a repulsion. The action is where we will take as before κ 2 = 8πG = 1, Λ = −6/l 2 . The isotropic ansatz for the metric of AdS d+1 can be expressed as follows We are interested in time-periodic solutions with harmonic boundary conditions such that φ(t, π/2) = ρ b cos(ω b t). Continuing with the same notation as in the complex case we will use ρ o = φ(0, 0) and ρ b = φ(0, π/2). The effective system of equations is conveniently expressed in terms of the following variables Working in the boundary gauge, we have that Π(t, π/2) = 0, while Φ(t, π/2) = 0 due to the asymptotic near-boundary expansion. With these definitions, the equations of motion becomeΦ The periodicity of the solution (7.1) instructs us to consider the most general Fourier series expansion. However, the structure of the equations of motion allows for a truncation to odd modes, which correspond to bifurcations emanating from single normal modes of AdS, Φ = Using the boundary conditions discussed previously (Φ(t, π/2) = Π(t, π/2) = 0), the values of the fields at x 0 are restricted to Φ k0 = Π k0 = 0. Hence, the unknowns to be fixed are ω b , ρ b , Φ ki and Π ki for k = 0, ..., K − 1 and i = 1, ..., N , while the equations are (D.3) and (D.4) evaluated at each of these collocation points (τ n , x i ). This gives 2KN equations for 2KN + 2 unknowns. We fix a numerical value for one of them, being it either ω b or ρ b , and add another equation of motion that sets the value of ρ o . Finally, the number of equations and variables match, 2KN + 1, and the discretized system yields an algebraic nonlinear system of equations that can be solved using a Newton-Raphson algorithm.
After obtaining time-periodic solutions we are interested in their linear stability. Consider linearized fluctuations of the form 16 φ(t, x) = φ p (t, x) +φ(t, x), Π(t, x) = Π p (t, x) +Π(t, x), δ(t, x) = δ p (t, x) +δ(t, x), where the subindex p stands for "periodic solution" and fields with a tilde are the perturbations. Inserting these ansatz into the equations of motion, and at first order in the amplitude, we obtain the equations for the perturbations where L(t) is a linear integro-differential operator (in x) constructed with the periodic fields. The periodicity of the background is inherited by the operator, L(t) = L(t + T ), T being the period of φ p . For x(t) = (φ,Π), the Floquet theorem establishes the existence of a solution of the form x(t) = e λt P(t) with P(t) = P(t + T ) . (D.15) The structure of (D.10)-(D.13) allows a truncated version of the Fourier expansion for P(t) in odd modes which rule the transitions of stability that we have found. It is convenient to express P(t) in two parts 16) where the time dependence is distributed as follows where we have definedF (i) ≡F (φ (i) ,Π (i) ) and used the time structure of the periodic fields (D.7). This system of equations has the property that given a solution λ, p (1) , p (2) another solution exists with −λ, p (1) , −p (2) , implying that an unstable mode exists if Re(λ) = 0.
On the other hand, if all λ ∈ iR, the periodic solution is linearly stable.
Solving (D. 19) through (D.22) is very similar to obtaining the time-periodic solutions φ p and Π p themselves. So the same techniques explained in the previous section are in order here. The differences lie in the fact that, in this case, the ansatz (D.16) has two more unknown functions (p (2) ) and the role of ω is played by λ. Namely, we discretize the problem as in (D.8) and set boundary conditions so as to study perturbations which don't modify the source (φ (i) (t, π/2) =Π (i) (t, π/2) = 0). Finally, by adding an additional equation to fix the amplitude ofφ (1) (0, 0) = 1, we obtain a system of 4N K + 1 nonlinear equations for our 4N K + 1 unknowns, which can be solved using a Newton-Raphson algorithm.