{'fraction_non_alphanumeric': 0.04063851335058755, 'fraction_numerical': 0.012490555725239924, 'mean_word_length': 4.390174161684429, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Phase transition into Instanton Crystal

We propose a class of models exhibiting instanton crystal phase. In this phase, the minimum of the free energy corresponds to a configuration with an imaginary-time-dependent order parameter in a form of a chain of alternating instantons and antiinstantons. The resulting characteristic feature of this state is that the average of the order parameter over the imaginary time vanishes. In order to study the model in a broad region of parameters of the model quantitatively, and prove the existence of the instanton crystal phase, we develop an efficient numerical scheme, suitable for the exact treatment of the proposed models. In a certain limit, results demonstrating the existence of the instanton crystal phase are obtained also analytically. The numerical study of the model shows that there is a phase transition between the instanton crystal and the state with the imaginary-time-independent order parameter.

We propose a class of models exhibiting instanton crystal phase. In this phase, the minimum of the free energy corresponds to a configuration with an imaginary-time-dependent order parameter in a form of a chain of alternating instantons and antiinstantons. The resulting characteristic feature of this state is that the average of the order parameter over the imaginary time vanishes. In order to study the model in a broad region of parameters of the model quantitatively, and prove the existence of the instanton crystal phase, we develop an efficient numerical scheme, suitable for the exact treatment of the proposed models. In a certain limit, results demonstrating the existence of the instanton crystal phase are obtained also analytically. The numerical study of the model shows that there is a phase transition between the instanton crystal and the state with the imaginarytime-independent order parameter. The standard way of describing a phase transition is based on the concept of an order parameter introduced by Landau [1]. This quantity equals zero in the disordered phase but is finite in the ordered one. The order parameter can be scalar, vector, tensor, etc. The beauty of this approach follows from the universality of the description because the critical behavior depends on the symmetry of the order parameter rather than on the details of the interaction.
Although the Landau theory is by construction applicable only near the critical point and not too close to it, so that the fluctuations can be considered small, the concept of the order parameter provides the means of description of the ordered phase for all the temperatures below the transition temperature. As such, in the case of Z 2 symmetry breaking, for example, the order parameter is real and the minimum of the free energy corresponds to the two possible values +1, −1 (if properly rescaled) of the order parameter.
The situation becomes more interesting in the quantum limit at low temperatures, where the tunnelling effects become important. The toy model which is usually employed to discuss this kind of phenomena is that of a particle moving in a double-well potential. The convenient way to study the thermodynamics is to use the Euclidean path integral formalism with the imaginary time τ . The corresponding Euclidean action of the toy model is introduced as where x is the coordinate of the moving particle. The lagrangian L can be interpreted as the classical lagrangian of the particle in the inverted potential One can solve for the minimum of the action (1.1) by putting its first variation to zero. The action is minimized by the trajectories x(τ ) ≡ corresponding to the particle sitting in either of the minima of the double-well potential. However, there are also additional imaginarytime-dependent solutions where τ 0 is an arbitrary time. These solutions are usually called 'instantons' or antiinstantons depending on the sign. They describe the classical trajectories connecting the two "humps" of the inverted potential. In the Hamiltonian language, this new solutions correspond to the tunnelling between the two "vacua" leading to the splitting of the ground state. Besides the single instanton or antiinstanton solutions, there are also the solutions consisting of multiple instantons and antiinstantons chained together. The action at stationary solutions x(τ ) ≡ ±1 equals zero. At the same time, its value at the instanton solutions is higher: S inst = 4/3. In the condensed matter setting, we should scale the action (1.1) by the volume of system V . As the result, the contribution of instanton solutions is exponentially small in the thermodynamic limit. We should note that this picture stays valid if we consider any potential with the shape similar to the one displayed in Fig. (1.1): there would be instanton solutions, however their action would be higher than that for the time-independent solutions.
It is interesting to note that the study of instanton physics has been pushed forward in the field of QCD (see e.g. [2]) in order to understand the structure of the "ground" state in that theory. The situation there is quite to similar to what we have just discussed, namely, there are many different minima of the Euclidean action corresponding to the different vacua of the theory, and there are instanton configurations connecting these vacua. However, the number of vacua is infinite, and the exponential supression of the instanton configurations caused by the increased action is compensated by the increased phase space factor. As a result, the "ground" state of the system is described by a non-trivial combination of instantons, which is referred to in the field of QCD as "instanton fluid".
In view of these interesting developments, it is natural to ask whether it is possible to formulate a model in the setting of Condensed Matter Physics, which would admit a thermodynamically stable state described by an order parameter consisting of a system of instantons and antiinstantons? Of course, one should expect that such a model should be consideraly more complicated than the simple toy model given by Eq. (1.1). Still, the question remains the same.
In this paper, we suggest a new model that allows us to obtain in some region of parameters the ground state with the imaginary-time-dependent order parameter. This order parameter can be visualized as a lattice of alternating instantons and antiinstantons, and therefore we coin for this phase the name "instanton crystal". The model contains both interacting fermions and boson modes. It does not contain any infinite or long-range interactions, thus this model is in all respects rather conventional for condensed matter physics. Using the meanfield approximation, we solve this model both analytically (in certain region of parameters) and numerically. The study reveals a competition between the phase described by the static time-independent order parameter and the instanton crystal phase with a transition between the two phases.
In spite of the popularity of the instanton physics in QCD, only few works have been done in the past on the investigation of a possibility of non-perturbative effects in imaginary-time representation in condensed matter physics. To the best of our knowledge, this problem has been first attacked in Refs.° [3][4][5][6] using a two-band model of interacting fermions. A solution with the chain of instantons and antiinstantons has been obtained. Unfortunately, it has finally been concluded that in this model the free energy for the imaginary-time-dependent configuration was always higher that for the static configuration. A similar model with a coordinate instead of the imaginary time τ had been used long ago as 1D mod-els of polymers, and a solution with the chain of kinks and antikinks (spatial analogue of instantons and antiinstantons) had been obtained. Ironically, it this case the energy of the "instantons" in the coordinate space could be lower than that for the homogeneous solution.
Recently, one more attempt has been undertaken [13] to obtain the instanton crystal. In that work, an additional interaction term has been added to the previous model of Refs. [3][4][5][6] with the hope that it might make the free energy of the instanton crystal lower than the static state. Indeed, the presence of this new term reduced the free energy but using the perturbation theory could not help proving the existence of the instanton crystal. Many other guesses remained just guesses due to the technical difficulties.
In the present paper, we modify the previous models by introducing bosonic current-like modes which are coupled to the fermions. In comparison to the previous work [13], we formulate an effective numerical scheme, which enables us to solve the mean-field equations in the general case. Moreover, in a certain region of parameters the analytical treatment of the model is feasible. Results of both the numerical and analytical study allow us to conclude that the instanton crystal can exist in the thermodynamic equilibrium. Actually, this work is the beginning of a systematic study of properties of the instanton crystal phase.
The paper is organized as follows. In Section II we introduce the model without discussing its origin. This is because we hope that the model is rather general. In Section III we minimize the effective Lagrangian and derive the mean-field equations. In Section IV, we solve the mean-field equations in a certain region of the parameters of the model and calculate the free energy, demonstrating the possibility of the instanton crystal phase. In Section V we describe the numerical scheme for the solution of mean-field equations, which we use in Section VI to conduct a detailed numerical investigation of the general case. In Section VII we discuss a possible origin of the model. Finally, in Section VIII, we discuss the results obtained and the perspectives of the future studies. The appendices contain technical details of the calculations.

II. GENERAL MODEL FOR THE
IMAGINARY-TIME CRYSTAL.
A. Hamiltonian of the model.
In this section, we formulate a rather general macroscopic model of interacting fermions and bosons without going into details of its possible origin. The latter will be done in Sec. VII but here we simply introduce the general HamiltonianĤ and discuss its structure.
The total HamiltonianĤ of the model consists of three In Eq. (2.1), HamiltonianĤ 0 stands for a system of noninteracting fermionŝ These fermions live in two bands 1 and 2. Fourcomponent vectors contain as components creation and destruction operators c s αp for the fermions from the bands s = 1, 2 with spin projections labeled by α = 1, 2 (actually, the spin variable α is not very important here). The vectors c + p are hermitian conjugated to c p and contain creation operators as components. The energies ε ± p are expressed in terms of the spectra ε 1,2 (p) in the bands 1, 2 as The operatorsΣ i , i = 1, 2, 3 are Pauli matrices acting in the subspace of the bands 1 and 2, whileǏ is the identity operator acting in the same subspace. The second termĤ int in Eq. (2.1) stands for the interaction between the fermions from different bandŝ where V is the volume of the system. The Hamiltonian H int contains contact attraction (first term) and repulsion (second term). Actually, the first term in Eq. (2.5) describes attraction of fermionic currents, while the second one stands for repulsion of charges. More information about the possible origin of the model and interpretation of the terms in the Hamiltonian is given in Sec. VII. It is worth emphasizing that the HamiltonianĤ int does not contain any long-range interactions. We should also note that a model with the HamiltonianĤ 0 +Ĥ int was considered previously in Ref. [13]) in a form adopted to a direct use of the mean-field theory. The third termĤ B describes a system of current-like modes labeled by different momenta q.
(2.6) WhereQ q andP q are conjugated coordinates and momenta of these modes satisfying the following relationŝ The modes could be, in principle, just phonons but the latter generate extremly low currents with respect to the fermionic ones (of the order m/M, where m is the electron mass, while M is of the order of atomic masses) that cannot efficiently interact with fermionic currents. The current-like modes are coupled to the vector po-tentialÂ q created by the fermions: For a typical electron-phonon interaction, the electrons are coupled to coordinates of phonons. In contrast, in our case, the fermions are coupled to the momenta of the modes. In other words, we include into consideration current-current interaction. This is an unusual feature, and it is crucial for our new results. However, we argue in Sec. VII that the existence of these modes and their interaction with the fermions is not unrealistic. We assume that all the coupling constants are not negative, namely As usual, in the limit of large volume V → +∞, one can replace the sum over the momenta by integrals using the standard replacement (d is dimension), which allows one to see thatĤ is proportional to the volume V , as it should be. In this paper we restrict ourselves to studying thermodynamic real-time-independent properties of the model. In order to study them, the partition function Z in grand canonical ensemble is introduced in the standard way where µ is chemical potential, and it is convenient in the following to absorb it into the definition of ε 1,2 (p). The HamiltonianĤ, Eqs. (2.1-2.8), describes a system of interacting fermions and bosons, and it does not contain any long-range interactions. We would also like to emphasize that neither the HamiltonianĤ nor the partition function Z contain any time-dependence, be it real time or imaginary time. As a consequence, there can be no doubts that the HamiltonianĤ describes a rather convential system in thermodynamic equilibrium.
In principle, one could proceed with the analysis of the model using the operator formalism. However, it is more convenient for our study to use the corresponding Lagrangian formulation based on rewriting the partition function in the form of the functional integral over commuting and anticommuting fields. This way, the imaginary time also enters the picture. At the same time, we found instructive to provide the explicit Hamiltonian of the system in operator formalism: it helps to stress the fact that we are considering a system in thermodynamic equilibrium without any pumping or relaxation. In the Lagrangian formulation, the partition function Z can be written in the form of a functional integral as Z =ˆexp −S χ, χ + , a DχDχ + Da. (2.12) where the action S [χ, χ + , a] contains anticommuting fermionic fields χ s α (τ ), χ s+ α (τ ) and commuting bosonic fields a q (τ ) and reads S χ, χ + , a = S 0 + S int + S B + S FB . (2.13) In Eq. (2.13), the terms S 0 and S int correspond to the termsĤ 0 andĤ int in the Hamiltonian approach, Eq. (2.1), respectively. At the same time, the terms S B and S FB correspond to the termĤ B . As usual, the imaginary time τ is defined for 0 ≤ τ ≤ β ≡ 1/T where T is the temperature. The fermionic fields χ s αp (τ ), χ s+ αp (τ ) obey standard antiperiodic boundary conditions 14) and have the structure identical to the vectors c p and c † p . In contrast, the bosonic fields obey the periodic boundary conditions a q (τ + β) = a q (τ ) . (2.15) The first term in Eq. (2.13), S 0 , is the action of noninteracting fermions The interaction term S int [χ, χ + ] has the form The term S B [a] in Eq.(2.13) stands for the action of the phonon-like modes where a q (τ ) are complex fields satisfying (a q (τ )) * = a −q (τ ). (2.19) The fields a q (τ ) correspond to the coordinates in the language of oscillator modes, and da q (τ ) /dτ correspond to their velocities. Finally, the coupling of the fermions to the current-like modes is described by the term S FB [χ, χ + ] in Eq. (2.13), which takes the form The terms S B [a] and S FB [χ, χ + ] constitute together the imaginary time Lagrangian corresponding to the Hamiltonian termĤ B . The Lagrangian formulation of the functional integral can be obtained writing the phasespace functional integral corresponding toĤ B and then integrating out the momenta.
In the model under consideration, the expression χ + p (τ )Σ 2 χ p+q (τ ) describes a current (see Section VII). Therefore, one can interprete the term S FB [χ, χ + , a] as the interaction of fermionic and bosonic currents.
In principle, one can integrate out in Eq. (2.12) either fermionic or bosonic fields just in the beginning of calculations. In order to compare the model described by Eqs. (2.12) and (2.13) with models studied in the previous works [3][4][5][6]13], it is helpful first to integrate out the bosonic fields. This leads to the following representation of the partition function Z, Z =ˆexp −S fermion χ, χ + DχDχ + , (2.21) where the effective fermionic action S fermion [χ, χ + ] takes the form The first and the second terms are the same as those in Eqs. (2.16) and (2.17), while the third one contains both attraction and repulsion due to a special form of the interaction kernel is the solution of the differential equation Putting in Eq. (2.22) U 2 = 0, one arrives at the model considered in Ref. [13]. Putting in additionŨ 0 = 0 one comes to the model studied in Refs. [3][4][5][6]. Both these models contain only the interactions local in imaginary time. On the contrary, in our case, the kernel K (τ − τ |ω q ) contains the additional repulsion term which is non-local in imaginary time. This term is new and very important for the present study. At this point, we would also like to emphasize that the fermion-fermion interactions remain short-ranged in the real space.

