Holographic thermalization with a chemical potential from Born-Infeld electrodynamics

The problem of holographic thermalization in the framework of Einstein gravity coupled to Born-Infeld nonlinear electrodynamics is investigated. We use equal time two-point correlation functions and expectation values of Wilson loop operators in the boundary quantum field theory as probes of thermalization, which have dual gravity descriptions in terms of geodesic lengths and minimal area surfaces in the bulk spacetime. The full range of values of the chemical potential per temperature ratio on the boundary is explored. The numerical results show that the effect of the charge on the thermalization time is similar to the one obtained with Maxwell electrodynamics, namely the larger the charge the later thermalization occurs. The inverse Born-Infeld parameter, on the other hand, has the opposite effect: the more nonlinear the theory is, the sooner it thermalizes. We also study the thermalization velocity and how the parameters affect the phase transition point separating the thermalization process into an accelerating phase and a decelerating phase.


Introduction
Over the last years the AdS/CFT correspondence [1,2,3] has proved to be a powerful tool to describe gauge theories at strong coupling regime, where standard perturbative methods fail to work, by mapping them to a dual gravitational theory in higher dimensions whose treatment is tipically much easier. For that reason, it has found a wide range of applications in different areas of theoretical physics going from condensed matter systems to quantum chromodynamics. In particular, if the strongly coupled gauge theory is at zero temperature, it is well known that the dual gravity description involves a purely AdS space, while a finite temperature field theory requires the presence of an asymptotically AdS black hole. An interesting implication of this equivalence will be investigated in the present work, namely, the process of black hole formation in AdS space as the analog of nonequilibrium dynamics that leads a system towards thermal equilibrium after a sudden injection of energy.
One of the most interesting scenarios where these ideas have been applied is to describe properties of the quark-gluon plasma (QGP) formed in heavy ion colliders such as the RHIC and LHC. Recent results suggest that the QGP behaves as an ideal fluid with a very small shear viscosity over entropy density ratio (η/s) [4]. This implies that the QGP takes place at a strong coupling regime and therefore is amenable to a dual gravity treatment. Indeed, holographic calculations using the prototypical N = 4 SU(N ) supersymmetric Yang-Mills (SYM) theory at finite T and its string theoretical AdS gravity dual show that there seems to be a small and universal lower limit for the ratio η/s for all theories with gravity duals [5,6]. This is one of the most prominent predictions of AdS/CFT at the moment. However, while the near-equilibrium dynamics (e.g., transport coefficients such as viscosity and electrical conductivity [7] and more general aspects of dissipative hydrodynamics [8,9]) of the QGP is well known, the far from equilibrium process of formation of QGP after a heavy ion collision, often referred to as thermalization, is not well understood. The thermalization time scale observed at RHIC is considerably shorter than expected according to perturbative techniques [10], reinforcing the need of a strong coupling description of the thermalization process.
Some attempts were made in the literature to address this problem from a holographic point of view as a dual process of black hole formation via gravitational collapse in AdS space [11,12,13,14,15,16,17,18,19]. A slightly different and simpler model for holographic thermalization was introduced by Balasubramanian et al. in [20,21], which despite its simplicity captures many important features of the thermalization process. It consists in the collapse of a thin shell of matter described by an AdS-Vaidya metric that interpolates between pure AdS space at early times and Schwarzschild-AdS black hole at late times. The authors used the dynamical background of the collapsing shell to study the time evolution of nonlocal thermalization probes of the boundary conformal field theory with well known dual gravity descriptions in terms of geometric quantities. They considered equal-time two-point correlation functions of local gauge invariant operators, expectation values of Wilson loop operators, and entanglement entropy which correspond in the gravity side to minimal lengths, areas, and volumes in AdS space, respectively. They found that the thermalization is a top-down process (i.e., UV modes thermalize first while IR modes thermalize later), in contrast to the predictions of bottom-up thermalization from perturbative approaches [22]. This has a clear and intuitive interpretation from the AdS/CFT perspective: UV modes correspond to small distance scales in the boundary of AdS space and, therefore, they do not capture much of the details of the collapse process happening deep into the bulk. IR modes, on the other hand, penetrate deeper into the bulk and for that reason are naturally more sensible to details of the bulk dynamics, therefore they should thermalize later. In addition, the authors found that the thermalization time scales typically as t therm ∼ /2, where is the characteristic length of the probe.
A natural extension of this model was proposed in [23] (see also [24]) to include the effect of a nonvanishing chemical potential µ, which is usually the case in real heavy ion collision processes. The authors considered the collapse of a thin shell of charged matter in the bulk described by an AdS-Vaidya-like metric leading to a thermal equilibrium configuration given by a Reissner-Nordström-AdS black hole. They argue that varying the charge of the final state black hole from zero to the extremal value allows to explore the full range of chemical potential per temperature ratio µ/T in the dual strongly coupled QGP. Their main conclusion was that as the charge is increased, the thermalization time for renormalized geodesic lengths and minimal area surfaces becomes larger. Further investigations were made later to include Gauss-Bonnet higher curvature corrections [25,26], angular momentum [27], noncommutative [28], Lifshitz and hyperscaling violating geometries [30,29], de Sitter boundary field theories [31], as well as more elaborated dynamics for the collapsing shell [32,33,34] and a discussion on spectral functions of boundary two-point correlators [35]. Related work can also be found in [36,37,38,39,40,41,42].
In this paper we propose a further study of holographic thermalization with a chemical potential using Einstein gravity coupled to Born-Infeld (BI) nonlinear electrodynamics in the bulk. This is a natural generalization of the discussion initiated in [23], since it accomodates more elaborated dynamics for the gauge field including (all order) higher-derivatives of A µ and, therefore, it may give rise to interesting effects on the chemical potential of the dual boundary theory that are not captured by the Maxwell description. Although BI electrodynamics has its origin a long time ago [43] as an attempt to obtain a finite self-energy of point-like charged particles, currently a renewed interest has been raised due to recent developments in superstring theory. In particular, it is well known that the low energy behavior of the vector modes of open strings is governed by the BI action [44,45], while the low energy dynamics of D-branes is given by a similar non-Abelian version of the BI action [46] (see also [47,48]). Consequently, BI electrodynamics provides a promising scenario to explore deviations from Maxwell electrodynamics, specially from the point of view of AdS/CFT calculations where string theory plays a prominent role (see [49,50,51,52,53,54,55,56,57,58] for an incomplete list of previous works in this direction).
The paper is structured as follows. In section 2 we review the black hole solutions of Einstein-Born-Infeld theory in AdS space and construct their Vaidya-like extensions modelling the collapse of a thin charged shell, to be used in the sequence. Section 3 is devoted to the holographic setup for the non-local observables chosen as probes of thermalization, namely, the equal-time two-point correlators and the expectation value of Wilson loops. In Section 4 we give details of the numerical calculations and present all the results with the effects of the chemical potential and BI parameter on the thermalization curves and velocities. Finally, Section 5 contains our concluding remarks.

