{'fraction_non_alphanumeric': 0.03897307898021038, 'fraction_numerical': 0.016188268853628097, 'mean_word_length': 4.25599700149925, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Iterative Fluiddynamics

An iterative scheme is presented to solve analytically the relativistic fluid dynamics equations. The scheme is applied to longitudinal expansion, transversal symmetric and transversal asymmetric (triaxial) expansion as well. Within this scheme it is possible to describe the dynamics of a strongly coupled (i.e. conformal) medium for parameters referring to heavy-ion collisions at LHC.

initial conditions. We test an iterative procedure which allows for solving fully analytically the fluid equations of motion up to a given order in the expansion of powers of the time coordinate. Such iterative procedures have been established in the context of general relativity (cf. [27]). In employing the AdS/CFT correspondence, the iterative procedure is exploited, e.g. in [28], to solve the five-dimensional vacuum Einstein equations with negative cosmological constant and to derive useful insights in the behavior of the four-dimensional finite-temperature, strongly coupled medium in the conformal limit. Here, we outline such a solver for the four-dimensional fluid dynamics in section 2 and make it explicit for ideal fluid dynamics in section 3. Various applications up to the triaxial expansion dynamics including eccentricity as basis of the elliptic flow for conditions referring to heavy-ion collisions at LHC are presented in section 4. We discuss the limitations of the method in section 5 and draw conclusions in section 6.

Description of the iterative scheme
The basic idea of the iterative scheme to be described in the following paragraphs is to transform the differential equations of interest (in our case: local energy-momentum conservation) into an infinite set of algebraic relations between the Taylor coefficients of the unknown functions, which can be solved iteratively. For achieving this, the unknown functions need to be formally Taylor expanded and afterwards the series is plugged into the differential equation. Comparing coefficients of the terms with the same power in the time-coordinate leads to the desired set of algebraic relations which is solved iteratively up to a truncation order. For this purpose, let us assume that the only fluid dynamical equations governing the dynamics of the sQGP are T µν ;ν ≡ ∂ ν T µν + Γ µ αν T αν + Γ ν βν T µβ = 0, where the Christoffel symbols Γ µ αβ take into account non-Cartesian coordinates mapping out the four dimensional Minkowski space-time. In the absence of conserved charges the symmetric energy-momentum tensor T µν is assumed to be related by a constitutive equation to energy density e, pressure p, fourvelocity u µ and additionally to gradients of these quantities in the form with ∆ µν = u µ u ν −g µν 1 being the projector orthogonal to u µ . The quantity Π µν uncovers the dissipative effects. An equation of state e = e(p) and the normalization u µ u µ = 1 closes (1) to a system of four evolution equations for four independent quantities, e.g. e and the spatial components u a of the fourvelocity. The Landau Lifschitz condition u µ Π µν = 0 provides a unique definition of the velocity (cf. [19,29] for that issue). Cauchy initial conditions are specified by giving e(θ 0 , ξ) and u a (θ 0 , ξ) on the hypersurface θ = θ 0 for a suitable time coordinate θ. For the subsequent calculations we use the coordinates ξ µ ∼ (Θ, η, x ⊥ ) being related to Cartesian coordinates (t, x) via Θ = ln 1 τ0 t 2 − (x 1 ) 2 , η = 1 2 ln t+x 1 t−x 1 with τ 0 being a reference proper-time; we choose it as the instant where the hydrodynamical evolution starts. The transverse coordinate vector x ⊥ = (x 2 , x 3 ) is equal to the Cartesian one. In these coordinates, the non-zero Christoffel symbols are and the metric tensor g µν is diag(τ 2 0 e 2Θ , −τ 2 0 e 2Θ , −1, −1). The constitutive equations can be used together with the equation of state and the normalization condition of the four-velocity to obtain the components of T µν as functions of four independent variables. We choose these variables to be T α0 , i.e.
In principle, four arbitrary functions can be used to parametrize the components of the energy-momentum tensor. Since we want to apply the iterative scheme to relativistic hydrodynamics, one could choose the energy density and the spatial components of the four-velocity as well. There are two reasons why we not do so. First, setting up the iterative scheme is most straightforward for this choice and second, we want to keep the discussion in this section general, so the iterative scheme is not only applicable for ideal hydrodynamics, but also for other systems obeying local energy-momentum conservation. However, there are some subtleties concerning energy-momentum tensors which contain derivatives as in the case of viscous hydrodynamics. For such a situation (4) needs to be modified (see appendix A). One important remark is that, when the setup is specified (e.g. by choosing to investigate ideal hydrodynamics with a certain equation of state), the functional form of f µν is fixed and all its derivatives with respect to T α0 are given. Equation (1) can be rewritten as Latin letters for indices refer to spatial coordinates ranging from 1 to 3. Since f µν depends solely via its arguments T α0 on the space-time coordinates, the spatial derivative of f µν must be transformed into spatial derivatives of T α0 . The iterative scheme we want to describe is based on the Taylor expansion where we have chosen without loss of generality coordinates for which the initial conditions are given at the manifold θ = 0. Analogue expansions hold for the Christoffel symbols and for f µν : Here f µν The time derivatives of f µν are then transformed according to the chain rule into time derivatives of T α0 (see appendix B). Inserting these expansions into (5) yields It is important to note that f µν (m) contains only time derivatives of T α0 up to order m. By comparing the coefficients, (8) can be decomposed into a series of algebraic recurrence formulas relating the Taylor coefficients: Accordingly, the (k + 1)th Taylor coefficient T µ0 (k+1) is algebraically to be calculated by lower-order coefficients. This makes an iterative solution possible. The initial conditions determine T α0 (0) . Using the equation with k = 0 one can calculate T α0 (1) from that. Applying the equation with k = 1 one can calculate T α0 (2) etc. In such a way, arbitrarily high Taylor coefficients are calculated iteratively. T α0 (0) = T α0 (θ = 0, ξ) is supposed to be a smooth (differentiable) function. The determination of highorder time derivatives is therefore via (9) reduced to calculating high-order spatial derivatives of T α0 (θ = 0, ξ). For the (Θ, η, x ⊥ ) coordinate system, the Christoffel symbols (3) are constants. Therefore only the terms with k = l contribute, and (9) reduces to 3 Ideal fluid dynamics The scheme described in the preceding section is, in principle, applicable to any physical system obeying local energy-momentum conservation. Therefore, it can be used to calculate solutions of the relativistic hydrodynamical equations. For ideal fluid dynamics, where Π µν in (2) is neglected, the relation (4) becomes explicitly with f = g 00 5A + 2g 00 T 00 2 − T 00 √ g 00 −5A + 4g 00 T 00 2 3g 00 T 00 + √ g 00 −5A + 4g 00 T 00 2 A , p = g 00 T 00 2 + A g 00 By definition f µ0 = f 0µ = T µ0 . The four-velocity components u µ are given by with h = 5A + 2g 00 T 00 2 − T 00 √ g 00 −5A + 4g 00 T 00 2 The above equations hold for the conformal equation of state e = 3p, but can be generalized to others in a straightforward manner. They are obtained by considering the (µ 0) components of (2). These four equations are solved for the pressure p and the spatial components of the four velocity u a applying the equation of state and the normalization of the four velocity. Therefore, one can express p and u a as functions of the components T α0 of the energy-momentum tensor (cf. (14) and (16)). The expressions p(T α0 ), u a (T α0 ) are plugged into (2) from which the algebraic relations (12) can be read off by applying (4). We note that the equations (11) -(17) apply for any time-orthogonal coordinate system.

Examples
The iterative scheme is implemented in MAPLE to calculate derivatives and consecutively the expansion coefficients.

Longitudinal expansion
Let us consider first a purely longitudinal expansion. Following [30] we employ for the initial conditions at Θ 0 = 0 a Gaussian energy distribution, e(Θ 0 , η) = e 0 exp{− η 2 2σ 2 } with σ = 3.8 for LHC energies and a synchronized flow y(Θ 0 , η) = η. The latter one is important for the finding, already anticipated in [30], that the energy density evolution can be approximated by e(Θ, η) ≈ e(Θ 0 , η) exp{− 4 3 Θ} and the flow velocity is only modified on the 2% level by a time dependent deviation from the Bjorken flow pattern characterized by a(Θ) according to y = η(1 + a(Θ)). This statement is made more explicit below. Figure 1 exhibits e(Θ, η)/e(Θ 0 , η = 0) (a), e(Θ, η)/ e(Θ 0 , η = 0) exp{ 4 3 Θ} (b) and y(Θ, η) − η (c) as a function of η for various time instants. The iterative scheme has been used up to order N = 48. In figs. 1 (a) and (c), the iterative results (curves) are compared to a numerical solution of the hydrodynamical equations (cf. eqs. (10) and (11) in [30]) obtained with a modified forward integration scheme (symbols). One observes agreement of both solution methods up to Θ = 3.2, i.e. τ = τ 0 e Θ = 24.5 τ 0 . At later proper times τ , the truncated iterative result shows some instability types, signaled by irregularities in e(η) accompanied by negative energy densities at large values of η. However, these can be separated from the physically interesting region. For larger values of Θ, e.g. Θ ≈ 4.0 (i.e. τ ≈ 55 τ 0 ), these instabilities even influence the energy density at the center. But at this proper time, the medium has cooled down so far that hadronisation is almost complete. In fig. 1 (c), the value of a(Θ) can be read off as the slope of the appropriate curve. Even for the largest depicted time Θ = 3.2 the slope ∆(y − η)/∆η being equal to a(Θ) is found to be small: 0.25/12 ≈ 0.02. Since the slope of y(η) increases monotonously with Θ, the absolute value of a(Θ) for 0 ≤ Θ ≤ 3.2 is always smaller than the above mentioned 2%. The conclusion for this simplified flow pattern is that the Bjorken type scaling of the energy density ∝ exp{− 4 3 Θ} and synchronized flow y − η ≈ 0 are excellent approximations for mid-rapidity |η| < 0.9, accessible in the ALICE experiment. Analytical corrections to the Bjorken behavior that can be obtained utilizing the iterative scheme are discussed below. Taking the initial temperature as T (Θ 0 , η ≈ 0) ≈ 430 MeV corresponds to the energy density e(Θ 0 , η ≈ 0) ∼ 65 GeV fm −3 (obtained with the s95p-PCE equation of state in [31]; for lattice data see e.g. [32]). Following the evolution until freeze-out at T f.o. ∼ 165 MeV, corresponding to the energy density of e f.o. ∼ 500 MeV fm −3 , requires to go up to Θ = ln(τ /τ 0 ) = ln(e(Θ 0 )/e f.o. ) 3/4 ≈ 3.7. Figure 1 demonstrates that the iterative scheme, truncated at N = 48 is stable in this time interval. To expose the importance of the initially synchronized Bjorken type flow pattern we chose y(Θ 0 , η) = (1 + a(Θ 0 ))η. For values of a(Θ 0 ) not to far away from the Bjorken flow (i.e. a(Θ 0 ) = −1 . . . 0.5), the deviation of the energy density from the reference system (a system with initially synchronized flow, i.e. a(Θ 0 ) = 0) can be parametrized by the series Here, e denotes the energy density of the system with initially not synchronized flow, and e ref is the energy density of a reference system with the same initial energy distribution, but Bjorken type initial flow. The parameters in this expansion can be directly linked to the Taylor coefficients of the energy density, by differentiating (18) and evaluating the result at Θ 0 . Apart from constant factors, the left hand side is the difference of the Taylor coefficients of e and e ref , which are calculable within the iterative scheme. The right hand side gives algebraic expressions for the parameters d n and b n . This method is quite powerful and can be used to transform the iterative solution given in terms of Taylor coefficients into an equivalent solution in terms of a different set of basis functions than powers of Θ (e.g. those given in (18)). For a system with non synchronized flow and initially Gaussian energy distribution, it turns out that the first non-vanishing terms already give excellent agreement: For general values of η, the expressions are quite cumbersome; at mid-rapidity they reduce to with a ≡ a(Θ 0 ) for ease of the notation. Figure 2 illustrates this approximation. One observes that at mid-rapidity a larger value of a(Θ 0 ) (i.e. a larger initial flow) leads to negative corrections with respect to an initially synchronized flow, meaning a more rapid decrease of central energy density. At large times one sees the difference approaching zero, which is a consequence of the energy density approaching zero. But if the initial flow is realistic (i.e. close to synchronized flow, a(Θ 0 ) ≈ 0), the difference decreases even  faster (i.e. with a larger decay rate in the exponential) than the energy density itself (which decreases approximately like the Bjorken energy density e BJ = e 0 e −4Θ/3 ). That means that the solution stabilizes around the solution with initially synchronized flow, and small perturbations have little effect on the long-time behavior. Since for realistic initial conditions (a(Θ 0 ) close to zero and σ ≈ 4 for LHC) already the leading-order correction is very good for all times until freeze-out, one can calculate the time and the value of the largest deviation from the reference solution. The maximum of e − e ref ≈ d 1 Θe −b1Θ is at Θ = 1/b 1 , where the deviation has the value of (d 1 /b 1 )e −1 . A similar analysis can be done to study the dependence of the energy density on the width of the initial energy distribution. It is found that the deviation of the energy density e obtained with a Gaussian shaped (standard deviation: σ) initial energy density from the Bjorken energy density e BJ can also be described with the series (18). In this case, the leading terms are The normalisation to the initial energy density is necessary to disentangle effects originating from gradients from the effect of different initial values at different spatial coordinates. Since the expressions for d 4 and b 4 at a generic η are quite extensive this coefficients are only given at η = 0: For LHC conditions the leading order term is much more important than the subleading order term (except at η ≈ √ 2σ, where it vanishes), i.e. the leading term gives already good agreement. In this case, the position and the value of the largest deviation e/e(Θ 0 ) − e ref can be calculated analytically, and one finds that the maximum value of (4d 2 /b 2 2 )e −2 is taken at Θ = 2/b 2 . In Fig. 3, the approximations are plotted at several values of η. One notes that the next-to-leading order approximation gives excellent agreement with the numerical data and the leading-order approximation is also good.

Transverse expansion: axisymmetric flow
Equipped with the finding of the previous subsection we now turn to the axisymmetrical transverse expansion. In doing so, we assume for the longitudinal flow y(Θ, η) = η consistent with e(Θ, η, r) = e(Θ, r) i.e. longitudinal scaling symmetry, which implies an infinitely large extension in longitudinal direction. The initial transverse energy density distribution is e(Θ 0 , r) = e 0 exp{−r 2 /(2κ 2 )} with the transverse coordinate r = x ⊥ 2 and κ being the standard deviation of the distribution. This transverse distribution is certainly to smooth and therefore underestimates the transverse gradients. In this respect it may give some lower limit on the transverse flow which is initially zero. Figure 4 shows the evolution of the scaled energy density and the transversal flow rapidity for a distribution with initial width κ = 7.5 fm corresponding roughly to the radius of lead nuclei used for the heavy-ion experiments at the LHC. The transversal expansion has a much larger effect on the dynamics than it was the case for the longitudinal expansion, since it is not dominated by the initial flow. For instance, at mid-rapidity, the scaled energy density drops by ≈ 5 % for Θ = ln(τ /τ 0 ) = 2.5 in the case of purely longitudinal expansion(cf. green dashed curve in Fig. 5), while for the transverse expansion, the scaled energy density is reduced by 58 % at the same time (cf. blue dotted curve in Fig. 5). The backbending of the rapidity curve in Fig. 4 (b) for Θ = 2.6 may be attributed to instabilities of the truncated iterative scheme for larger times. However, the energy density at the corresponding space-time region is already below the freeze-out energy density.

Transverse expansion: asymmetric flow
Transport properties of the sQGP can be extracted from the evolution of initial asymmetries. The most intensively studied quantity that encodes such information is the elliptic flow v 2 , which is used to measure, e.g., the shear viscosity of the sQGP. Therefore it is interesting to check whether the iterative scheme can cope with asymmetric initial conditions and if the expected change of eccentricity can be observed in the results. For this purpose, initial conditions e(Θ 0 , x 2 , x 3 ) = e 0 exp{−x 2 2 /(2σ 2 2 ) − x 2 3 /(2σ 2 3 )} with different widths σ 2 = 7.5 fm and σ 3 = 3.75 fm, respectively, along the two transversal axes were used. The choice of σ 2 : σ 3 = 2 : 1 corresponds to an impact parameter of b = 6/5R ≈ 9 fm for lead nuclei (in the hard sphere approximation). As in the previous example the initial flow is purely longitudinal with y(Θ 0 ) = η, and the longitudinal scaling symmetry is kept during the expansion implying independence of the energy distribution on the longitudinal coordinate η. In Fig. 6, the resulting eccentricity with e(x 2 , x 3 ) from a calculation up to order N = 43 is plotted (solid red curve).

Triaxial expansion
Even triaxial expansion can be treated within the iterative scheme. Since quantities like the elliptic flow are measured at mid-rapidity and since in the vicinity of η = 0 the energy density does not change much in the longitudinal direction, no significant change is expected for the eccentricities at η = 0. The dashed green curve in Fig. 6 is obtained from a calculation with triaxial initial energy density distribution e(Θ 0 , η, x 2 , x 3 ) = e 0 exp{−η 2 /(2σ 2 ) − x 2 2 /(2σ 2 2 ) − x 2 3 /(2σ 2 3 )} with σ = 3.8, σ 2 = 7.5 fm and σ 3 = 3.75 fm and initially synchronized Bjorken flow, i.e. y(Θ 0 , η) = η and v ⊥ = 0. It follows nicely the curve obtained for longitudinal scaling expansion, which shows that the spatial eccentricity and thereby the elliptic flow at mid-rapidity is not sensitive to the details of the longitudinal flow. However, at Θ 1.9 instabilities due to truncation make a comparison difficult. Curing this either by calculating at higher truncation order, optimizing the code or switching to an expansion similar to (18) remains to be done and is left to further studies.

Limitations
If it would be possible to calculate the Taylor coefficients up to infinite order, the solution would be exact, at least for analytical initial conditions and a smooth equation of state. However, due to limited computation power the iteration must be truncated at finite order N . As a result, only a truncated series is obtained and errors occur that are at least one order higher than the truncation order. This is not necessarily a problem, since in case of heavy-ion collisions one is only interested in a result covering a limited time span, for which such corrections may be small. In Fig. 5, the energy densities at the center of the fireball for the examples discussed above are plotted together with the freeze-out region for typical LHC initial conditions. One observes two facts. First, the central energy density of the purely longitudinal expanding system does not differ much from the Bjorken solution. This is seen already in Fig. 1 (b). Second, the inclusion of transversal expansion leads to a much faster cooling of the medium since it expands in more directions. Figure 5 shows that the iterative scheme gives reasonable results for the energy density at the center of the fireball. However, more problematic than the center are the regions where the gradients are large. To make this statement explicit, in Fig. 7, for each of the examples discussed above, the energy density at the edges is plotted. As representative for the edge, (one of) the point(s) is chosen, where the initial conditions have the largest gradients. One notes that the truncated iterative solutions get unstable quite close to the freeze-out region. This is an artefact of the truncated iterative scheme, which demands for calculation with a low freeze-out temperature of ∼ 130 MeV a higher truncation order. It is interesting to find out which properties of the initial condition improve the iterative result and which have a negative effect on it. Of course, the result gets better if one can calculate very high orders. To achieve this one needs to use initial conditions that can be differentiated most easily. For this reason, Gaussian energy distributions were used. An example for initial conditions that are not well suited for the iterative method are the ones deduced by Gubser [22]. To understand the computational limits, one can estimate the number n k of terms involved in the kth Taylor coefficient if the initial conditions are of Gaussian type. One observes that the Taylor coefficients have the structure of an exponential multiplied with a polynomial containing basically all spatial variables and all parameters in all powers up to k. This leads to the estimate, that n k grows like k m1+m2 with m 1 being the number of spatial variables and m 2 being the number of parameters. The total number of terms involved in the Taylor series truncated at order N is therefore ∼ N k=1 k m1+m2 ∼ N m1+m2+1 . A typical operation in the iterative evaluation is the multiplication of two of such series. For such a multiplication every term in one series needs to be multiplied with one in the other, leading to N 2(m1+m2+1) terms before simplification. For the transversal asymmetric system discussed in section 4.3 this results in a number of terms which is in the order of 43 2·3 ≈ 6 × 10 9 . Therefore, considering additional dimensions or letting parameters general leads to a rapidly increasing demand of computation power and reduces the calculable order dramatically. A second way to obtain results that are close to the complete solution is looking for initial conditions that lead to a quickly converging Taylor series. In this case, less orders need to be calculated, which obviously helps the computation. Testing different initial shapes of the energy profile we realized that the convergence is better if the initial conditions are more flat. Steep initial profiles lead to rapid oscillations in iterative results which can only be pushed to larger times by computing considerably higher orders. A typical example for this behavior is shown in fig. 8. The kink that denotes the point, where instabilities of the truncated scheme modify the energy density at the central point, occurs the later the broader the initial Gaussian is.

Conclusions
In this paper we present a scheme that can be used to construct solutions of the equations of relativistic fluid dynamics for differentiable initial conditions. The scheme is based on a formal Taylor expansion of the components of the energy-momentum tensor which is used to reformulate the time evolution equations as an infinite set of algebraic relations between the Taylor coefficients that can be solved iteratively. In the context of AdS/CFT similar methods are used to obtain solutions of the Einstein equations close to the boundary of (asymptotic) anti-de Sitter space [27,28], but to our best knowledge it was not applied to hydrodynamics or heavy-ion collisions in the literature so far. The scheme is utilized to address various cases of physical interest. It was shown that for the case of purely longitudinal expansion the results obtained with the iterative scheme are in excellent agreement with the numerical solution obtained with a standard numerical integration routine. One can parametrize the difference between a solution and a reference solution in terms of a special set of basis functions {d n Θ n e −bnΘ }. The parameters are related to Taylor coefficients obtained within the iterative scheme and it was found that already the leading order terms give very good agreement with the numerical solution obtained with a standard integration scheme. For transversal expansion it is possible to find solutions for the axial symmetric as well as for the axial asymmetric case. In the latter case, the change of the shape of the energy density due to elliptic flow has been considered. Even triaxial expansion needed for more realistic simulations is accessible within the iterative scheme. At mid-rapidity the longitudinal shape of the medium has almost no influence on the transversal eccentricity, but at finite η it is expected to have an effect.
Concluding we can state that the iterative scheme is a promising tool to investigate the properties of the sQGP produced in heavy-ion collision at the LHC via signatures built up during the hydrodynamic part of its evolution. With the computational power available for a personal computer at present, its application to fluctuating initial conditions or other systems with initially large spatial gradients is somewhat limited. Nevertheless there are several applications thinkable, e.g. the calculation of the background field in which a jet is moving (cf. the use of hydrodynamics in [34]) or the emission of real and virtual photons as in [35]. Even an application as a benchmark is possible, since the result is exact up to the truncation order and should therefore match the results of other hydro-codes.
A Energy-momentum tensor with derivatives Several systems are characterized by an energy-momentum tensor containing derivatives of the basic physical degrees of freedom. One important example is viscous relativistic hydrodynamics, e.g. formulated in the Landau-Lifschitz frame. For such an energy-momentum tensor it is not possible to express all components as functions of T α0 only. One possible way to cope with that problem is to allow also derivatives of T α0 as arguments, i.e.
with ; denoting the covariant derivative. From (30) one can follow the steps outlined in section 2. Another possible way would be to construct a similar iterative scheme directly for the basic physical degrees of freedom, i.e. in the case of viscous relativistic hydrodynamics for the energy density and the spatial components of the four-velocity.

B Faà di Brunos formula
The Taylor expansion (7) contains time derivatives of f µν , evaluated at the initial time. Since f µν does not depend explicitly on the coordinates the only time dependence originates from the time dependence of its arguments T α0 . With Faà di Brunos formula [36] it can be decomposed into derivatives of f µν with respect to T α0 and derivatives of T α0 with respect to θ: D n (f • g) = n! k 1 ! · · · · · k n ! (D k1+···+kn f ) • g n m=1 D m g m! km . (k1,...,kn)∈Pn To keep the notation as simple as possible, f and g are taken to be scalar functions, but the formula can be generalized to g being a vector and f being a tensor function. The set P n denotes all n-tuples of nonnegative integers which fulfill 1 · k 1 + · · · + n · k n = n and D is the differentiation operator. Since m is less or equal n, the nth time derivative of f • g contains only derivatives of g up to order n. Analogously the nth time derivative of f µν (T α0 ) contains only time derivatives of T α0 up to order n.