Sparse modeling approach to obtaining the shear viscosity from smeared correlation functions

One of the most important results obtained by RHIC (Relativistic Heavy Ion Collider) experiment is an elliptic flow of QCD particles Adams:2005dq ; Adcox:2004mh . The elliptic flow data using a Boltzmann-type equation for gluon scattering indicates a large cross section above the phase transition temperature ($T_{c}$) Molnar:2001ux . Furthermore, this elliptic flow is well explained by the hydrodynamics of the quark-gluon-plasma (QGP) state Teaney:2000cw ; Kolb:2000fha ; Huovinen:2001cy ; Teaney:2001av ; Hirano:2001eu ; Hirano:2002ds ; Teaney:2003kp ; Teaney:2009qa . Thus, the QGP state is not described by a free gas system of perturbed gluons. It might be natural since around $T_{c}$ the interaction force between the quarks and gluons becomes strong. On the other hand, the numerical simulation of the relativistic liquid system gives that the upper limit of the ratio between the shear viscosity ($\eta$) to the thermal entropy ($s$) is very small, $\eta/s\leq 0.4$ Teaney:2009qa . It turns out that the QGP has the perfect liquid property. The result gives us a non-sticky image of QGP states even though the QCD particles are strongly interacting. To make these curious properties of QGP clear, a considerable number of studies have been conducted on the determination of the shear viscosity from the first principle calculation in the pure SU($3$) gauge theory Aarts:2002cc ; Nakamura:2004sy ; Meyer:2007ic ; Moore:2008ws ; Meyer:2008dq ; Meyer:2009jp ; Meyer:2011gj ; Mages:2015rea ; Astrakhantsev:2017nrs ; Astrakhantsev:2018oue ; Pasztor:2018yae .

The shear viscosity is given by the first-differential coefficient of the spectral function at zero frequency ($\eta(T)=\pi d\rho(\omega)/d\omega$ at $\omega=0$). The spectral function is defined by the Euclidean correlation function of the renormalized spacial energy-momentum tensor (EMT) in the static state as follows:

$\displaystyle C(\tau)$

$\displaystyle=$

$\displaystyle\frac{1}{T^{5}}\int d\vec{x}\langle T^{R}_{12}(0,\vec{0})T^{R}_{12}(\tau,\vec{x})\rangle=\int^{+\infty}_{-\infty}d\omega K(\tau,\omega)\rho(\omega).$

(1)

Here, $K(\tau,\omega)$ denotes the kernel of an integral transform. In the lattice simulations, $C(\tau)$ is measured by generated configurations using the first equality sign, and then $\rho(\omega)$ is estimated through the second one.

During these calculations, there are at least three difficulties to obtain $\rho(\omega)$: (i) How to define the renormalized EMT on the lattice (ii) how to improve a bad signal-to-noise ratio of the correlation function of EMT (iii) how to estimate $\rho(\omega)$ from the limited number of the data $C(\tau)$ on the lattice.

As for the first and second problems, a promising method, called the gradient flow method, has been provided in Ref. Suzuki:2013gza ; Asakawa:2013laa . By using the UV finiteness of the flowed operators in the gradient flow Luscher:2010iy ; Luscher:2011bx , we can obtain the renormalized EMT without hard numerical calculation of $Z$-factor for quenched QCD. Furthermore, thank the smearing property of the gradient flow, the statistical uncertainty is suppressed in this method. On the other hand, it is notable that we have to analyze the data only in the fiducial window of the flow-time to obtain the thermodynamic quantities using the gradient flow method in Ref. Asakawa:2013laa . The range is given by $2a\ll\sqrt{8t}\ll T^{-1}=N_{\tau}a$, where both a strong lattice discretization and over-smearing corrections are suppressed.

The third problem is well known as an ill-posed inverse problem. Especially in the finite temperature system, the number of sites in the temporal direction is very limited, then it makes harder to solve the inverse problem.

