Holographic discommensurations

When the system with internal tendency to a spontaneous formation of a spatially periodic state is brought in contact with the external explicit periodic potential, the interesting phenomenon of commensurate lock in can be observed. In case when the explicit potential is strong enough and its period is close to the period of the spontaneous structure, the latter is forced to assume the periodicity of the former and the commensurate state becomes a thermodynamically preferred one. If instead the two periods are significantly different, the incommensurate state is formed. It is characterized by a finite density of solitonic objects -- discommensurations -- on top of the commensurate state. In this note I study the properties of discommensurations in holographic model with inhomogeneous translational symmetry breaking and explain how one can understand the commensurate/incommensurate phase transition as a proliferation of these solitons. Some useful numerical techniques are discussed in the Appendix.


Introduction
Different kinds of spatially modulated structures which break translation symmetry either explicitly or spontaneously are abundant in condensed matter system. It all starts from the crystal lattice, which forms a basis for any theoretical description of the condensed matter and breaks translations and rotations down to the discrete groups associated with Bloch momentum and crystal symmetry. On top of that some most interesting systems, including high temperature superconductors, demonstrate the spontaneous formation of the superstructures in the form of charge and spin density waves. Clearly, the interplay between these explicit and spontaneous mechanisms of translation symmetry breaking is of interest.
In holographic models of condensed matter systems (AdS/CMT) the status of translation symmetry breaking is different. This is an additional ingredient which one has to introduce on top of the Lorentz invariant theory of gravity. Firstly, the spontaneous translation symmetry breaking has been considered in [1][2][3][4][5][6] and later on the explicit potentials have also been introduced [7][8][9][10][11]. Only recently the interplay between these different mechanisms attracted some attention in the context of pinning of the spontaneous superstructure, pseudo-Goldstone modes and phonons [12][13][14][15][16][17][18]. But more importantly for the present study, the commensurate lock in between explicit and spontaneous structures have been studied in [19][20][21]. In [21] it has been shown that the commensurate state in holographic model provides a description to the Mott insulator, and the doping can be understood as departure from the commensurate to incommensurate or higher order commensurate state. In this note I will study in detail the mechanism behind commensurate/incommensurate phase transition in the holographic setup of [21] The paper is organized as follows. In the rest of the Introduction I'll discuss a toy model, which will help to set up the necessary notions and intuition. In Sec.2 the holographic setup will be introduced and the commensurate state will be discussed. Sec.3 is devoted to the discussion of the peculiarities in obtaining the nonlinear incommensurate solutions. In Sec.4 I will study the features of the discommensurations and address their role in commensurate/incommensurate phase transition. Two Appendices describe the numerical techniques and precision control.
The effect of commensurate lock in between two periodic structures is well known in physics [22][23][24][25]. The simplest system which demonstrates it is the Frenkel-Kontorova model, which considers a set of point particles, connected with springs, lying on top of the periodic lattice potential. In case the potential is absent, the springs assume the normal state and the particles form a periodic structure with spontaneous period λ p 0 . Clearly if one turns on the potential with exactly the same period λ k = λ p 0 , the particles will just fall in the minima of the potential and the springs will not be deformed. The state where the resulting periods of the explicit and spontaneous structures are equal λ k λp = 1 is called the lowest order commensurate. What will happen if the potential has a different period λ k > λ p 0 ? Then there are several possibilities. If the potential is strong enough, then the springs will be forced to stretch, in order that all particles fall into minima, and assume the modified period λ p = λ k (Fig.1, bottom left). This state is commensurate and the spontaneous structure is called commensurately locked in by the lattice. If instead the springs are very strong, and potential relatively weak, then some particles will acquire the additional potential energy, keeping the springs from stretching (Fig. 1, top left). This is incommensurate state since the resulting period of the spontaneous structure has no relation to the lattice spacing λ p = λ k . This incommensurate state is characterized by the fact that at any large enough sample of the system there are more particles then the minima of the potential.
If the potential in the incommensurate state becomes stronger, more and more particles fall in the minima, locally deforming the springs. Eventually the mismatch in the number of particles gets localized in the single unit cell of the lattice (Fig.1, middle left). The resulting state looks like the commensurate state everywhere, except around a local defect where two particles fall in the same potential well. This defect is a discommensurations -the soliton on top of the commensurate state, which accounts for the mismatch of the periods of two structures. Clearly, the discommensuration can be defined only if the state is not too far from commensurate point: λ p /λ k − 1 1, 1 otherwise the density of them gets so high that one can not  Figure 1. Cartoon of discommensurations. Left panel: λ k > λ p 0 , there are more particles then potential minima, positive discommensuration is formed. Right panel: λ k < λ p 0 , fewer particles then minima, negative discommensuration is formed. Note in both cases the discommensuration changes the staggered "red-blue" order in comparison with the "parent" commensurate state.
define a "parent" state. If the particles are charged, then the discommensuration, having one excess ball in the unit cell, bear positive charge. In complete analogy one can consider the negative discommensuration, which arises when λ k < λ p 0 and there are fewer particles then the potential wells ( Fig.1). This one is seen as one empty potential well and bears a unit negative charge. Interestingly, if we now introduce the additional Z 2 quantum number and assume that the dynamical system tends to form "anti-ferromagnetic" order, alternating the charge of the neighboring particles, then both types of discommensurations will act as a domain walls in this staggered order.
In what follows I will study the similar features of the discommensurations, which arise in the holographic model with spontaneous and explicit periodic structures.