Vaidya AdS Black Hole solutions in Einstein-Born-Infeld theory
The starting point is the (d + 1)-dimensional Einstein gravity action with a negative cosmological constant Λ = −d(d − 1)/2l 2 (being l the AdS curvature radius) minimally coupled to Born-Infeld electrodynamics where L BI (F ) is given by The constant β is the BI parameter with dimension of mass. 1 It is defined in such a way that the limit β → ∞ corresponds to the standard Maxwell Lagrangian. We choose units in which 16πG = 1, G being the Newton's constant in (d + 1) dimensions.
The charged black hole solution to the equations of motion coming from the action (1), first obtained in [59] (see also [60]), reads where dΩ 2 d−1 denotes the metric on the unit sphere S d−1 and In the above equation 2 F 1 (a, b; c; x) is the hypergeometric function and M, Q are integration constants related to the ADM massM and chargeQ of the black hole via 2 1 In the context of string theory, β appears tipically in terms of the string parameter α via β = 1/2πα . 2 For simplicity, we will keep referring to M and Q hereinafter simply as "mass"and "charge"parameters of the black hole without any risk of confusion. ω d−1 being the volume of the S d−1 . There is also a purely electric gauge field given by where Φ is a constant corresponding to the electrostatic potential at r → ∞, which will be related to the chemical potential in the dual gauge theory according to the AdS/CFT correspondence. It is defined such that the gauge field vanishes at the horizon, i.e., The electric field associated to (5) is finite at the origin r = 0, which is a key feature of BI theories. The black hole function (4), on the other hand, is in general singular at the origin. Such a singularity is hidden behind an event horizon provided the free parameters are chosen so that the equation V (r h ) = 0 admits a real positive solution. We should also mention that taking the limit β → ∞ in (4) gives the well known Reissner-Nordström-AdS black hole studied in [61].
The solution (3) has the topology of R × S d−1 at the AdS boundary r → ∞. In the context of the AdS/CFT correspondence it is interesting to consider the limit where the AdS boundary is R 1,d−1 instead, since one is often interested in dual gauge theories living on flat space. This procedure is known in the literature as the "infinite volume limit", and it arises only due to the presence of a negative cosmological constant [61]. The idea is to introduce a dimensionless parameter λ (which will be set to ∞) and rescale all dimensionful quantities as r → λ 1/d r, t → λ −1/d t, M → λM, Q → λ (d−1)/d Q, β → β, l → l while at the same time blow up the S d−1 as i . This leaves the (t, r) block of the metric almost invariant (except for the contribution of the constant term in (4)). Finally, taking λ → ∞ yields where U (r) ≡ V (r) − 1. Notice that now the horizon, defined by U (r h ) = 0, is planar instead of spherical, so we should refer to (7) as a black brane instead of a black hole. In order to avoid the coordinate singularity at r = r h it will be interesting to express the metric in Eddington-Finkelstein coordinates by introducing a new time coordinate v defined by dv = dt + dr/U (r), and also it will be convenient to work with an inverse radial coordinate z = l 2 /r such that the AdS boundary stays at z = 0 while the singularity r = 0 sits at infinity. The resulting metric is where we have defined Notice that f (z) → 1 near the AdS boundary z = 0.
The Hawking temperature of a black hole in the context of AdS/CFT can be viewed as the equilibrium temperature of the dual field theory living on the boundary. It is obtained as usual by continuing the black hole metric to its Euclidean version via t = −it E and demanding the absence of conical singularities at the horizon. This results in a periodic Euclidean time t E whose period is identified with the inverse Hawking temperature. For the AdS Einstein-Born-Infeld black brane (7) this calculation gives This expression reduces to the Hawking temperature of the Reissner-Nordström-AdS black hole in the Maxwell limit β → ∞ [23]. When T = 0, the black brane is called extremal. If we think of all the parameters but the charge as fixed, then we can characterize the extremal black brane solution by a maximal value of charge Q ext given by According to the AdS/CFT dictionary, the asymptotic value of the time component A t (r) of the gauge field at the AdS boundary r → ∞ (namely, the constant Φ in equation (5)) corresponds to the chemical potential µ in the dual quantum field theory, µ ∼ lim r→∞ A t (r). Actually, the precise relation should include some scale ξ with length units since the chemical potential must have energy units (or [length] −1 ) while A µ as defined by the action (1) is dimensionless. Hence, the chemical potential per temperature ratio of the boundary field theory is given by with Φ and T given by expressions (6) and (10), respectively. A remarkable feature is that if the horizon radius r h and the BI parameter β are kept fixed, then by varying the charge from Q = 0 (vanishing Φ) to Q = Q ext (vanishing T ) it is possible to explore the whole range of values of the ratio µ/T in the dual field theory, i.e., from µ/T = 0 to ∞.
A Vaidya-like extension of the BI AdS black brane metric (7) can be constructed by promoting the mass and charge to arbitrary functions M (v) and Q(v) of the advanced time v. The resulting dynamical metric has the same form as in (8) but now with f (v, z) instead of f (z) due to the time dependence introduced on the mass and charge. The same also holds for the gauge field (5), which now becomes A µ (v, z). Such a spacetime describes the collapse of a thin-shell of charged dust from the boundary of the AdS space towards the bulk interior.
Of course such a metric is not a solution of the action (1) anymore: there must be some external matter action S m sourcing the time variation of M (v) and Q(v). If we take this external contribution into account, the Einstein-BI equations of motion become (we restore the factors of G for a moment): The Vaidya-BI-AdS metric above-mentioned is a solution to these equations provided the external sources satisfy where the dot denotes ∂ v . We notice that there is no β dependence on J ν (m) above, and indeed this is exactly the same current found in [23] in the Vaidya-Reissner-Nordström-AdS case. T (m) µν , on the other hand, differs from the corresponding one in the Reissner-Nordström case due to the hypergeometric term (but naturally reduces to it in the β → ∞ limit, since 2 F 1 [a, b; c; 0] = 1).