If we use the gradient flow method to obtain the correlation function, taking the data only in the fiducial window makes the third problem worse, since the available number of the data is further limited. To see the origin of the demerit, we show the smeared correlation function $C(t,\tau)$ in Fig. 1. We find that the slope of $C(t,\tau)$ in the short (and long because of the periodicity) $\tau$-regime shows the discrepancy from the non-smeared data presented in Ref. Meyer:2007ic . The discrepancy becomes large if the distance of two operators ($\tau$) gets shorter than the smeared length ($\sqrt{8t}$) by the gradient flow since the smeared regime of one operator reaches at the position of the other operator in the correlation function.

Now, a question arises: whether we can straightforwardly apply ordinary estimation methods of $\rho(\omega)$ to such a smeared correlation functions. Until now, the fitting with some ansatz (e.g., the Breit-Wigner ansatz) and the Bayesian methods (e.g., the maximal entropy method Nakahara:1999vy ; Asakawa:2000tr ) have been utilized for the estimation of the spectral function in the case of the non-smeared correlation function, $C(\tau)$. For the smeared EMT correlation function obtained by the gradient flow method in $N_{f}=2+1$ QCD, the fitting analyses with the Breit-Wigner ansatz and the hard thermal loop ansatz for $\rho(\omega)$ have been applied Taniguchi:2019eid . The authors have tried to fit the data of $C(t,\tau)$ only in the limited $\tau$-regime to eliminate the over-smeared data. However, the analysis for the eliminated $C(t,\tau)$ using the same functional form of $\rho(\omega)$ is questionable since originally the non-smeared $C(\tau)$ and $\rho(\omega)$ are related via the integration equation Eq.(1) and $C(\tau)$ at each $\tau$ gets the corrections from all $\omega$-regime. Thus, the lack of $C(t,\tau)$ at several $\tau$ may change the functional form of the spectral function. As the other directions, several methods to reconstruct the spectral function with smearing have been also proposed in Refs. Burnier:2013nla ; Hansen:2017mnd ; Tripolt:2018xeo ; Hansen:2019idp ; Kades:2019wtd ; Bailas:2020qmv .

In this work, we propose the sparse modeling method to estimate the spectral function at finite flow-time, which is based on Bayesian statistics. Here, we need not introduce any functional form of the spectral function. The sparse modeling method can be applied to both non-smeared and smeared correlation functions. It has been utilized in board topics, e.g. MRI MRI and taking a picture of a large blackhole Blackhole (see also a recent review Ref. review-sparse ). As for the estimation of spectral function, it has been succeeded for several analyses in the many-body system Shinaoka2017a ; Shinaoka2017b . Utilizing the intermediate-representation basis (IR basis) Shinaoka2017a ; Nagai2018 ; Chikano2018a ; Chikano2018b ; Li2019 and a few reasonable constraints, the sparse modeling makes it possible to estimate the real-frequency representation of $\rho(\omega)$ from a sparse imaginary-time correlation function.

We formulate the sparse modeling analysis for $C(t,\tau)$ only in the fiducial window. Here, we carry out the singular value decomposition (SVD) of the kernel $K(\tau,\omega)$ in the IR basis, which is independent of the Monte Carlo data. The singular values of the kernel decay exponentially or even faster, then we drop the degree of freedom for the correlation function and the spectral function in the IR basis associated with the sufficiently small singular values. Then, we find the spectral function to minimize the square error given by the difference between the input correlation function and the output one reconstructed by $\rho(\omega)$. During this process, the $L_{1}$ regularization term is added to be consistent with the truncation of the degree of freedom in the IR basis. Technically, we utilize the ADMM algorithm to solve the optimization problem with such a regularization term, which is recently proposed by S. Boyd et al. ADMM1 ; ADMM2 . In this work, we demonstrate that the sparse modeling analyses using very limited data for the quenched QCD at finite temperature.