The holographic setup
Consider the model of [26], used in [21]. It includes dynamical gravity with negative cosmological constant in 3+1 dimensions, U (1) gauge field and a pseudoscalar field, axion, which is coupled to the θ-term. The action reads [5,27]: here.
Here F = dA is the field strength of the U (1) gauge field A. Following [3,5,21,26,27], the couplings are chosen to be 2) Note that in these conventions the cosmological constant is Λ = −3 and the mass of the scalar is m 2 = −2.
It was shown in [3,5,26,27] that due to the ϑ-term this model develops an instability at low temperature, which evolves into the spatially modulated ground state, which breaks translations spontaneously and has a periodic charge density (CDW) with wavelength λ p 0 = 2π/p 0 2 , where p 0 (T ) is a temperature dependent thermodynamically preferred momentum. This is in complete analogy to the system with springs discussed in the Introduction. Noteworthy, this inhomogeneous state features the oscillating diamagnetic currents J y ∼ cos(xp 0 /2) on the boundary. The features of this state depend on the value of the coupling c 1 in (2.2) and I'll focus on c 1 = 17 used in [21].
Furthermore, one can introduce the background lattice potential, which would break translations explicitly. Following [8,28,29] it can be done by turning on a spatially modulated chemical potential, Featuring two spatial modulation scales: p 0 and k, this setup is sufficient to address the interesting physics of commensurability. At large temperature in absence of the explicit potential the ground state of (2.1) is the translational invariant Reissner-Nordström (RN) black hole The conformal boundary is located at z = 0 while the black hole horizon is at z = 1. Without loss of generality, one can set µ 0 =μ. I will express the dimensionful parameters of the model, denoted up until now in bold script, in units ofμ by making the replacements It can be shown that the ansatz with all unknown functions dependent on the holographic coordinate z and the boundary coordinate x, is sufficient to obtain the spatially modulated solutions of interest.
Following the standard AdS/CFT prescription, the boundary values of the holographic fields are dual to the one-point functions in the boundary theory: where µ is a spatially modulated chemical potential, ρ is a charge density and J ydiamagnetic current. We search for finite temperature solutions, which means that in the near horizon all functions must be regular. In addition, the equations of motion require Q tt (1, x) = Q zz (1, x), which in turn implies that the surface gravity is constant and given by (2.6), see e.g. [8].
It is relatively easy to construct the commensurate state. One starts from the spontaneous solution at certain finite temperature T 0 (in this note I will consider fixed T 0 = 0.1) with thermodynamically preferred momentum p = p 0 (T 0 ) and then turn on the explicit lattice with exact same momentum k = p. Then the two structures will be commensurate by construction. And one can use the standard DeTurck method in the periodic computation domain with period λ = 4π/k 3 . As it has been shown in detail in [21], this state is the holographic analogue of the Mott insulator. This will be the "parent" commensurate state, mentioned in the Introduction. The features of this 1/1 state are mostly similar to those of the pure spontaneous crystal: the staggered diamagnetic currents are seen and the charge modulation is present as well [21].
The focus of this note is thermodynamic stability of this state and transitions to the incommensurate one. But in order to study thermodynamic stability, we need to construct the competing incommensurate solutions first and compare the spatially averaged thermodynamic potential, given by where is the energy density, and s -the entropy density.