Holographic thermalization
In this section we are going to study the thermalization process of a strongly coupled quantum field theory whose bulk gravity dual corresponds to the Einstein-Born-Infeld system presented in the previous section. According to the AdS/CFT correspondence, a zero temperature state on the d dimensional boundary theory is dual to pure AdS d+1 in the bulk, while a thermal state corresponds to the BI-AdS black brane (8). Therefore, any dynamical bulk spacetime which interpolates between these two situations would be a natural candidate to holographically model the nonequilibrium process leading to thermalization of the boundary theory after a rapid injection of energy.
Following [20], we choose our dynamical bulk metric to be the Vaidya-BI-AdS metric discussed in Section 2, namely (we set the AdS radius l = 1 hereafter) with the mass and charge functions given by (see [62] for an interesting discussion on the corresponding bulk null energy condition) Clearly, for v → −∞ we have pure AdS (M (v) = 0 = Q(v)) and for v → ∞ we have M (v) = M and Q(v) = Q which is the BI-AdS black brane (8). Indeed, equations (18)- (19) are just smooth versions (convenient for the numerical analysis) of the step functions which represent a shock wave (a zero thickness shell of charged matter suddenly forming at v = 0). The constant v 0 represents a finite shell thickness and for v 0 → 0 we go back to the step function.
The next step is to choose a set of observables to use as probes of thermalization. Since local observables in the boundary such as expectation values of the energy-momentum tensor are not sensitive to the thermalization process, one needs to consider extended non-local observables. 3 In this work we shall focus on equal time two-point correlation functions and expectation values of rectangular Wilson loops, which have well known holographic descriptions in the bulk in terms of renormalized geodesic lengths and minimal area surfaces, respectively. A third observable that could be used is the entanglement entropy of boundary regions, which have a very similar description in terms of minimal volumes of codimension-two surfaces in the bulk. However, as the results of entanglement entropy lead essentially to the same conclusions, we will not show them here in order to avoid unnecessary repetitions.