This paper is organized as follows: In §. 2, we give a review of the sparse modeling method. In §. 3.1, we explain how to calculate the correlation function of the renormalized EMT for the quenched QCD in the gradient flow method. In §. 3.2, we give a standard description to obtain the shear viscosity, where we assume that the correlation function in the whole $\tau$-regime is available for the sparse modeling analysis. Then, we modify the standard description to analyze the smeared correlation function in which the over-smeared data are eliminated in §. 3.3. Section 4 presents the simulation setup of this work. We show the result of the sparse modeling analysis using the central value of the smeared correlation function in §. 5. We find that it works well since the reconstructed correlation function from the obtained spectral function is almost consistent with the input data. §. 6 contains the error estimations. We investigate several systematic errors; $\omega_{cut}$, $\tau$-regime, and the flow-time dependence in our analysis, and also estimate the statistical uncertainty.. Although the statistical error of our result is sizable, we show that the bootstrap analysis is promising for estimating the statistical error of this analysis, since an expected relationship between the correlation function and the spectral function in the IR basis is satisfied. In this work, we have only $2,000$ configurations for the measurement of the correlation function, and it is very few in comparison with $800,000$ and $6$ million configurations in Refs. Nakamura:2004sy ; Pasztor:2018yae . Judging from such a poor statistic, we will give up the precise determination of the shear viscosity and focus on the feasibility of the sparse modeling method to obtain the spectral function combing with the gradient flow method. The last section is devoted to the summary.

Here, we give a brief review of the sparse modeling method in the IR basis following Refs. Shinaoka2017a ; Shinaoka2017b (see also a review paper review-sparse ).

The input is the imaginary-time Green’s function $C(\tau)$. Its Fourier transform $C(i\omega_{n})$, where $\omega_{n}$ denotes a Matsubara frequency, is related with the spectral function $\rho(\omega)$ as $\rho(\omega)=\frac{1}{\pi}\mbox{Im}C(\omega+i0)$ by replacing $i\omega_{n}$ with $\omega+i0$. Then, the input $C(\tau)$ is written by the integral form of the spectral function,

$\displaystyle C(\tau)=\int^{+\infty}_{-\infty}d\omega K(\tau,\omega)\rho(\omega),$

(2)

where $0\leq\tau\leq 1/T$. The kernel $K$ in this work is given by

$\displaystyle K(\tau,\omega)=\frac{\cosh\left(\omega(\frac{1}{2T}-\tau)\right)}{\sinh(\frac{\omega}{2T})},$

(3)

but in this section the discussion does not depend on the details of the kernel. Here, we assume that the kernel is described by some exponential function as in many systems.

For simplicity, we denote Eq.(2) using the dimensionless vectors as

$\displaystyle\vec{C}\equiv K\vec{\rho}.$

(4)

The component of the vector $\vec{C}$ is $C_{i}\equiv C(\tau_{i})$, where $\tau_{i}$ labels the temporal site on the lattice with $0\leq\tau_{i}\leq N_{\tau}-1$. The kernel $K$ becomes $N_{\tau}\times N_{\omega}$ matrix, while $\vec{\rho}$ denotes a vector whose component is $\rho_{j}\equiv\rho(\omega_{j})$ with $j=1,\cdots,N_{\omega}$.

Our goal is to find $\vec{\rho}$ to minimize the square error

$\displaystyle\chi^{2}(\vec{\rho})=\frac{1}{2}\|\vec{C}-K\vec{\rho}\|^{2}_{2}.$

(5)

Here, $\|\cdot\|_{2}$ stands for the $L_{2}$ norm defined by $\|\vec{\rho}\|_{2}\equiv(\sum_{j}\rho^{2}_{j})^{1/2}$. Note that the vector $\vec{C}$ is sparse, so that the amount of the information in $\vec{C}$ is much smaller than that in $\vec{\rho}$. The fact leads to the instability of $\vec{\rho}$.

We have only the Monte Carlo data of $C(\tau)$, which is not the same with the true value of $\vec{C}$ and the data have the statistical uncertainty. Here, we call the discrepancy between the Monte Carlo data and the true value of $\vec{C}$ a “noise”.

Equation (4) gives simultaneous linear equations, where the number of equations is $N_{\tau}$ while one of the unknown variables is $N_{\omega}$. We can solve the equation if $N_{\tau}=N_{\omega}$ and the kernel is given by a regular matrix, while the unique solution does not exist and several possible solutions are allowed in the case of $N_{\tau}<N_{\omega}$.