A. Mean-field action
Starting with the effective fermionic action S fermion [χ, χ + ], we could, in principle, develop the perturbation expansion in the coupling constants of the interaction terms. However, the framework of the perturbation theory is not suitable for studying the type of problems we consider in this paper, which is typical for strongly correlated systems. The most common alternative is to perform the change of variables from the fermionic degrees of freedom to the bosonic collective degrees of freedom. This transformation is convenient for analytical studies and is absolutely neccessary for the numerical ones.
The standard way to facilitate this change of variables is to decouple the interaction terms with the help of the Hubbard-Stratonovich transformation. In our case, this procedure gives us a model of fermions interacting with auxiliary bosonic fields b q (τ ) and b 1q (τ ) (corresponding toΣ 2 andΣ 1 fermionic terms respectively). Then, the resulting integral over the fermionic fields χ, χ + can be calculated exactly to obtain the final representation of the partition function Z in the form of a functional integral over the fields b q (τ ) and b 1q (τ ) The fields b q (τ ) and b 1q (τ ) in Eqs. (3.1-3.5) satisfy the constraint analogous to the one in Eq. ((2.19)) and obey the periodic boundary conditions The function K −1 (τ − τ |ω q ) is the inverse of the interaction kernel K(τ − τ |ω q ) and equals (see Appendix A) B. Minimum of the action and mean-field equations.
Although Eqs. (3.1-3.2) can serve as a direct calculation procedure, explicit computation of the functional integral (3.1) is still difficult even numerically. This is rather typical problem in study of strongly correlated systems. A standard way to overcome this problem is to start with developing a proper mean-field approximation. In many cases, the mean-field theory allows one to understand properties of new models and figure out what are the possible states, phase transitions between the states, etc. After these first properties are understood, one proceeds with studying fluctuations. Very often they are not so important, at least qualitatively, but it may happen that they lead to significant changes of the mean-field picture. However, starting with the mean-field approximation is the first step that is worth doing.
The mean-field approximation corresponds to the calculation of the functional integral, Eq. (3.1), using the saddle point method. Within this technique one should find the minimum of the action S final [b, b 1 ] and approximate the free energy F as The proportionality of b q=0 (τ ) and b 1,q=0 (τ ) to √ V is typical for condensate functions, and b (τ ) and b 1 (τ ) play the role of order parameters. In the Hamiltonian language one can say that, below the phase transition temperature, a macroscopic number of bosons is located at the state with q =0. At the same time, the fact that b (τ ) and b 1 (τ ) may depend on τ signals about completely new phase transitions and thermodynamic states.
The fields at non-zero q correspond to the fluctuations around the saddle point of the action. We will not consider them in this paper. Neglecting the fluctuations, one can introduce the free energy functional Here,ȟ The equations for the minimum of the free energy functional can obtained by putting to zero its first variation In Eqs. (3.14) and (3.15), the Green functionǦ p (τ, τ ) satisfies the following equatioň We should note that putting U 2 = 0 in Eqs. (3.14) and (3.15), we come to the mean-field equations of Ref. [13].
Equations (3.14) and (3.13) admit both the imaginarytime-independent and the time-dependent solutions. In this paper, we are going to show that in some region of parameters, an imaginary time-dependent solution is energetically more favorable. However, in order to determine the favorable configuration, one not only needs to obtain the different solutions of the mean-field equations but to also calculate the corresponding values of the free energy functional and compare them with each other. As it appears, this is quite a non-trivial task. In Section IV, we study analytically the limiting case of U 2 U 0 . The general case can only be tackled numerically. Thus, in Section V, we formulate a suitable computational scheme to treat the case of general parameters, while the applications of the scheme are discussed in Section VI.
It is worth mentioning that the origin of the interesting physics is the existence of the non-zero imaginary-timedependent order parameter b(τ ). At the same time, the field b 1 (τ ) plays rather a supporting role helping to decrease the free energy for the time-dependent configurations of the field b(τ ). As a consequence, we are going to neglect the field b 1 (τ ) in our analysis of the model to simplify the calculations. On the other hand, the inclusion of this field might be important for realistic description of experiments.
In the caseŨ 0 = 0, U 2 = 0, the exact solutions of Eq. (3.14) are known. However, the free energy of the imaginary-time-independent configuration happens to be the lowest. In the limit U 2 U 0 , we can treat the non-local term in the inverse kernel K −1 (τ − τ |ω 0 ) (see Eq. (3.7)) perturbatively making expansion in U 2 /U 0 . As a result of these procedure, one can obtain analytically the first order corrections to the exact solutions for U 2 = 0 as well as to the corresponding free energies. As we will show, one can identify the region of parameters for which the free energy of the time-independent configuration gets pushed above the corresponding energy of the time-dependent configurations.
If one putsŨ 0 = 0 and U 2 = 0 in Eq. (3.14), it gets transformed into Equations (4.1) and (4.2) has static solutions b(τ ) ≡ ±γ T . The parameter γ T here is determined by the selfconsistency equation Of course, there is also a trivial solution b(τ ) ≡ 0, however we assume that the parameters of the system are such that a non-trivial static solution exists and is more energetically favorable than the trivial one.
The interesting fact is that, besides the static solutions b(τ ) ≡ ±γ T , there is also a whole family of oscillating solutions consisting of the instanton-antiinstanton pairs (see [3][4][5][6]13]), bouncing back and forth between the two static solutions ±γ T . These class of solutions can be written exactly in terms of Jacobi elliptic function sn (x|k): For a solution corresponding to m instantonantiinstanton pairs, the parameters k and γ should satisfy the system of equations The parameter κ p is given by In Eqs. (4.6) and (4.7), K(k) is the complete elliptic integral of the first kind, while Π(n, k) is the complete elliptic integral of the third kind (see, for example, [14,15] to read more about elliptic integrals and elliptic functions). The equation (4.6) tells that the integer number of instanton-antiinstanton pairs should fit onto the interval [0, β]: one of the periods of sn(x|τ ) is 4K(k). Equation (4.7) obtained from the condition of the minimum of the free energy is actually the self-consistency equation in the mean-field theory. In the case of large periods of instanton-antiinstanton pairs for which k → 1, Eq. Since we treat the non-local part of the kernel K −1 (τ − τ |ω 0 ) as a perturbation, it is convenient to separate the corresponding term in the free energy functional, Eq. (3.11). In the limit U 2 /U 0 1, we can neglect U 2 when it appears in combination U 0 + U 2 and write The perturbative expansion can be obtained if we substitute the ansatz into the gap equation Analogously, the corrections to the free energy can be obtained substituting the ansatz (4.11) into Eq. (4.10). The first order correction in the expansion of b (τ ), Eq. So, we write the free energy F in the form C. Comparison of the free energies of the instanton-antiinstanton configurations with the free energy of the static configuration.
Let us consider the static configuration b static (τ ) ≡ γ T and a configuration consisting of m instatonantiinstanton pairs b m,k inst (τ ) = kγ sn(γτ |k) (parameter γ is fixed by the choice of the parameters m and k according to Eqs. (4.6), (4.7), (4.8) and (4.9)). We will denote the corresponding free energies as F static and F m,k inst . Using Eq. (4.13), we can write To simplify the calculations, we consider the limit k → 1. In this limit, instanton-antiinstanton configurations spend almost all the time in the vicinities of the static configurations ±γ 0 . Correspondingly, if one neglects the non-local interaction term in Eq. (3.14), the difference between the action of the instanton-antiinstanton configuration and the action of the static configuration is proportional to the number of the instanton-antiinstanton pairs. Thus, we can write the difference of the free energies without the non-local interaction as Here, we used Eq. (4.6) using the fact that one can take the parameter γ to be equal to γ T for k → 1. Since the limit k → 1 also corresponds to the limit of the zero temperature, we take γ 0 = γ T =0 = γ 0 . The constant ∆S 0 is the action difference for a single instanton-antiinstanton pair: (4.16) The correction to the free energy of the static configuration is evaluated to be The correction to the free energy of the instantonantiinstanton configuration can be calculated in the limit ω 0 /γ 0 1/K(k) using the Fourier expansion for the Ja-cobi elliptic function (see Appendix C). This gives us Combining everything together, we obtain In the limit k → 1, the elliptic integral K (k) diverges, (This corresponds to a large period of the instanton-antiinstanton lattice). So, it is possible to choose k sufficiently close to unity so that Then, if we keep ω 0 sufficiently small, the term quadratic in ω 0 cannot change the sign of the free energy. As a result, the instanton-antiinstanton configuration can really be energetically more favourable. This striking result that already very small coupling constants U 2 U 0 can make the state with instatons more favorable is based on the strong sensitivity of the second term in Eq. (4.10) to whether the solution b (0) is static or consists of the instanton-antiinstanton pairs in the imaginary time τ . In the former case this term can be large, while in the latter case its value can considerably be reduced.
So, we have demonstrated here analytically that there is a region of the parameters of the model, where the instanton crystal exists. Calculations in a more broad region, as well as phase transitions between the states can be studied only numerically and this will be done in the next section.