Full backreacted solutions
In order to study the thermodynamically preferred phases one has to construct fully backreacted nonlinear solutions corresponding to the coexisting ionic lattice and spontaneous crystal structures. There is a peculiar technical difficulty, which arises as soon as one addresses the nonlinear solutions. In [20] the spontaneous striped instability was considered in the perturbative regime. At the linear order we were able to introduce the continuous "Bloch momentum", characterizing the spontaneous structure which had not be proportional to the lattice period in any way. In case of finite amplitude of the striped structure this can not be done any more. The technical reason is that the e ipx multipliers can not be factored out from the nonlinear equations of motion and one has to rely on the position space representation, where the spontaneous structure is characterized by a certain period λ p , which is not necessarily proportional to the period of the ionic lattice λ k . In order to set up the numerical PDE solver procedure one has to specify only one scale, corresponding to the the size of the computation domain with periodic boundary conditions. At this point we see, that in practice we can only access the values λ p which are a rational multiples of λ k : In this case one can choose the computation domain of the size N k λ k equal to the integer number of lattice periods, which would simultaneously accommodate N p periods of stripes. 4 We see here that the accessible range of spontaneous structure wave-vectors k is now discrete and the density of them is limited by the maximal size of the computation domain, which we can handle in our numerical analysis. In what follows I will use the computation domains including up to N k = 20 periods of the lattice or up to N p = 20 periods of the spontaneous CDW 5 , which allows to achieve reasonably dense mesh in our plots for the thermodynamic potential of these solutions (see Fig.2). One might object that this technicality immediately renders it impossible to access the true incommensurate solutions in the mathematical sense, where λ p /λ k must be irrational number. But I would stress that due to the density of the rational numbers, in physics there is no way to distinguish between high order rational and irrational number. The practical definition of the commensurate (and higher order commensurate) state in this case will be the state where both numbers N p and N k are small integers. To certain extend in reality the maximum value of N k is restricted by the quality of the crystal, i.e. the size of the patch of a crystal lattice without defects, or the correlation length of the spatial order. In this regard having N k ≈ 20 gives us a reasonable approximation to physically incommensurate numbers.