Here, we would like to choose a stable solution against the noise among these possible solutions. For this purpose, we take an efficient basis called an IR basis. By introducing the singular value decomposition (SVD) of the matrix $K$ defined as

$\displaystyle K=USV^{\dagger}.$

(6)

We rewrite Eq.(4) as

$\displaystyle\vec{C}\equiv USV^{\dagger}\vec{\rho},$

(7)

where $S$ is an $N_{\tau}\times N_{\omega}$ diagonal matrix, and $U$ and $V$ are unitary matrices of size $N_{\tau}\times N_{\tau}$ and $N_{\omega}\times N_{\omega}$, respectively. The new vectors in the IR basis

$\displaystyle\vec{\rho}^{\prime}\equiv V^{t}\vec{\rho},~{}~{}~{}~{}~{}~{}~{}~{}\vec{C}^{\prime}\equiv U^{t}\vec{C}.$

(8)

imply the square error Eq.(5) as

$\displaystyle\chi^{2}(\vec{\rho})=\frac{1}{2}\|\vec{C}^{\prime}-S\vec{\rho}^{\prime}\|^{2}_{2}=\frac{1}{2}\sum_{l}(C^{\prime}_{l}-s_{l}\rho^{\prime}_{l})^{2}.$

(9)

Thus, at the minimum point of the square error the $l$-th element almost satisfies

$\displaystyle[\vec{C}^{\prime}]_{l}=s_{l}[\vec{\rho}^{\prime}]_{l}.$

(10)

It turns out that the contribution of $\rho^{\prime}_{l}$ to $\chi^{2}(\vec{\rho}^{\prime})$ is weighted by the corresponding $s_{l}$. The singular values $s_{l}$ with $l=1,2,\cdots$ of this type of kernels decay exponentially or even faster. Then, for example, if double precision numbers are used, $C^{\prime}_{l}$ can not include the information about $\rho^{\prime}_{l}$, where $s_{l}/s_{1}$ is less than $10^{-16}$ on a computer. It does not depend on simulation details since a kernel $K$ is fixed. This fact naturally explains the reason why a tiny noise of $C(\tau)$ leads to a large difference of $\rho(\omega)$. Therefore, a tiny fluctuation in the left-hand side can affect the “best” inference of $\rho(\omega)$, and several “likely” solutions satisfy Eq. (2).

In actual calculations, the component $s_{l}\rho^{\prime}_{l}$ with sufficiently small singular values gives a negligible contribution to $\chi^{2}(\vec{\rho})$. Then, by introducing some threshold $s_{cut}$, we can drop such a component $l>l_{cut}$, where $s_{l}<s_{cut}$. This truncation leads us to obtain a stable solution which is robust against the noise of $C(\tau)$.

Now, the vectors in the IR basis are reduced to the $l_{cut}$-component vectors. Therefore, we would like to find the solution $\rho^{\prime}_{l}$ in Eq. (9), where many components of $\vec{\rho}^{\prime}$ takes zero to be consistent with the truncation. To search for such a solution $\rho^{\prime}_{l}$, we consider the cost function with an $L_{1}$ regularization term,

$\displaystyle F(\vec{\rho}^{\prime})\equiv\frac{1}{2}\|\vec{C}^{\prime}-S\vec{\rho}^{\prime}\|_{2}^{2}+\lambda\|\vec{\rho}^{\prime}\|_{1},$

(11)

where $\lambda$ is a positive constant and $\|\cdot\|_{1}$ denotes the $L_{1}$ norm defined by

$\displaystyle\|\vec{\rho}^{\prime}\|_{1}\equiv\sum_{l}|\rho_{l}^{\prime}|.$

(12)

This $\lambda$ plays the role of the Lagrange multiplier. If the value of $\lambda$ is smaller than the contribution of the noise for $\vec{C}$, then the obtained $\vec{\rho}$ must be consistent with the true spectral function within the uncertainty coming from the noise.