Renormalized geodesic lengths and two-point functions
We start with the equal-time two-point correlation functions of local gauge invariant operators O(t, x) of conformal dimension ∆. The AdS/CFT correspondence provides a simple geometrical way to compute it in the bulk gravity dual when the operator O is "heavy", i.e., when ∆ 1 [63]. Namely, where L stands for the length of the bulk geodesic between the points (t, x) and (t, x ) located on the AdS boundary. 4 Actually, one should be careful when doing such an approximation because the geodesic length above is divergent due to the contribution of the AdS boundary (because the AdS metric itself diverges at the boundary z = 0). In order to extract a meaningful quantity we need to introduce a cutoff z 0 . It turns out that the divergent part of L is universal and equals to 2 ln(2/z 0 ) [20]. Then, we define a (finite) renormalized geodesic length L ren = L − 2 ln(2/z 0 ), which will be related to the renormalized two-point function O(t, x)O(t, x ) ren just as in equation (20).
We choose the coordinate axes for the boundary directions (t, x) such that the spatial separation = |x−x | between the points lies entirely over the x 1 direction, i.e., we consider space-like geodesics between the boundary points (t, x 1 = − /2, . . .) and (t, x 1 = + /2, . . .), where the ellipsis denotes the remaining coordinates (x 2 , . . . , x d−1 ) which are the same for both points. Then, by symmetry the geodesics cannot depend on coordinates other than x 1 , and we can use x 1 as the geodesic 3 Holography provides a geometric intuition for why local operators are insensitive to details of the progress towards thermalization: being local, they are only sensitive to phenomena happening in their vicinity near the AdS boundary. Thus they are not aware of the details of phenomena occurring near the thermal scale. We need observables dual to AdS quantities that probe deeper into the bulk in order to see signals of thermalization. 4 If there is more than one geodesic we should sum over them on the right-hand side.
parameter (we call it simply x hereinafter). The solutions of the geodesic equations are then given by a pair of functions v(x) and z(x). The boundary conditions at the AdS boundary z = z 0 are The length functional between the referred points follows immediately from the line element (17) as being It clearly depends on the path taken from one point to the other. The geodesic corresponds to the functions v(x) and z(x) that minimize the length and can be found by standard methods. The geodesic length (which we will call L) is just the value of the length functional L evaluated at the geodesic solution.
The variational problem simplifies by noticing that the integrand in (22) does not depend explicitly on x and, therefore, there is a conserved Hamiltonian H given by Using the conditions on the turning point of the geodesic, arising from the fact that the geodesic must be symmetric with respect to x = 0, the conservation equation simplifies to The advantage of working with the boundary conditions (24) instead of the original ones (21) is that we can use the conservation equation above to write the geodesic length (22) as Hence, after solving the equations of motion and finding the geodesic functions v(x) and z(x) for a given pair of initial conditions (v * , z * ) it is a trivial task to use the relations (21) to read the corresponding values of boundary separation and time t, as well as obtaining the corresponding renormalized geodesic length L ren = L − 2 ln(2/z 0 ) by means of expression (26).
The hard part is to find the geodesic, which means to solve the Euler-Lagrange equations for z(x) and v(x). With the help of the conservation equation (25) they can be written, respectively, as This is a set of coupled, highly nonlinear differential equations and for that reason it is quite hard to handle with analitical methods. However, for a given pair (v * , z * ) it is possible to find a numerical solution subject to the boundary conditions (24). Indeed, by solving for sufficiently many pairs of initial conditions (v * , z * ) (carefully chosen in order to give the same boundary separation and different times), we can track time after time the whole evolution of the geodesics in the Vaidya-BI-AdS spacetime. In particular, if we calculate the renormalized geodesic length of each of these solutions using (26) we will be able to see the full time evolution of L ren towards thermalization. This will be done in Section 4, where we provide a detailed explanation of the numerical procedure as well as the choice of parameters and show our results.