V. NUMERICAL MINIMIZATION OF THE FREE ENERGY FUNCTIONAL.
At U 2 = 0 andŨ 0 = 0, the solutions of Eq. (3.14) with the different number m of instanton-antiinstanton pairs are topologically distinct. As U 2 is gradually turned on, the topological classes are kept intact. As a result, the solution with m pairs gets deformed yet the period of the configuration W = β/m is preserved.
Alternatively, we can access different solutions of Eq. (3.14) if we minimize the free energy functional (3.11) in the classes of configurations corresponding to different fixed periods W = β/m. In order to turn this recepy into a numerical scheme, we just need to formulate a suitable discretization of the expression (3.11) for the free energy functional, and find a way to enforce the restriction on the period of configurations.
We should also note, that we are going to neglect the field b 1 (τ ) for simplicity. However, the resulting numerical scheme can be easily adapted to take this extra field into account.

A. Transformation of the free energy functional
The free energy functional, Eq. (3.11), consists of two parts: one is the part which is purely quadratic in the fields b(τ ); another part is the fermionic part which originates from the integration of the fermionic degrees of freedom. Discretization of the quadratic part of the free energy functional is a rather straightforward matter. On the contrary, it is impossible to directly discretize the fermionic part of the free energy as it is written in Eq. (3.11). However, it can be recast into the form suitable for the numerical treatment by replacing the functional trace with the expression involving a time-ordered exponential of the energy operatorȟ p (τ ). This can be achieved with the help of the standard trick in the field of Determinant Monte Carlo (see, for example, Ref. [17]). We describe in details the transformation of the fermionic part of the free energy functional in Appendix B.
Using the same regularization as in Appendix B, we write the quantity of interest in our numerical studies as and (see Eq. (B11)) Let us consider configuration b(τ ) with period W such that m periods of the configuration fit into the interval [0, β]. We shall rewrite Eqs. (5.1-5.4) in such a way that only the dependence on values of b(τ ) for τ ∈ [0, W ) explicitly enters the equations. This allows us to fix the period constraint for the optimization procedure.
Notice that the expression for the averaged kernelK 0 (τ − τ |ω 0 ) is identical to Eq. (3.7) for the kernel K 0 (τ −τ |ω 0 ) with the only difference being that β is replaced by W .
Finally, let us defině With the help of this definition, we can rewrite Eq. (5.4) as The goal of discretizing the expression for the free energy is achieved if we replace the function b(τ ) with its values at the discrete set of imaginary time points b i = b(τ i ). For the discretization scheme with N points, we are going to take τ i = (i − 1)∆τ = (i − 1)W/N , where i runs through integer values from 1 up to N . It is a straightforward task to write up the discretized versions of Eqs. (5.5-5.9): is the discrete approximation of the timeordered exponential: Substituting this expressions into Eq. (5.1), we obtain the discretized version of the full free energy functional: To run the optimization procedure, we also need the formulas for the gradient of the free energy The explicit expressions for the matrices appearing in Eqs. (5.13) and (5.18) are (5.20) Here, the parameter κ ip is The scheme we just introduced can be implemented in the programming language of the choice. Since the scheme provides the expressions both for the discretized free energy functional and its gradient, it can be plugged into any first order optimization algorithm.
The analytical solutions of Eq. (3.14) without the nonlocal part of K 0 (τ − τ |ω 0 ) can be used as the initial conditions for the optimization procedure. For fixed W = β/m, this requires to solve the system of equations (4.6-4.9) to determine the parameters k and γ. Since we neglect the non-local part and do not put U 2 = 0, one should replace U 0 by U 0 + U 2 in Eq. (4.7). Then, the initial condition is defined as b (0) i = kγ sn(γτ i |k). Alternatively, one can use instead b where γ T is the solution of the static gap equation (4.3) (but with U 0 + U 2 instead of U 0 ) and k is determined from the condition W = 4K(k)/γ T . In the end, both choices of the initial conditions lead to the same results of the optimization procedure.
C. Variation of the scheme in the limit T → +0.
The numerical scheme we introduced in the previous subsection can be adapted to treat the limiting case of zero temperature, which is equivalent to the limit m → +∞. The only expressions that need to be adjusted are Eqs. (5.12) and (5.17) which define the discretized version of the fermionic part of the free energy functional and its gradient.
Suppose that one calculates the matrixǓ p [b i ] for some specific value of p. The diagonal decomposition of this matrix is given by where we assume that λ 1p is the eigenvalue with the largest absolute value. Then, one can write Eq. (5.12) in the limit m → +∞ as Analogously, Eq. (5.17) transforms into Here, θ(x) is the Heaviside function, while Ǎ i,j denotes the matrix element i, j of the 2 × 2 matrixǍ.