An intuitive picture of the role of the $L_{1}$ regularization is following. The $L_{1}$ regularization ($\|\vec{\rho}^{\prime}\|_{1}=$const.) gives a linear constraint for the components of $\vec{\rho}^{\prime}$, and is graphically drawn by dotted and solid lines in Fig. 2. The minimize problem of the cost function turns to tune $\lambda$ to minimize the sum of two constants: $\|\vec{C}^{\prime}-S\vec{\rho}^{\prime}\|_{2}^{2}$ is a constant and $\|\vec{\rho}^{\prime}\|_{1}$ is the other constant. In Fig. 2, we see that the open-circle point becomes the most favorable. Note that at the open-circle point, $\rho_{2}$ vanishes and the number of the effective components is reduced.

Actually, according to the example analyses in Ref. Shinaoka2017b ; review-sparse , the spectral function becomes featureless for $\lambda>\lambda^{opt}$, while artificial spikes appear for $\lambda<\lambda^{opt}$. In other words, the former case corresponds to the under-fitting, where the $L_{1}$ regularization term is too strong and the number of components $\rho^{\prime}_{l}$ is too reduced. On the other hand, the latter case does to the over-fitting, where the $L_{1}$ term is too weak and the vector $\vec{\rho}^{\prime}$ has redundancy.

The optimization problem with these $L_{1}$ and/or $L_{2}$ regularization is called the LASSO (Least Absolute Shrinkage and Selection Operators) problem LASSO . It allows obtaining the global minimum regardless of initial conditions ADMM1 ; ADMM2 . As an additional constraint, we will add the non-negativity condition of the spectral function,

$\displaystyle\rho(\omega)\geq 0,$

(13)

in this work. Furthermore, we can also use the sum rule,

$\displaystyle\sum_{j}\rho_{j}=1,$

(14)

as an additional condition, while this constraint is not used in this work. The numerical algorithm to solve the optimization problems with these constraints has been developed ADMM1 ; ADMM2 , and it is called the alternating direction method of multipliers (ADMM) algorithm. See Appendix A for the detail of the algorithm.

In this section, we explain how to obtain the shear viscosity for the quenched QCD in the lattice simulation. Firstly, we explain the calculation strategy of the correlation function of EMT using the gradient flow. Secondly, we give a standard formula for the sparse modeling method. Here, we assume that the correlation function in the whole $\tau$-regime is available for the sparse modeling analysis. Then, we modify the standard formula for the smeared correlation function at finite flow-time, where we eliminate the over-smeared data.

We consider the Wilson plaquette gauge action under the periodic boundary condition at $\beta=6.93$ on $N_{s}^{3}\times N_{\tau}=64^{3}\times 16$ lattices. The lattice parameter realizes the system at $T=1.65T_{c}$. We use the relation between $a/r_{0}$, where $r_{0}$ denotes the Sommer scale, and $\beta$ in Ref. Guagnelli:1998ud . Then, $T/T_{c}$ is fixed by the resultant values of $Tr_{0}=(N_{\tau}(a/r_{0}))^{-1}$ using the result at $\beta=6.20$ in Ref. Boyd:1996bx . The gradient flow method in Ref. Asakawa:2013laa gives the thermal entropy, $s/T^{3}=4.98(24)$, after $a\rightarrow 0$ and then $t\rightarrow 0$ extrapolations.

Gauge configurations are generated by the pseudo-heatbath algorithm with the over-relaxation. We call one pseudo-heatbath update sweep plus several over-relaxation sweeps as a “Sweep”. To eliminate the autocorrelation, we take $200$ Sweeps between measurements. The number of gauge configurations for the measurements is $2,000$. Note that it is quite a small number of configurations in comparison with $800,000$ in Ref. Nakamura:2004sy and $6$ million in Ref. Pasztor:2018yae .