Ω-Ωlat
Np/Nk At this point it is worth mentioning that because the striped structure is spontaneous, unlike the externally sourced lattice, one does not have a direct control over it's shape, including the number of periods and the relative phase (shift) between the stripe and the lattice. In principle the energetically preferred values are set by the dynamics and should be achieved by following the trajectories prescribed by the equations of motion. This is true for the relative shift, which is a continuous parameter. We observe that if we start computation with the different seeds corresponding to different alignments of the stripe and lattice structures the numerical procedure always converges to the same solution, corresponding to the preferred value of the relative shift. But due to the above mentioned effect, the set of accessible values of the wave-vector is discrete and the system can not smoothly propagate between them. By choosing the initial seed solution with a given periodicity, regarding the numerical relaxation procedure (see Appendix A) as adiabatic process, we ensure that the system will converge to the state with prescribed number of the spontaneous periods. This state though is only a local minimum of the thermodynamic potential and may be a false vacuum, so one has to construct the solutions with all possible number of spontaneous periods and compare their thermodynamic potentials in order to find the true ground state. See Fig.2.
Practically, in order to construct the solution with N p CDW periods on top of the N k lattice cells with period λ k and amplitude A we first find the spontaneous stripe solution with specific period λ p from (3.1) on top of the translationary symmetric background. Then we catenate N p copies of these stripes fitting them in the enlarged calculation domain. At this point we turn on N k periods of the background lattice by slowly changing the boundary condition for the chemical potential, eventually achieving the desired value of A in (2.3). This adiabatic process preserves the initial number of the CDW periods what we check numerically at every stage by counting the number of zeros of the oscillating A y field at the horizon (see Fig.3).
In complete analogy with the perturbative study of [20], in order to explore the phase diagram at given temperature we first choose the period and the amplitude of the explicit ionic lattice. Then we construct a set of nonlinear solutions, corresponding to the spontaneous structures with different wave-vectors λ p on top of this lattice. We calculate the thermodynamic potential (2.14) for these solutions w(λ k ) and we find the one which is thermodynamically preferred. The sample of the w(λ k ) curves, which we get, is shown on Fig. 2. As an example throughout this note I will use two values of the lattice momentum: k = 1.06 and k = 1.53. At the temperature under consideration T = 0.01 the spontaneous momentum of the CDW is p 0 (T ) = 1.28. Therefore the former lattice has longer wavelength then the CDW, promoting positive discommensurations and the latte has shorter wavelength, promoting negative discommensurations.
There are several features on Fig. 2, worth noticing. Firstly, at A = 0 the solutions follow the curve, which one would obtain in the study the spontaneous striped solutions on the homogeneous RN background [3,5,27]. I've checked that for c 1 = 9.9 the results coincide with Fig. 2 in [27]), which is a valuable check of the numerical method.
As the amplitude of the lattice is increased, the Ω(λ k ) curves start to deviate smoothly from the RN case. Interestingly, on both of the plots "one half" of the curve at finite A is missing. In case k < p 0 it is the part of the curve with 2πλ −1 p < k, in case k > p 0 -the one with 2πλ −1 p > k. The technical reason is that as we turn on finite amplitude of the chemical potential for a prescribed number of stripes in the seed solution, the numerical procedure does not converge, or converges to the solution with different number of stripes. More precisely, if for k = 1.06 (left panel) we choose a seed with 10 CDWs and turn on 11 periods of the lattice (this would correspond to the point on the "missing" left shoulder N p /N k = 10/11 <! on the plot), we will see that the resulting solution has 10 CDWs (now sitting on the right shoulder N p /N k = 12/11 > 1). Effectively, one period of the CDW is created dynamically and we can not construct the desired solution with N p = 10. In principle this kind of behavior is allowed, given that the stripes are completely spontaneous, and the phenomenon is quite robust: we observe the disappearance of one half of the curve everywhere in the parameter space, which we study. I will address the physical reason behind it in the following Sections.
As one can see, even though we have access only to the discrete set of values, they lie on the smooth curves which have a well defined minima. There are two distinct possibilities, the minimum thermodynamic potential is achieved for the period of the stripe close to it's spontaneous value λ p 0 , or for the period which is proportional to the lattice spacing 6 . The former possibility defines the incommensurate state, the latter -commensurate lock in. One can see that as the amplitude rises the minimum smoothly shifts from the incommensurate to commensurate point. Henceforth by rising the amplitude we observe the smooth, at least second order phase transition.