VI. NUMERICAL ANALYSIS OF THE MODEL.
A. Zero temperature.
In order to perform the actual numerical simulations, one needs to specify the fermionic dispersion ε 1,2 (p) (see Eqs. (2.2) and (2.4)). As we explain in Section VII, the model introduced in the present paper originates from the Spin-Fermion model with Overlapping Hotspots (SF-MOHS) studied in references [18][19][20]. As a consequence, we have chosen the fermionic dispersion in the same form as it appears in SFMOHS: where µ is the chemical potential. We also introduce an energy cutoff Λ limiting the width of dispersion: We should note, that in all the computations we neglected the field b 1 (τ ), which is equivalent to setting U 0 = 0.
In Fig. 2, we display the dimensionless difference between the free energies of the instanton-antiinstanton and of the static configurations ∆F/γ 0 = (F instanton − F static ) at zero temperature as the function of the dimensionless period of the instanton lattice γ 0 W and of the dimensionless parameterω 0 /γ 0 which corresponds to the modified frequency of the current-like mode (see Eq. (3.8)). The energy scale γ 0 , which we use to make the physical quantities dimensionless, is the solution of the static gap equation (4.3) at zero temperature and in the abscence of the nonlocal repulsion term. In addition to that, the same equation In each of the subplots of Fig. (2), one can clearly identify the regions where the free energy of the instanton configurations becomes less than the free energy of the static configuration. As a consequence, in these regions the instanton crystal phase should be the one which is thermodynamically stable.
In the instanton crystal phase, the actual period of the lattice is determined by the minimum of the free energy at fixed value ofω 0 /γ 0 . In Fig. (2), the blue curves show the positions of the minima of ∆F/γ 0 as the functions ofω 0 /γ 0 . These minima were extracted by interpolation from the same data used to plot the surfaces. As the value of the parameterω 0 /γ 0 grows, we observe the transition from the instanton crystal phase to the phase with imaginary-time-independent order parameter. In Fig. (2), this transition is marked by the blue points at the end of the blue curves.
Just below the transition, the period and the amplitude of the instanton lattice has finite values, as a result, the transition should be accompanied by an abrupt change in the order parameter. Thus, the transition must be of the first order. In order to prove this point, we display in Fig. (3) the slope d(∆F/γ 0 )/d(ω 0 /γ 0 ) = d∆F/dω 0 of the dimensionless free energy ∆F/γ 0 as the function ofω 0 /γ 0 in the instanton crystal phase at the transition point for the same values of the ratio U 2 /U 0 as in Fig. (3). In each of the four cases, the slope has the finite value in the instanton crystal phase, while it is zero in the static phase, from which one can conclude that in each of the four cases the transition is accompanied by the jump in the first order derivative of the free energy dF/dω 0 = d∆F/dω 0 .