Judging from such a poor statistic, we give up a serious estimation of the statistical error and focus on the feasibility of the sparse modeling method to obtain the shear viscosity combing with the gradient flow method. In §. 5 we explain how to estimate the spectral function using the central value of the correlation function. Then, in §. 6, we discuss the systematic and statistical errors in the analysis. We demonstrate the bootstrap analyses to estimate statistical uncertainty.

The flowed gauge field is obtained by solving the ordinary first-order differential equation (Eq. (16)). Numerically, we utilize the third-order Runge-Kutta method in which the error per step ($t\rightarrow t+\varepsilon$) is ${\cal O}(\varepsilon^{5})$. We take $\varepsilon=0.01$, and confirm that the accumulation errors are sufficiently smaller than the statistical errors. The gauge action of the flow is the Wilson plaquette gauge action.

Now, we demonstrate the sparse modeling analysis using the central value of the measured correlation functions. Table 1 shows the technical parameters in the analysis in this section.

The reason why we take these values of $s_{cut}$ and $a\omega_{cut}$ will be explained below and §. 6.1. Fixing $[\lambda_{min},\lambda_{max}]$, $N_{\lambda}$, and $(\mu,\mu^{\prime})$ are correlated with each other. This set is a possible choice to find an optimal value of $\lambda$.

In Step 1’ listed the end of §. 3.3, we measure the correlation function using the gradient flow method in the range of the flow-times $0.50\leq t/a^{2}\leq 2.50$ with the interval $\Delta t/a^{2}=0.10$. Take $t/a^{2}=0.50,1.00,1.50,2.00$, and $2.50$ for example in this section. The data of $C(t,\tau/a)$ appears in Fig. 3.

We see from Fig. 3 that the longer flow-time data have a smaller statistical error. Actually, in these poor statistics, namely $2,000$ configurations, the correlation function of the non-smeared $U_{12}$ operator is quite noisy and its central values at some $\tau/a$ often take negative values because of the large fluctuations. However, all flowed data shown in the left panel in Fig. 3 (except for only one data-point at $\tau/a=8$ at $t/a^{2}=1.00$) indicate correctly positive values. It is a great advantage of the usage of the gradient flow method to obtain the correlation function.

On the other hand, there is a demerit of the usage of the gradient flow for $C(t,\tau/a)$. Thus, the smearing changes the slope in the short $\tau$-regime because of the over-smearing. To see the merit and demerit clearly, we rescale the correlation function by multiplying the flow-time, $(t/a^{2})C(t,\tau/a)$ (right panel in Fig. 3). In $\tau/a\lesssim 3$, the slope of the correlation function strongly depends on the flow-time. It indicates that the kernel term, which is proportional to hyperbolic cosine function, is not available for the whole $\tau$-regime in the flowed correlation function (see also Appendix B). On the other hand, around $\tau/a=N_{\tau}/2$ the curve of the correlation function does not change, while the statistical errors are reduced thanks to the smearing. A similar property of the smearing has been reported in the case of the APE smearing for the Wilson loop (see Fig. $5$ in Ref.Bilgici:2009kh ). To take only the merits of the gradient flow, it is better to take the fiducial $\tau$-regime at finite flow-time, $\sqrt{8t}<\tau$, for the estimation of $\rho(\omega)$. Here, we fix the $\tau$-regime as $3\leq\tau/a\leq 13$ and investigate the flow dependence of the results.

The next step, Step 2’, is the SVD decomposition of the kernel matrix. We introduce the cutoff $\hat{\omega}_{cut}=4.0$ and discretize $a\omega$ in the range of $-4.0\leq a\omega\leq 4.0$ into $N_{\omega}=3001$ data. The singular values of the SVD decomposition appear in the left panel of Fig. 4.