Discommensurations
Let's now consider the incommensurate state. As I mentioned earlier, the numerical computation in this case is technically more involved, as the numbers of periods in (3.1) can become large. It is instructive to start the discussion by focusing on the solution which is closest to commensurate Np N k = 1 1 value. Given that I will choose N k = N p − 1 and maximal N p reachable by my numerics N p <= 20.
As we learned in the previous sections, in the commensurate state one finds one period of the lattice potential per one period of the spontaneous CDW. One can say that the incommensurate solution with 20 CDW's per 19 lattice periods would have exactly one excess period of spontaneous CDW structure per 19 unit cells as compared to the commensurate state on top of the same lattice. By inspecting this solution (Fig.3) we see, that the solution profile coincides with the commensurate state almost everywhere except from the finite size region in the core, where this excess of 1 period of the spontaneous structure is accounted for. We can also study the TD potential and charge density of such solutions as compared to the pure commensurate "parent" state, Fig.4, which shows clearly that this  The excess/lack of one period of the spontaneous structure is clearly seen near horizon (top). Outside the core of a disc. the profile of the solution is identical to that of the commensurate state.
incommensurate solution can be seen as a commensurate state with one localized discommensuration(disc.) on top of it. Similarly, the solution with 18 CDW's per 19 unit cells includes a single discommensuration with deficiency of 1 CDW. As seen on Fig.5, the size of a disc. does indeed decrease when the lattice gets stronger as anticipated from our toy model. In this point of view the discommensuration is a soliton on the commensurate background with a topological charge ±1, coinciding with the number of missing/excess periods of the spontaneous charge modulation in the domain. We observe both of these types of discommensurations in our model. The positive disc. appears when one considers wavelengths of the ionic lattice, larger then the wavelength of the spontaneous crystal λ k > λ RN p , and the preferred commensurate fractions N p /N k are larger then 1. The negative discommensuration is seen when N p /N k < 1. Both types are direct analogues of discommensuration studied in the context of charge density waves in [23].
From Fig.4 it is apparent that the disc. carry a manifestly positive (negative) electric charge. The important difference between the holographic discommensuration and our trivial example is that its charge is not fixed and varies smoothly, see Fig.6. To complete the comparison of the holographic discommensuration with the toy model expectations, let's consider the current profile of these solutions, Fig.7. Negative discommensuration Figure 5. Size of the discommensuration, obtained as a width at 0.1 height of the charge density, depending on the amplitude of the lattice. As the lattice becomes stronger, the discommensuration localizes, in complete analogy with the toy model.
Clearly, disc. serves as a domain wall in the staggered current order. At this moment we can understand the absence of one shoulder in the w(λ p ) curves, observed in Sec. 3. Indeed, starting from the commensurate point one can move to higher stripe wave vectors by adding charge +1/2 disc., or to lower wave vectors by adding charge −1/2 disc. For any given set of background parameters only one of them is dynamically stable. This can be understood from the simple energy argument: for k > p 0 the CDW in commensurate state has higher wave vector then the spontaneous one, hence lowering the wave vector by negative charge disc. would lower the potential energy. From the other hand, the shape of the disc. as a localized object contributes to kinetic energy. The balance between these contributions would stabilize the soliton with charge −1. If now one would consider the positive charge  Figure 6. Charge of a single discommensuration measured with respect to the "parent" commensurate state, depending on the amplitude of the lattice. There is no preferred values, as opposed to the naive expectations.
Negative discommensuration disc. on the same background, one would see that it rises the stripe wave vector, hence the potential energy is also rising, while the contribution from the kinetic energy is always positive. Hence there is no way the different energy contributions can be balanced in the soliton and it is not stable.
One can see that the further deviation from the commensurate point 1, according to (4.1), is achieved by rising the density of discommensurations (disc.) of the particular type, i.e. considering one disc. per fewer lattice periods. In this regard we can re-analyze the data which we obtained in Sec. 3. Figure 8 shows the thermodynamic potential of the state, measured w.r.t. the thermodynamic potential of the commensurate state, as a function of the density (n) of disc.. While the density is low enough, and the separation between the solitons is much larger then their size, the apparent linear dependence of w(n) just follows from the fact that the isolated soliton has a certain fixed "mass", which is given by the slope of the w(n) curve at n → 0. The field profiles in this regime look as the domains of commensurate solutions separated by the evenly spaced lattice of discommensurations. It is also evident that the "mass" of a soliton depends on the parameters of the solutions and can be negative as well as positive, see Fig 9. The positive mass would mean that the creation of disc. costs energy and the commensurate state is therefore energetically stable. The negative mass, on the contrary, signals the instability of the commensurate state. The solitons start to proliferate and their density rises, forming the increasingly dense discommensuration lattices and driving the state further from the commensurate value of 1/1. As the density rises, the distance between the solitons becomes smaller and they start to interact. The repulsion between disc. limits the energetically preferred density from above and in this way the system assumes the incommensurate stripe wavelength. p Ω-ΩMott p Ω-ΩMott Figure 8. Thermodynamic potential as a function of the density of discommensurations, measured with respect to the commensurate state at different amplitudes of the lattice (same data as on Fig.2). For temperature T = 0.01 the spontaneous momentum is p 0 (T ) = 1.28. Left panel shows the case when the lattice momentum is smaller: k = 1.06 < p 0 , and positive discommensurations increase the total momentum of the incommensurate solution. Right panel -the lattice momentum is larger k = 1.53 > p 0 and negative discommensurations reduce the momentum of incommensurate solution. The slopes of the dashed lines correspond to the mass of a single discommensuraton. When this slope becomes negative, the commensurate state is unstable and the discommensurations proliferate.