B. Finite temperatures
Besides the quantum phase transitions at zero temperature, it is also interesting to study how the model may enter the instanton crystal phase as the temperature is varied. Thus, we also calculated the dependence of the free energy of the instanton crystal configuration as the function of the inverse temperature.
The results of the calculations are presented in Fig. 4. There, we plot the dimensionless free energy of the instanton crystal configuration and the dimensionless free energy of the imaginary-time-independent configuration as the functions of the dimensionless inverse temperature γ 0 β. The free energies are determined with respect to the free energy of the normal metal configuration (b(τ ) ≡ 0). For the calculations, we used the same parameters of the fermionic dispersion and the same set of ratios U 2 /U 0 as we did for the calculations at zero temperature. In all the cases, the value of the parameterω 0 /γ 0 was fixed: ω 0 /γ 0 = 0.08. We should also emphasize that γ 0 and (U 0 + U 2 )/γ 0 were determined by the static gap equation (4.3) at zero temperature, so that their values stayed constant as the temperature was varied.
For the cases U 2 /U 0 = 0.5, U 2 /U 0 = 1.0 and U 2 /U 0 = 2.0 which correspond to the plots Fig. 4 (a,b,c), we observe that as the temperature is lowered (equivalently, as the inverse temperature grows), the system first undergoes a second order transition to the static phase at temperature T static . As the temperature is lowered even further, there is a transition into the instanton crystal phase at temperature T inst , which corresponds to the intersection of the two free energy curves on the plots (marked by black dots). Since the slopes of the free energy curves at the intersection point are different, this transition is accompanied by a jump in the first derivative of the free energy, henceforth it is of the first order. On the contrary, for the case U 2 /U 0 = 4.0 which correspond to the plot Fig. 4 (d), we observe that the system undergoes transition to the instanton crystal phase whithout ever entering the static phase. Overall, the picture observed in Fig. 4 suggests that as the value of U 2 /U 0 grows larger, T inst moves closer to T static until the point where they coincide. For larger values of U 2 /U 0 , T inst > T static and the transition to the instanton crystal phase happens without the intermediate static phase. The instanton crystal phase has a rather peculiar feature. As the temperature is lowered, more and more instanton-antiinstanton pairs can fit onto the interval [0, β]. As a result, there is a series of the first order transitions characterized by the change in the number of the periods of the instanton lattice m by 1. In Fig. 4, these transitions are marked by red dots.
Finally, we would like to discuss the order of the transition from the normal metal to the instanton crystal phase in the case where is no intermediate static phase involved. This transition is of the second order, which can be understood from the following argument. Let us consider the configurations with a single instanton-antiinstanton pair. Let us also assume for a moment that there is no nonlocal repulsion, so that the instaton-antiinstanton configuration is described by equations (4.5-4.9). The instanton-antiinstanton configuration has a minimal period W 0 which is finite. This period corresponds to the solution of the system (4.6-4.9) with k = 0: W 0 = 4K(0)/γ k=0 . The instanton-antiinstanton configuration has the form b(τ ) = γk sn (γτ |k), thus, as the period gets close to W 0 , k gets close to zero and the amplitude of the configuration vanishes in the limit W → W 0 + 0. Now, if we take into account the non-local repulsion term, its effects should vanish together with the amplitude of the instanton-antiinstanton configuration. From this we can conclude that the minimal period of the configuration would stay the same, and the amplitude of the instanton-antiinstanton configuration would vanish in the limit W → W 0 + 0 as before. The transition to the instanton crystal phase happens when 1/T = W 0 . As the order parameter vanishes at the transition, we expect it to be of the second order. A. Previous models.
The idea of investigating the present model, Eqs. (2.1-2.8), originates from the previous studies of superconducting cuprates using the so-called spin-fermion (SF) model. This phenomenological model has been proposed in order to enable analytical study of low energy physics of cuprates [21][22][23][24]. The philosophy underlying this approach is based on integrating out the high energy degrees of freedom (of the order of the bandwidth) and writing an effective model containing only the low energy excitations. Of course, after such an integration one obtains a very complicated effective Lagrangian that can hardly be treated analytically. In this situation, the only thing that can be done is to simplify the resulting effective model by reducing it to a form containing a small number of different types of the excitations. It is important to have a sufficienly simple form of these excitations and of their interactions.
Originally, Spin Fermion model [21] was introduced as an effective model containing the fermions in the vicinity of the Fermi surface interacting with bosonic antiferromagnetic waves propagating with vector Q close to the antiferromagnetic vector Q AF . The latter are assumed to be the remnants of the parent insulating AF state. A weak interaction between the fermions and the antiferromagnetic waves is most efficient at 8 points of the Fermi surface that can be connected by the vector Q AF (hot spots). The resulting interaction is strongly peaked at the wavevector Q 0 = (π, −π) corresponding to the antiferromagnetic order with vector Q AF and is described by a propagator In Eq. (7.1), v s is the spin velocity and ξ is the correlation length which is supposed to diverge at the antiferromagnetic transtion. It is important to note that the fermions and bosonic spin waves actually have the same origin. The spin waves in the effective model are some complicated collective spin excitations of the bare interacting holes constituting the original microscopic model (we consider the hole-doped cuprates). At the same time, many details of the microscopic model are not so important for the investigation of universal phenomena such as phase transitions, symmetry of the phases, etc.
Having integrated out the high energy fermions, one loses the detailed information about the structure of the lattice. At the same time, one can use Fermi-liquid-like arguments to conclude that one can still use the basic shape of the Fermi surface (see Fig. 5). So, identifying the vector Q 0 with the antiferromagnetic vector Q AF , one can play with the model with 8 hotspots and obtain many interesting resuls at low energies. On the other hand, one can not really exclude the possibility of additional lowlying collective excitations since the Spin Fermion model was introduced on phenomenological basis in the first place. The use of the propagator for the fermion-fermion interaction, Eq. (7.1), is motivated by the proximity to AF quantum critical point (QCP) where ξ AF → ∞. Then, one can consider only small δp ∼ 1/ξ AF vicinities of the 'hot spots' as strongly affected by the interaction. However, at temperatures relevant for e.g. the pseudogap state [25][26][27] this argument does not have to hold because the experimentally reported correlation lengths [28,29] are indeed rather small. Moreover, ARPES experiments [27] show that the effects of the pseudogap extend well beyond the 'hot spots' to the Brillouin zone edges (π, 0), (0, π) without being significantly weakened.
A different version of the SF model (spin-fermion model with overlapping hot spots (SFMOHS)) has been introduced in Refs. [18][19][20]. As ξ AF becomes smaller, the 'hot spots' expand and can eventually overlap and merge forming two 'hot regions' (see Fig. 5). For the latter to occur the fermionic dispersion in the antinodal region should be shallow, which is supported by the experimental data [30,31]. This (SFMOHS) model differs from previously used spin-fermion models with 8 hot spots [21][22][23][24]32] by the assumption that the hot spots on the Fermi surface are not isolated, may overlap and form antinodal 'hot regions'. This can happen when the fermion energies are not far away from the van Hove singularities in the spectrum of the cuprates, which corresponds to the results of ARPES study [30,31,33,34].
The two hot regions 1 and 2 are centered at the middles of the edges of the Brillouin zone and can be connected by the vector Q 0 . Then, one comes to a description in terms of the fermions located in two bands with the interaction between the bands. Due to proximity to the van Hove singularity one can write the spectra of the fermions in the bands near points (π, 0) and (0, π) as 2) The momenta p are counted now from the points (π, 0) and (0, π). In Eq. (7.2), µ is the chemical potential and α and β are constants.
This discussion makes clear the origin of the bare part of the HamiltonianĤ 0 , Eqs. (2.2-2.4) and of the action S 0 [χ, χ + ], Eq. (2.16). For investigating interaction of the fermions via the bosonic spin mode D 0 (ω, q), Eq. (7.1), the field theoretical formulation is more convenient than the Hamiltonian one. As we are intrested in studying possible states different from the aniferromagnets, it is convenient to single out the pairs with assumed strong correlations. For example, one could consider particle-particle pairs when studying superconductivity, or particle-hole pairs for studying charge density waves (CDW) and loop currents [13,[18][19][20][21].
Here, as in Refs. [13,[18][19][20][21], we are interested in loopcurrent channel and charge density channel, both in the vicinity of the vector Q 0 . The channels with the nontrivial spin structure of the pairs correspond to weaker interactions and, thus, are not important.
So, singling out the most interesting pairs, we write the interaction term in the action in the form S int χ, χ + → S (current) χ, χ + + S (density) χ, χ + .
stands for attration of the fermionic loop currents, while stands for the repulsion of the fermion densities oscillating at wave-vector Q 0 in space, λ is the coupling constant of the interaction between the spin mode the spins of the fermions, and X = (τ, r) are vectors containing as components the imaginary time and space coordinates. The fermionic fields χ, χ + have already been intoduced in Sec. II. The signs of the interactions are unambigously determined by the SF interaction. The precise momentum dependence of the propagator D 0 (ω, q) is not important for us. Therefore, to simplify the model, we replace the propagator by a constant in frequency-momentum space or, equivalently, by δ-function in imaginary time and real space. Thus, we arrive at Eq. (2.17) for the interaction part of the action. The inclusion of the long-range part of the Coulomb interaction renormalizes the coupling constant in S (density) [ξ, ξ + ]. As a result, we keep the couplingsŨ 0 and U 0 in Eq. (2.17) as independent constants.
The model only with the interaction S (current) [χ, χ + ] was studied in Ref. [20], and proposed earlier DDW state [35] was obtained in the mean-field approximation. This state corresponds to the static loop current modulated with the vector Q AF and flowing around the CO 2 elementary cells. In this case, the instantonantiinstanton solutions for the order parameter could be obtained, however, the free energies of these configurations were higher than that of DDW state.
In Ref. [13], it was argued that a stable instantonantiinstanton crystal state can be obtained if one adds an extra interaction term S (density) [χ, χ + ] to the model. However, this result was based on an analytical computation invlolving perturbation expansion up to the second order in the abscence of a "small parameter". In order to get a conclusive proof of the stablity of the instantonantiinstanton crystal state, we designed the numerical scheme described in Section V. Unfortunately, although the inclusion of S (density) [χ, χ + ] term lowered the free energy of the instanton-antiinstanton configurations, numerical analysis revealed that DDW state was still energetically more favourable.
On the other hand, in sections IV and VI, we have seen that introducing interaction of the fermions with an additional bosonic mode can drastically change the situation and make the state with the imaginary-timedependent order parameter thermodynamically stable.
B. Interaction of fermions with the bosonic current-like mode.
Integrating out the high energy modes in microscopic models leaves a lot of possibilities for low energy modes obtained after this procedure. This is especially true for cuprates which have a very complicated structure. The only possibility to model the emerging low energy modes is to introduce them phenomenologically.
According to this philosophy, the integration of the high energy degrees of freedom in Spin Fermion model is assumed to lead to the appearance of the bosonic spin modes which are identified with the surviving lowenergy antiferromagnetic fluctuations near the antiferromagnetic quantum critical point. In particular, large values of the propagator D 0 (ω, q) at the antiferromagnetic vector Q AF follow from this fact. These bosonic modes couple to the spins of the fermions corresponding to the magnetic origin of the interaction. Given the phenomenological character of the Spin Fermion model, it is plausible to assume the possibility of the appearance of additional low-energy bosonic modes which couple instead to the magnetic moments induced by the fermion loop currents modulated in space at the wave-vector Q AF .
The bosonic modes described by the Hamiltonian (2.6) are precisely of this type. There, we have a system of oscillators withP q and coordinatesQ q labeled by different wave-vectors q. It is assumed that these wave-vectors are counted from Q AF , and their lengths are not large. The interaction of these modes with the fermionic currents is included in a gauge-invariant manner by coupling the former to the vector potential A q created by fermions (see Eq. (2.8)).
In the Lagrangian formulation, one writes the action S B , Eq. (2.18), with the coordinate field a q (τ ), while the interaction of the bosonic and fermionic currents is given by the term S FB [χ, χ + , a] , Eq. (2.20).
It is worth emphasizing that the Hamiltonian of the model and the corresponding action describe a system of particles with different p and q and, in this sense, do not differ from the standard many body models. The interesting effects show up when one studies the condenstate of particle-hole pairs with dominating contribution at q = 0 (see Eq. (3.10)).

VIII. DISCUSSION AND OUTLOOK.
We have proposed a new thermodynamic model of interacting fermions and bosonic current-like modes. It does not contain any special features like long-range or infinite-range interaction. Using the methods of the field theory we introduced collective boson degrees of freedom and integrated out the fermionic ones, reducing the model to a system of interacting bosons. The model cannot be solved exactly and, as usual, one starts with developing mean-field approximation. In our formulation, the mean-field equations are just equations for the minimum of the bosonic action. This fact considerably simplifies both analytical and numerical study.
We have demonstrated that the system can be (in addition to normal metal phase) either in the stationary phase with a conventional imaginary-time-independent order parameter or in the instanton crystal phase. The numerical investigation at zero temperature, performed in Section VI, reveals the existence of a quantum phase transition between these phases. It is also shown that the derivative of the free energy ∂F/∂ω 0 experiences a jump at the transition indicating that it is a first order transition.
In addition to that, in Section VI, we also performed the numerical investigation of the temperature dependence of the free energy of the model. The results of our calculations indicate that, as the temperature is lowered, the transition from normal metal phase to the instanton crystal phase can happen either via an intermediate stationary phase or directly. In the former case, the transition from the normal metal phase to the stationary phase is of the second order, while the subsequent transition into the instanton crystal phase is of the first order. In the latter case, the direct transition into the instanton crystal phase is of the second order. As the temperature is lowered further, there is a series of first order transi-tions corresponding to the change of the number of the periods of the instanton lattice.
As the results have been obtained using the mean-field scheme, it is important to understand how fluctuations near the saddle-point solution affect the results. For this purpose, one should make expansion of the effective action up to the second order and check the eigenvalues of the corresponding quadratic form. The stability of the long-range order is endangered by the fluctuations associated with the gapless zero modes, and in our case, there is a zero mode originating from the translational invariance of the instanton lattice. On the other hand, we consider the model with at least two spatial dimensions, and the order parameter in the abscence of instantons corresponds to discrete Z 2 symmetry breaking. Also, the imaginary time acts as an extra dimension, which helps to reduce the effect of the fluctuations. Overall, the role of the fluctuations is an open question which we plan investigating on in the future.
It is hard to speak about the possible experimental observation of the instanton crystal phase at this stage. In this paper, we considered only the equilibrium properties of the system. As a result, we can only suggest to look for the discontinuity in the derivative of the free energy. The significant amount of information about the system, on the other hand, can be obtained by measuring its response to various probes. As a consequence, another important direction of our future studies is the analysis of the real-time correlation functions. This task is challenging and deserves a special treatment. Therefore, we decided not to touch the subject in this paper.

ACKNOWLEDGMENTS
Financial support of Deutsche Forschungsgemeinschaft (Projekt EF 11/10-1) and of the Ministry of Science and Higher Education of the Russian Federation in the framework of Increase Competitiveness Program of NUST "MISiS"(Nr. K2-2017-085") is greatly appreciated.
Appendix A: Inversion of the integral operator describing the effective fermion interaction The interaction kernel K(τ − τ |ω q ) was defined in Eq. (2.23) as where the function is the greens function for a certain differential operator (see Eqs. (2.24) and (2.25)) As it turns out, that is all we need to construct the inverse of the integral operator described by kernel K(τ − τ ). Let us denote the kernel for the inverse operator as Let us seek K −1 (τ − τ ) in the form where A is a yet unknown parameter and g(τ − τ ) is an unknown function. In order to determine them, let us plug the ansatz (A5) into Eq. (A4).
It seems reasonable to put A = (U 0 + U 2 ) −1 . Thus, we obtain the following integral equation for the unknown function g(τ − τ ): or, equivalently, Let us apply the differential operator from Eq. (A3) to the both sides of Eq. (A8). This way, we get We can further rewrite it as This equation has the same functional form as the equation (A3) which K 0 (τ − τ |ω q ) satisfies. As a result, we can write the solution right away: Appendix B: Transformation of the fermionic part of the free energy functional.
Let us focus on the part of the free energy functional (3.11) which originates from integrating out the fermionic degrees of freedom: Here, we neglect the field b 1 (τ ), so thať Here, Tr τ,s stands for combined trace in the subspace of anti-periodic functions and in the subspace of the bands 1 and 2, while det τ,s stands for the combined determinant in the same subspaces. In the following, it is also convenient to regularize F ferm by subtracting the constant term corresponding to the normal metal configuration b(τ ) ≡ 0. This way, we obtain (B4) It is hard to work directly with the functional determinants. However, we can re-express the ratio of two functional determinants in terms of the time-ordered exponentials of the corresponding operators (It is the standard trick in the field of Determinant Monte-Carlo. For proof of this relation, see, for example, [17]): where T is the time ordering operator and we used the explicit expressions forȟ p (τ ) andȟ 0p (τ ). We should note that this particular form withǏ + T exp[· · · ] is attributed to the fact that we considered the functional determinants of the operators with anti-periodic boundary conditions. The expression in Eq. (B5) can be further simplified. First, for 2 × 2 matrices, one can show by direct substitution that det Ǐ +Ǎ = 1 + detǍ + TrǍ (B6) Secondly, sinceǏ commutes with any 2 × 2 matrices, we can write in Eq. (B5) and Finally, we should note, that for the time-ordered operator in Eq. (B7), which follows from Liouville's theorem and from the fact that Pauli matrices are traceless. Substituting Eqs. (B6), (B7), (B8) and (B9) alltogether into Eq. (B5), we obtain . (B10) Now, we can rewrite Eq. (B4) as Instead of the normal metal configuration, we could have used a static configuration b(τ ) ≡ γ to regularize the fermionic part of the free energy functional. In this case, one can write Appendix C: Correction to the free energy of the instanton-antiinstanton configurations due to the non-local repulsion.
Let us consider the integral appearing in Eq. (4.18), which describes the correction to the free energy of the instanton-antiinstanton configuration due to repulsive interaction: Here, b m,k inst (τ ) = γk sn (γτ |k) is the configuration with m instanton-antiinstanton pairs and the parameter γ is fixed by the choice of the parameters m and k according to Eqs. (4.6), (4.7), (4.8) and (4.9).
In order to evaluate the integral, it is convenient to use the known Fourier decomposition of the snoidal Jacobi function (see Ref. [14]): where the complementary modulus is k 2 = 1 − k 2 and we used the fact that mT = γ/(4K(k)). In addition to that, we need the Fourier decomposition for the kernel K 0 (τ − τ |ω 0 ): where Ω n = 2πT n are bosonic Matsubara frequencies.
Subsituting Eqs. (C2) and (C3) into Eq. (C1), one ob-tains ∆F repul inst = In the limit k → 1 − 0, we can replace γ by the corresponding value γ 0 = γ T =0 for the static configuration. Also, in this limit K(k ) = K(0) = π/2. In Eq. (C4), the first factor is naturally cut off at |n| ∼ K(k)ω 0 /γ 0 . At the same time, the nonlinearity of hyperbolic sine kicks in for |n| ∼ K(k). If we assume that ω 0 /γ 0 1, the convergence of the series is then determined by the first factor, thus one can safely replace the hyperbolic sine by its argument. As a result, Eq. (C4) transforms into In this Appendix, we would like to mention several details which are important for the speed and stability of the implementation.
The most crucial part is the calculation of the timeordered exponentialǓ p [b i ]. In the case of large period of the configuration W , it is easy to overflow the exponent of the floating-point numbers used to store the matrix elements. This problem may be overcome ifǓ p [b i ] is calculated in extended-precision arithmetics. A complementary solution to this problem is to use the different regularization of the fermionic part of the free energy functional. In Eq. (5.1), we subtracted the constant corresponding to the fermionic part of the free energy of the normal metal configuration. Instead, we could subtract F ferm [γ T ] where γ T is the solution of the static gap equation (4.3) with U 0 replaced by U 0 + U 2 . (The idea here is not to put U 2 = 0 but to neglect the non-local quadratic part of the free energy functional.) In this case, equations (5.12) and (5.17) can be rearranged in such a manner, that, instead ofǓ p [b i ], one needs to calculate which happens to be much more numerically stable. The calculation of non-local part of free energy and its gradient, Eqs. (5.11) and (5.16) requires the computation of matrix-vector products where the matrixK i−j is In fourier space, matrix-vector product of this type reduces to the element-wise multiplication of vectors. As a result, these products can be efficiently computed via the following sequence of steps: calculate fast fourier transform of b j ; multiply it element-wise by the precomputed fourier transform ofK i,j ; make another fast fourier transform.
For the numerical calculations in this paper, we implemented the numerical scheme using the programming language Julia [36]. This language combines the fast prototyping of Python, Matlab and Mathematica with the speed of Fortran, C and C++. For optimization, we used L-BFGS algorithm [37,38] implemented in Optim.jl library [39].

Supplementary Material
for the paper "Phase transition into Instanton Crystal." All equation numbers, figure numbers and reference numbers without prefix "S" refer to the respective numbers in the main text.