Note that the vertical axis takes a log scale. The figure tells us that a large hierarchy more than $10^{10}$ exists between $6$-th and $7$-th singular values if we take $3\leq\tau/a\leq 13$. Furthermore, we find that the number of $s_{l}$ larger than $10^{-16}$ is the same with the number of independent $\hat{\tau}_{i}$ under the periodic boundary condition. Refer to several condensed matter studies Chikano2018a ; Chikano2018b ; Li2019 , it will saturate to a specific value even though the number of the data points for $C(\tau)$ increases. It indicates that the number of $\tau_{i}$, here at most $7$ in taking $2\leq\tau/a\leq 12$ under the periodic boundary condition, to analyze this type of kernels is very sparse ⁴⁴4 We can also estimate a lattice size in the temporal direction to fully analyze the kernel (see Appendix C).. To minimize the square error in Eq. (9), the contributions of $s_{l}\rho^{\prime}_{l}$ with a sufficiently small $s_{l}$ are negligible. Therefore, we introduce the threshold of the singular value $s_{cut}$ as $s_{cut}=10^{-10}$ in the present analysis.

The most important point is that this property of the singular values is determined only by the kernel matrix Eq. (41) and is independent of the simulation details. Here, we just assume that the correlation function can be described by a hyperbolic cosine function.

The right panel in Fig. 4 depicts the $N_{\omega}$ dependence of the singular values using $3\leq\tau/a\leq 13$. We find that $s_{l}$ is almost independent of $N_{\omega}$ if we include $\sqrt{\Delta\omega^{\prime}}$ in the kernel matrix as Eq. (41).

It is worth to see the correlation function in the IR basis,

$\displaystyle C^{\prime}_{l}=\sum_{\hat{\tau}=3}^{13}U^{t}(l,\hat{\tau})C(\hat{\tau}).$

(44)

Figures 5 depicts $C_{l}$ as a function of $l$. In comparison with the right panel where $4\leq\tau/a\leq 12$, the values of $C_{l}$ is large. We take $3\leq\tau/a\leq 13$ and will discuss the results using $4\leq\tau/a\leq 12$ in §. 6.2. We see that $C_{l}$ with $l\gtrsim 5$ is also sufficiently small, so that it is consistent with the truncation of $s_{l}$ for $l\geq 7$.

Step 3’ is the $\lambda$ optimization. The $\lambda$ dependence of the square error (Eq. (5)) is a monotonically increasing function of $\lambda$ Shinaoka2017b . The optimized $\lambda$ is given by $\lambda^{opt}=1.1\times 10^{-7}$ – $1.9\times 10^{-6}$ in our analysis. To find this, we scan the value of $\lambda$ in the range of $10^{-15}\leq\lambda\leq 10^{2}$ by changing its exponent with $1/N_{\lambda}$ interval. Once we obtain $\lambda^{opt}$, then we check whether it shows roughly an inverse-rescaling with $(\mu,\mu^{\prime})$.

The results of the spectral function are shown in Fig. 6.

First of all, we find the tails of spectral function for all flow-times approach to zero. It tells us that $a\omega_{cut}=4.0$ is a reasonable choice in this analysis. Secondly, we see the integration of $\tilde{\rho}(t,a\omega)$ in terms of $a\omega$,

$\displaystyle{\mathcal{N}}=\int_{-a\omega_{cut}}^{a\omega_{cut}}\tilde{\rho}(t,a\omega)d(a\omega),$

(45)

monotonically decreases as a function of flow-time. In more detail, the spectral functions in large $|\omega|$ regime are highly suppressed in the longer flow-time. The gradient flow gradually reduces the degree of freedom with high-frequency and can be interpreted as a renormalization group flow. We can see that the results of the sparse modeling analysis give a good account of such an intuitive picture.

The left panel of Fig. 6 is an enlarged plot, which focuses on $\omega\approx 0$. The curvatures of $\tilde{\rho}(t,a\omega)$ at $\omega=0$ has an opposite sign between $t/a^{2}=0.50,1.00$ and $t/a^{2}=1.50,2.00,2.50$. It remind us that the gradient flow smears the data in $\tau/a<\sqrt{8t}/a$, then the data at $\tau/a=3$ (and $13$) are over-smeared in $t/a^{2}\leq 1.12$. Thus, the results in $t/a^{2}=0.50,1.00$ may suffer from the over-smearing corrections. To conclude this, we have to investigate the statistical uncertainty of the input $C(t,\tau/a)$. We will discuss this point in §. 6.