Conclusion
In this note I studied in detail the new solitonic object, which appears in the holographic model as soon as one considers the strong interplay between explicit and spontaneous symmetry breaking. I demonstrated that in many regards these holographic discommensurations conform the naive expectations obtained from the classical toy model: they are localized object, which are responsible to the mismatch between the periods of the explicit and spontaneous structure, they carry charge and they realize the domain walls in the staggered order parameter. The important Negative discommensuration Figure 9. Mass of an isolated discommensuration, depending on the amplitude of the lattice. At small amplitude the mass is negative signaling that the commensurate state is unstable. At large amplitude discommensurations become massive and the pure commensurate state stabilizes.
difference though is the absence of the preferred value of charge, which would be quantized in the classical model. Discommensurations play an important role in the commensurate/incommensurate phase transition. The commensurate state is stable as long as the mass of disc. is positive. Once the strengths of the lattice is lowered, the mass goes negative and discommensurations proliferate, forming the incommensurate state.
The involved numerical analysis is needed to study discommensurations. Given that they are solitons on top of the nonlinearly constructed ground state, it would be interesting to find out whether any analytic control over them is possible. To some extend they are similar to the Abrikosov vortices in the superconducting condensate. The latter can be analyzed perturbatively near the critical temperature, when the condensate itself is small. It would be interesting to analyze the similar approach to disc. and by the Foundation for Research into Fundamental Matter (FOM) and partially by RFBR grant 15-02-02092a.
The numerical calculations have been performed on the Maris Cluster of Lorents Insititute.