Minimal area surfaces and Wilson loops
We now study a second class of thermalization probes, namely the expectation values of Wilson loop operators in the boundary field theory. The Wilson loop is a non-local gauge-invariant observable defined as the path-ordered integral of the gauge field over a closed path C: where N is the number of colors and A µ is the non-abelian gauge field. Wilson loops contain useful information about the non-perturbative behavior of non-abelian gauge theories, such as whether they exhibit confinement or not. It is possible, in principle, to express all gauge-invariant functions of A µ in terms of Wilson loops by appropriate choices of the path C, but unfortunately they are in general hard to compute.
The AdS/CFT correspondence again provides an elegant way to compute the expectation value of Wilson loops of a strongly coupled gauge theory with a gravitational dual in terms of a geometrical quantity in the bulk [64]: where α is the inverse string tension, Σ 0 denotes the minimal area bulk surface whose boundary is the original contour C, and A(Σ 0 ) is the area of that surface. Σ 0 will be a solution of the bosonic part of the string action (the Nambu-Goto action), which is nothing but the area of the classical world-sheet with C as its boundary. Indeed, equation (29) has its origin in a saddle-point approximation of the string theory partition function around the classical solution, which is only valid when the dual gauge theory has strong coupling.
Now we proceed to compute the minimal area surfaces in the Vaydia-BI-AdS spacetime as described before just as we did for the geodesic lengths. We will focus on a spacelike rectangular Wilson loop on the boundary. The rectangle C can always be chosen to lie on the x 1 -x 2 plane, centered at the origin, with sides on the x 1 direction and R on the x 2 direction. One also assumes translational invariance along x 2 , such that the shape of the bulk surface depends only on x 1 and again we can use x 1 ≡ x to parametrize the functions v(x) and z(x) that characterize the surface. The boundary conditions at the AdS boundary z = z 0 are again given by equations (21).
Using the Vaidya-BI-AdS metric (17), the Nambu-Goto action (or area functional divided by 2π) becomes Notice that an obvious consequence of our assumption of translational invariance is that the length R factorizes. Since we are interested just in the dependence, we can study A N G /R instead of A N G itself and forget about R in what follows. As for the geodesics, the pair of functions (v(x), z(x)) that minimizes the Nambu-Goto action will be the minimal surface Σ 0 . The on-shell value of the Nambu-Goto action (i.e., A N G [Σ 0 ]), which we call A, will be our object of interest.
The subsequent calculation is closely analogous to the geodesics case. There is again a conserved Hamiltonian associated to (30) and we can introduce the alternative boundary conditions on the turning point of the minimal surface, which are the same as equations (24), to find an expression similar to (25) for the conservation equation. Replacing this back into (30), the on-shell Nambu-Goto action becomes which can be used to easily obtain the minimal area surface once we have solved the equations of motion and found the functions v(x) and z(x) for a given pair of initial conditions (v * , z * ). The corresponding values of boundary separation and time t can be read from the original conditions (21) as well. Here again we have to face the problem that the area A diverges due to the contribution near the AdS boundary, but we can regularize the divergent part using again a cutoff z 0 and subtract it from A to define the renormalized minimal area A ren = A − R/πz 0 [20].
The functions z(x) and v(x) are found by solving the Euler-Lagrange equations coming from the action (30). The simplified equations of motion (after using the conservation equation) are, respectively As before, for a given pair (v * , z * ) we can solve this set of equations numerically subject to the boundary conditions (21). Thus, for some chosen length of the Wilson loop, by iterating for enough pairs (v * , z * ) we can track the whole time evolution of the minimal surfaces and, in particular, of their renormalized areas A ren towards thermalization. This will be the aim of Section 4.