A Numerical techniques
In the present study we rely heavily on numerical analysis of the holographic nonlinear solutions. Moreover, in order to study the phase diagram and cover the parameter space we have to obtain several thousands of the solutions to equations of motion, with some of them requiring quite large calculation grids in the spatial direction. This situation puts a very strict requirements on the numerical techniques, which we use and the precision and accuracy of the results. In this Appendix we review the key features of the numerical setup, used in the present study.
Roughly speaking, the process of the numerical solution of the system of nonlinear PDEs consists of a few key steps [31] For nonlinear problem one has to evaluate these coefficients anew at every one or a few steps. This operation requires as many evaluations as the number of nodes and thus deserves optimization.
3. At this moment the problem can be written as system of linear algebraic equations. One can solve it either exactly by computing the inverse of the linear operator matrix (this corresponds to the Newton-Raphson method), or approximately, using some kind of "relaxation" scheme: the various options include (preconditioned) Richardson relaxation, Gauss-Seidel iterations or ILU decomposition.
4. Once the inverse is computed, the increment in the functional variables can be evaluated, which is used to construct the next approximation to the solution.
5. The iterations continue until some criterion of the accuracy or precision of the current approximation is fulfilled.
In order to build an efficient solver we need to choose carefully our tactics at every step.
Starting with the discretization scheme of step 1 one usually chooses between pseudospectral collocation [32] and FDD (finite difference derivative) of a certain order. The advantage of pseudospectral scheme is its improved accuracy, hence one needs much less grid points in order to approximate the solution well. In some cases this method also demonstrates the exponential convergence. Pseudospectral discretization works perfectly in the periodic spatial direction, where the solution is smooth and is well approximated by a Fourier series. The drawback is the efficiency of the scheme on the finite interval of the holographic coordinate. Here one has to use the smooth Chebyshev polynomials, which are not suitable for approximation of the non-analytic behavior near the UV boundary and IR horizon. This mismatch thwarts the exponential convergence and may lead to the breakdown of the whole scheme. The FDD approach doesn't suffer from this drawback but requires much more grid points to reach the comparable accuracy. This has undesirable side effects on steps 2 and 3, when the number of evaluations of coefficients and the size of the linear matrix to be inversed are increased, correspondingly. It is worth mentioning though, that FDD scheme leads to much sparser linear matrix on the step 3. If one chooses to inverse the matrix exactly, this task is next to impossible for sizable grids in the pseudospectral case, where the differentiation matrices are dense. The compromise would be to use several patches along the holographic direction [28], using pseudospectral approximation in the interior and FDD near the boundaries. After experimenting with all three options, see below, we've chosen a single patch pseudospectral scheme in the holographic direction, which proved to be quite robust in practice.
We use Wolfram Mathematica [33] in order to implement the numerical algorithm.
It may have the disadvantage in speed, when one compares it to the lower level computing languages like C++ or FORTRAN, but as long as one uses the high level efficient precompiled routines like: LinearSolve, NDSolve'FiniteDifferenceDerivative, and sparse matrix operators the difference in speed becomes less obvious. The elementwise operations required at step 2 can be efficiently compiled with Compile, which brings up a spectacular acceleration. In the end of the day, the most important limitation of Mathematica is the necessity to work with MachinePrecision numbers in the compiled function, which eventually limits the precision of our results, as discussed below.
As we've discussed already, the direct inversion (Newton-Raphson method) in case of pseudospectral discretization is extremely demanding for the large grid which we use. Moreover, since we are solving the nonlinear problem, which inevitably requires an iterative procedure, obtaining very precise result for the matrix inverse at step 3 doesn't make much sense. That's why we use a relaxation scheme instead. The speed of convergence of a relaxation scheme is defined by the highest eigenvalues of the linear operator matrix. In order to make the process more efficient one uses the preconditioning, i.e. one multiplies the operator by the preconditioner matrix, which brings all the eigenvalues to the same scale. The ideal preconditioner is the inverse of the operator itself, but a reasonable approximation to it will also work fine. We use the differential operator evaluated in the low order FDD scheme, as a preconditioner. It approximates the highest eigenvalues very well and is relatively easy to invert, being very sparse. The result is a nonlinear Richardson relaxation with Orszag preconditioning (See Sec. 15.14 and eq.(15.115) in [31]). One can also view this approach as a pseudo-Newton method, where we use the approximation to Jakobian instead of the exact one. It should be stressed here, that even though we use a low order FDD scheme in the construction of the linear operator, the equations to be solved use the full pseudospectral approximation to the derivatives, thence we achieve the pseudospectral accuracy by iteratively inverting the sparse operator on a relatively coarse grid. The relaxation scheme requires, in principle, more iteration steps than the Newton-Raphson, so we effectively exchange the memory consumption with the CPU time. Altogether for the hardware which we used the relaxation procedure turned out to be an order of magnitude faster then the analogous Newton-Raphson scheme.
As we mentioned already, the high eigenvalues of the linear operator are well approximated using the low level FDD preconditioner. This is not true for the lower eigenvalues, which define the long-wavelength errors with slowest relaxation rate. One can fight these ones and further improve the efficiency of the numerical scheme by using multigrid technique [34]. In practice we found that without multigrid our solutions converge already after ∼ 20 steps. In this situation the overhead of transitioning between fine and coarse grids becomes relatively significant and we found no improvement of the overall efficiency in full multigrid method.
In the end of the day we managed to optimize the calculation scheme to the extent when it takes about half an hour to obtain the precise solution on our largest grid of size ∼ 330 x × 80 z (pseudospectral) using a single core of a laptop CPU (Intel Core i7-5600U @ 2.60GHz ) and about 3 Gb of RAM.

B Precision control
As one can see from our results, the difference between the free energies of the solution with spontaneous structure and without it is just of order of few percents of the free energies themselves. That means that in order to reliably study this difference, we need to evaluate the free energies with accuracy of at least 10 −4 . This puts a challenge to our numerical scheme and renders the strict precision control a necessity.
The accuracy of the numerical results depends substantially on the size of the grid. On one hand, it is clear that the denser grid delivers higher accuracy, i.e. the closer approximation of the real result of continuum limit. On the other hand, very dense grids bring up numerical precision issues, when the rounding errors play a dominant role in the calculation of the derivatives. Thus we have to optimize our grid size maximizing the accuracy while keeping the rounding errors below the certain value. This study also helps to choose between different methods of discretization in the holographic direction, discussed earlier.
In our 2-dimensional problem the grid size in two directions: N x and N z , has to be optimized. As an example we use a sample solution describing the commensurate stripe on top of the lattice with amplitude A = 1., lattice wave-vector k = 1.5 at temperature T = 0.01. We perform a relaxation of this solution on the set of grids with N x = {9, 17 Each time the iteration procedure is run until max |df | hits the MachinePrecision bound ∼ 10 −11 , while the values of the functions are of order one. Due to the numerical rounding errors the precision of the solution can not be improved after this bound is reached. By studying the value of dw/w at this point one can estimate the precision of the physical observables, obtained on a given grid. On Fig. 10 the relative precision of the thermodynamic potential is shown for various grid sizes in pseudospectral approach. One can see that for dense grids N y > 80, the precision is decreasing. No similar effect is seen for the dense N x grid. This important observation tells us that with given MachinePrecision there is no reason to use the grids with N y > 80 and also that the numerical precision of the thermodynamical potential, which we calculate is not be higher then 10 −8 .
▼ N x =129 Figure 10. Precision of the thermodynamic potential limited by the rounding errors for different Chebyshev grid sizes.
In principle, the accuracy of the numerical results must increase with increasing grid size, giving the exact match in continuum limit. In order to estimate the exact result w ∞ we evaluate the thermodynamic potential for a set of increasing grid sizes and then extrapolate to infinity. The accuracy for a given grid is then given by |w − w ∞ |/w ∞ . As one can see from Fig. 11, the optimal accuracy reached at N y = 80 is of order 10 −7 . As to the spacial grid resolution, one can see that already for N x = 33 the result is close enough to w ∞ , i.e. already at N x = 33 the accuracy is controlled by the holographic axis resolution N z . Wn-W∞ Figure 11. Accuracy of the thermodynamic potential depending on the size of grid in holographic (N y ) and spacial (N x ) direction.
In the end of the day we observe that for a single patch Chebyshev grid the maximum N y resolution is limited by the rounding errors at N y = 80. The accuracy of the thermodynamical potential for a grid of this size is about 10 −7 . We use this value as a numerical error estimate throughout the present study and it proves to be quite enough for our main results. We also see that the x-axis resolution can safely be set to N x = 33 points per one period of the stripe, without affecting the accuracy of result. We should note here that the comparable accuracy in FDD approach is reached for N y ≈ 320. Choosing the grid of this large size is disadvantageous on the other stages of calculation, so we choose rely on the pseudospectral approach. We've also checked that the patching technique doesn't bring any significant improvement of the accuracy.
One should keep in mind that in the numerical procedure we solve DeTurck equations, so it must be checked that the Einstein equations are satisfied, which is done by two independent measures: max |G µ µ | and max | ξ 2 |. For not very low temperatures T > 0.005 at the grids which we use these values are both of order 10 −7 , which is quite satisfactory. However, at lower temperatures they increase and the use of higher N y -resolution grid is necessary. As stated above, our N y -resolution is bounded by the rounding errors and MachinePrecision, henceforth with the present numerical scheme we can not reliably access extremely low temperatures. Nonetheless, this drawback has little impact on the results of the present study.