Renormalized geodesic lengths
In this section, we numerically solve the geodesic equations of motion (27) in order to find how the geodesic length evolves with time. Afterwards, we explore how the charge of the black hole and BI parameter affect the thermalization time. For the latter it will be convenient in the numerical calculations to use an inverse BI parameter b = 1/β instead of the original β, such that the Maxwell limit is b → 0 and increasing b accounts for increasingly nonlinear electrodynamics.
First of all, we fix the free parameters. We will take the shell thickness and AdS space UV cut-off to be v 0 = 0.01 and z 0 = 0.01, respectively. Since the effect of the number of spacetime dimensions and boundary separation on the thermalization probes has already been analyzed in previous works [20,21,23], we focus here in the case d = 4 (namely, AdS 5 space, which is dual to a 4−dimensional gauge theory) and a fixed boundary separation = 4. The mass M of the final state black brane can be expressed in terms of the radius of its event horizon using the definition of r h , i.e., the largest solution of U (r h ) = 0. Then, if we choose to fix the horizon at r h = 1, 5 the mass is given by which of course only holds provided that b and Q take values consistent with the existence of an event horizon. For a given b this means that Q is allowed to take values from Q = 0 to the extremal value (11), which now reads Q ext (b) = 2 + 3b 2 .
As we pointed out in Section 2, for each b, considering values of charge 0 ≤ Q ≤ Q ext (b) we can study all the range 0 ≤ µ/T ≤ ∞ in the dual gauge theory. It is instructive to illustrate this fact here by looking at the form of expression (12) in the present case (we choose the scale ξ = 1 for simplicity), namely, A plot of this as a function of the charge for distinct values of b is shown in Figure 1, showing that indeed all the range of µ/T is covered. We also notice that for small values of charge (up to ∼ 0.5) the BI parameter has no effect on the chemical potential since all the curves agree, so we will concern only about charges above this value in what follows.
Now we proceed to investigate the effects of Q and b on the thermalization process. We choose the test values b = 0, 0.4, 1 and Q = 0.5, . . . , Q ext (b) for the numerical analysis. The procedure is the following: for a given (b, Q) pair we solve the geodesic equations (27) subject to the boundary conditions (24) characterized by a pair of values (v * , z * ). We do this recursively for various pairs of initial conditions and collect from these just those yielding a boundary separation = 4. 6 Each of these collected solutions correspond to a different stage of the motion of the geodesics, as we may check by computing the corresponding boundary time t via equation (21). We then calculate the renormalized geodesic length L ren of each of these collected solutions and construct a list of points (t, L ren (t)) which contains all the information about the time evolution of the renormalized geodesic length. Actually, it will be convenient to divide all the lengths by in order to obtain a dimensionless, −independent quantityL = Lren . In addition, we subtract from this the final (thermal) value just to force all thermalization curves to end at zero, so that the quantity to be plotted isL −L thermal versus t.
Before showing the thermalization curves, in Figure 2 we present an intuitive simple view of the effect of the BI parameter on the thermalization. It consists of a sequence of snapshots of the time evolution of geodesic profiles as well as the shell of charged dust described by the Vaidya-BI-AdS metric to form a black brane at z h = 1 at late times, for different values of the inverse BI parameter b and fixed charge Q = 1. Each column, top to bottom, corresponds to the time evolution for a given value of b. It is found that at the early stages of the evolution, up to t ∼ 1, the value of b has little effect on the dynamics. After that, b plays a crucial role in the evolution. We see that the bigger b is (i.e., columns to the right), the sooner the black hole is formed. This is clear from the bottom line of the picture, corresponding to boundary time t = 1.82, where the b = 1 black brane has just formed while the b = 0.4 one is about to form and the b = 0 one still needs some time to do so. This is a hint that the thermalization of the dual boundary field theory occurs sooner as b increases, what indeed will be confirmed in Figure 4.
In Figures 3 and 4 we show the thermalization curves for the renormalized geodesic lengths with varying charge at fixed b and varying b at fixed charges, respectively. Instead of plotting point by point all the results obtained as described above, we find more instructive to fit those points using some polynomial function and plot the resulting curve. Details about the fits will be given below. We use dashed curves in all the plots hereinafter to highlight the Maxwell limiting case studied in [23]. The thermalized (final) state corresponding to the completely formed black brane is reached in each case when the curve touches the zero point of the vertical axis. The effect of the charge Q on the thermalization is clear from Figure 3. As Q grows, the thermalization time increases, meaning that the dual field theory thermalizes later. This had already been pointed out in [23] for the case of Maxwell electrodynamics (b = 0), and now we show that the same holds for BI nonlinear electrodynamics. Since the charge corresponds to the chemical potential, this means that the smaller the chemical potential is, the faster is the pair production and the screening effect takes over in an easier way. This is compatible with lower dimensional models, where screening effects are known to prevail over confinement [65,66]. The second, more interesting result, is the effect of the inverse BI parameter b shown in Figure 4. As one can see, increasing b decreases the thermalization time, which means that the more nonlinear the bulk theory is, the sooner its dual field theory thermalizes. This confirms our intuition coming from the analysis of the geodesic profiles and shell motion in Figure 2. Such a behavior is similar to the effect of the Gauss-Bonnet parameter on the thermalization reported in [25]. As we discuss in the conclusions, this seems to be a general feature of introducing extra derivatives in the bulk theory. The numerical values obtained for the thermalization times are summarized in Table 1.  Having smooth fit functions for all the sets of numerical data we can use them to study the thermalization velocities d dt (L −L thermal ) aiming for more details of the nonequilibrium process. These are plotted in Figure 5 (only for the cases Q = 1 fixed and b = 0.4 fixed, respectively, to avoid unnecessary repetitions). We notice from the velocity curves the existence of a phase transition point at the middle stage of the thermalization, which divides the process into an accelerating and a decelerating phase. Furthermore, we see that the phase transition point is shifted depending on the values of b and Q. Figure (4.1) shows that increasing the value of b causes a delay in the phase transition point, meaning that the accelerated phase lasts longer for the b = 1 theory. This is to be contrasted with the fact that the b = 1 theory is the first to thermalize, indicating that the dynamical process in this case consists of a slowly accelerating phase followed by a quick deceleration towards the equilibrium state. On the other hand, Figure (4.1) shows that the charge has the opposite effect, i.e., as Q increases the phase transition point arrives earlier. In other words, for large values of Q (or µ/T in the boundary field theory) the thermalization process consists of a quick accelerating phase followed by a slowly decelerating phase to the final state.  It should be stressed that in this work, in contrast to the authors of [26], we find no evidence for a negative thermalization velocity at initial times. They argue that the velocity should be negative in the very beginning of the evolution corresponding to a"quantum"stage of the nonequilibrium process, which soon becomes"classical"once the velocity becomes positive. In our case ( Figure  5) all the thermalization velocities start from zero and increase monotonically until the phase transition point, indicating that nothing particularly odd seems to happen at the initial stages of the thermalization process. This is also clear from the comparison of the numerical results with the fitting functions presented in Figure 6. The zoomed region shows that the numerical points sit all over a horizontal line for initial times (up to ∼ 0.2) and therefore there is no reason for a nonvanishing slope at such stage. For that reason, we use for our fit functions degree 9 polynomials f (t) = 9 n=0 α n t n with the first powers of t (α 1 , α 2 ) set to zero in order to ensure the strictly constant behavior f (t) ∼ α 0 up to t ∼ 0.2. This allows us to make an accurate fit of the whole set of numerical points, which after all are the ones carrying the physical information.

Renormalized minimal area surfaces
In this section, we numerically solve the equations of motion (32) for the minimal area surfaces in order to track their time evolution. The strategy will follow closely that of the renormalized geodesic lengths done in the previous subsection, so we will not repeat all the details on the numerical procedure as well as the fixing of free parameters since they are essentially identical.
Again we choose the test values b = 0, 0.4, 1 and Q = 0.5, . . . , Q ext (b) for the numerical analysis. The procedure is the same as before, i.e., for a given (b, Q) pair we solve the geodesic equations (32) subject to the boundary conditions (24) characterized by a pair of values (v * , z * ). We repeat this for various pairs of (v * , z * ) and collect just those yielding a Wilson loop with side = 2 (keep in mind that the other side R does not influence in our analysis). 7 Each of the collected solutions will correspond to a different stage of the time evolution determined by calculating the boundary time t via equation (21). Then we integrate each of the collected solutions using equation (31), subtract the universal divergent part to obtain A ren , and construct a list of points (t,Ã −Ã thermal ) to be plotted against t. HereÃ ≡ A ren /(R π −1 ) is a dimensionless quantity independent of the dimensions of the boundary Wilson loop andÃ thermal is the corresponding thermal value.
A sequence of snapshots of the time evolution of the minimal area surfaces as well as the shell of charged dust described by the Vaidya-BI-AdS metric is shown in Figure 7 for different values of the inverse BI parameter b and a fixed charge Q = √ 2. Each column, from top to bottom, follows the time evolution for a given value of b. Again we can see that up to t ∼ 1.0 the value of b has little effect on the dynamics, while at the final stages of the evolution b plays a decisive role. This is clear from the bottom row at t = 1.54, where we see that the b = 1 black brane has already formed while the other two are about to form. That illustration suggests that increasing the value of b makes the black hole form earlier, which indeed will be confirmed below.
The thermalization curves for the renormalized minimal area surfaces are shown in Figures 8  and 9 for varying Q at fixed values of b and varying b at fixed charges, respectively. All the curves are polynomial fits of the numerical points (see below) and the zero point of the vertical axis corresponds to the final state of the process, i.e., the static Einstein-BI black brane fully formed. We immediately notice that all the effects are less evident than those displayed before for the geodesic lengths (the reason why we show the insets in Figure 9) due to our choice of = 2 as the characteristic scale in the boundary, in contrast to the = 4 used for the geodesics. This just illustrates the argument made in the beginning that the holographic thermalization is a top-down process. Figure 8 shows that for a given b the effect of the charge Q is to delay the thermalization process as it is increased, that is to say, as the chemical potential in the dual field theory grows, the time needed to reach the thermal state also raises. This reinforces the conclusions drawn from the analysis of the renormalized geodesic lengths in the previous subsection. The effect of the inverse BI parameter b for fixed charges can be inferred from Figure 9, namely, the larger b is, the shorter the thermalization time is. This implies that the boundary field theory is easier to thermalize in the nonlinear case, which is in perfect agreement with our results from the previous subsection.
The numerical values obtained for the thermalization times are summarized in Table 2.  Just as we did before, we also use the fit functions to study the thermalization velocities for the minimal area surfaces, d dt (Ã −Ã thermal ), which are plotted in Figure 10 for the cases Q =  and b = 1 fixed, respectively. The plots confirm the results obtained from the geodesics in what concerns the different stages of the dynamical process. Namely, there is always a phase transition instant at the middle stage of the evolution where the process changes from an accelerating phase to a decelerating phase. Such a phase transition point is reached later as we increase b for a fixed charge (see Figure 4.2), or sooner as the charge Q is increased for a given b (see Figure 4.2). In other words, the thermalization of expectation values of Wilson loops in the dual field theory consists in a slow (quick) accelerating phase followed by a quick (slow) decelerating phase towards the final state as the nonlinearity parameter b (the charge Q, or chemical potential µ/T ) is increased.  Once more we point out that our velocity plots do not provide any evidence of a negative thermalization velocity at initial times, in contrast to [26]. This can also be seen from Figure 11, where we show a comparison of the numerical results with the fitting functions. The zoomed region shows that at the initial times, up to t ∼ 0.1, all the numerical points lie over a horizontal line and thus the curve must have vanishing slope in that region. We should mention that the fitting functions used here were degree 7 polynomials f (t) = 7 n=0 α n t n with the linear term (α 1 ) set to zero to ensure the strictly constant behavior f (t) ∼ α 0 up to t ∼ 0.1 as required by the numerical data.

Conclusions
The effect of the chemical potential and the inverse BI parameter on the thermalization time of the dual boundary field theory is studied using the Vaidya-like toy model of a collapsing shell of charged dust. The bulk spacetime is dynamical, constructed to interpolate between a pure AdS space at initial times and a charged Einstein-Born-Infeld AdS black brane at late times. We use as thermalization probes the equal-time two-point functions and expectation values of Wilson loops, which have well defined dual gravity descriptions in terms of renormalized geodesic lengths and minimal area surfaces in the bulk. Another class of nonlocal observables, the entanglement entropy of boundary regions, which also have well known holographic descriptions in terms of minimal volumes of codimension-2 surfaces in the bulk, can also be used as a third thermalization probe. However, as the results are similar, we focus only on the first two since they are enough to capture all the relevant effects.
We conclude that as the charge (or, equivalently, the chemical potential) grows, the thermalization time is also increased. The inverse BI parameter, on the other hand, has the opposite effect, i.e., the larger is the value of b, the shorter the thermalization time is. At the initial stages of the thermalization we find that Q and b have little effect on the time evolution, becoming important only from the middle stage on. We arrive at the same results independently for both the renormalized geodesic lengths and minimal area surfaces. In each case they can be seen just by looking to a sequence of snapshots of the motion profiles of the geodesics and minimal surfaces as the time goes by or, alternatively, from the thermalization curves obtained from the full numerical analysis. The effect of the charge happens to be the same found in [23] using Einstein-Maxwell theory and [26], where Gauss-Bonnet curvature corrections were included. The outcome of introducing a nonvanishing b, on the other hand, is a new result. In particular, the effect is the same as that of the Gauss-Bonnet parameter reported in [26]. Since BI electrodynamics consists essentially of higher order derivatives of the gauge field, we suggest that this may be a general feature of introducing extra derivatives in the bulk, i.e., that boundary gauge theories whose gravity dual carry more than two derivatives tend to thermalize sooner than two-derivative theories after a sudden injection of energy. It would be interesting to explore this idea with other theories possessing this characteristic in order to have definite answers.
Moreover, by fitting the numerical data with smooth functions we were also able to study the thermalization velocities associated to each curve. Although this is not rigorous, it reveals some interesting features of the dynamical process of thermalization and how they are affected by Q and b. We notice the existence of a phase transition point at the middle stage of the thermalization, which divides the process into an accelerating and a decelerating phase. The phase transition point is shifted depending on the values of b and Q. Namely, as Q increases, the phase transition point arrives earlier. This indicates that for large values of chemical potential the thermalization process consists of a quick accelerating phase followed by a slowly decelerating phase to the final state. Increasing the value of b, oppositely, causes a delay in the phase transition point. This is to be contrasted with the fact that nonlinear theories thermalize first, indicating that the dynamical process for non-vanishing b consists of a slowly accelerating phase followed by a quick deceleration towards the equilibrium state. We also show from the velocity plots that the thermalization process is monotonic and the velocity is always positive at the initial stages, in contrast to the claim by the authors of [26] that there should be a negative thermalization velocity at the very beginning of the evolution.