Numerical Study of the Thermodynamic Uncertainty Relation for the KPZ-Equation

A general framework for the field-theoretic thermodynamic uncertainty relation was recently proposed and illustrated with the $(1+1)$ dimensional Kardar-Parisi-Zhang equation. In the present paper, the analytical results obtained there in the weak coupling limit are tested via a direct numerical simulation of the KPZ equation with good agreement. The accuracy of the numerical results varies with the respective choice of discretization of the KPZ non-linearity. Whereas the numerical simulations strongly support the analytical predictions, an inherent limitation to the accuracy of the approximation to the total entropy production is found. In an analytical treatment of a generalized discretization of the KPZ non-linearity, the origin of this limitation is explained and shown to be an intrinsic property of the employed discretization scheme.

production in terms of mean and variance of an arbitrary current, provided the system is in a non-equilibrium steady state (NESS), for recent reviews see, e.g., [4,5]. Specifically, the TUR product Q of entropy production ∆s tot and precision 2 of a current, both defined precisely below, obeys Q ≥ 2. In [6], we have proposed a general framework for formulating a field-theoretic thermodynamic uncertainty relation. To demonstrate this framework, we have analytically shown the validity of the TUR for the one-dimensional Kardar-Parisi-Zhang equation (KPZ) [7]. As a central result, we found that the TUR product Q is equal to 5 in the limit of a small coupling constant, i.e., the TUR is not saturated in this case. Since its introduction, the KPZ equation has been studied extensively and evolved into one of the most prominent examples in non-equilibrium statistical physics. An overview of the progress made can be found in, e.g., [8,9,10,11]. Regarding more recent theoretical developments on the aspect of the KPZ probability density functions and their respective universality classes we mention, e.g., [12,13,14]. Another active area of theoretical work deals with different types of correlated noise, see, e.g., [12,15,16,17]. Regarding experimental studies of the KPZ equation via liquid crystal turbulence, we refer to [18,19,20]. Recent numerical treatment of the KPZ equation is shown in, e.g., [21,22]. From a mathematical point of view, in [23] a space-time discretization scheme for the equivalent Burgers equation has been proposed and its convergence, albeit in a weak distributional sense, has been rigorously proven. In the present paper, we perform a numerical study of the KPZ-TUR to confirm the analytical results from [6] via a direct numerical simulation based on finite difference approximation in space [24,25,26]. Due to the poor spatial regularity of the KPZ equation the discretization of the non-linearity (∂ x h) 2 is not straightforward. There exists a variety of different procedures, which lead to significantly differing results regarding expressions like the surface width (see, e.g., [24,25,26]). In [27] a generalized discretization of the KPZ non-linearity has been introduced, which covers most of the above mentioned schemes. This generalization uses a real parameter 0 ≤ γ ≤ 1 to tune the explicit form of the respective discretization. For γ = 1/2, one obtains the so-called 'improved discretization' introduced in [24], which distinguishes itself by remarkable theoretical properties (see, e.g., [24,27] and subsection 3.1). For the equivalent Burgers equation a closely related scheme is used in [23], which results in the so-called Sasamoto-Spohn discretization of the Burgers non-linearity [23,28]. We perform the numerical simulations in section 5 with this improved discretization scheme (γ = 1/2). Moreover, in subsection 6.5, we also use the boundary cases of γ = 0, 1 as these turn out to mark the lower and upper bound, respectively, of the numerical (discrete) TUR product (see subsubsection 6.3.4). While the discretization with γ = 1/2 is the best approximation to the results from [6], we find numerically in subsection 5.3 that there is still a deviation of roughly 10% between the numerical results and [6]. Based on an idea presented in [29], we can explain this systematic deviation by an analytical test of the generalized discretization of the KPZ non-linearity. We show that this deviation is an inherent property of the generalized nonlinear discretization operator (see section 6), which we think is an intriguing finding in itself. The paper is organized as follows. The basic notions for formulating the field theoretic TUR are introduced in section 2. In section 3, we explain the spatial and temporal discretization of the employed numerical scheme. Our way of approximating the TUR constituents and their theoretically expected scaling according to [6] is presented in section 4. In section 5, we show our numerical results; section 6 is devoted to the analytical treatment of the generalized discretization of the KPZ non-linearity. We draw our conclusions in section 7.

Basic Notions and Problem Statement
We start with a brief sketch of the underlying continuum problem and the notions needed to formulate the TUR. Consider the (1 + 1)-dimensional KPZ equation modeling nonlinear surface growth with h = h(x, t) the surface height on a finite interval in space, x ∈ [0, b], b > 0, In (1) we employ periodic boundary conditions, h(0, t) = h(b, t), and vanishing initial condition, h(x, 0) = 0, x ∈ [0, b], i.e., we start from a flat surface. The KPZ equation from (1) is subject to Gaussian space-time white noise with zero mean, η(x, t) = 0 and covariance given by where ∆ 0 measures the noise strength. The parameter ν describes the strength of the diffusive term (surface tension) and λ is the coupling constant of the non-linearity that models surface growth perpendicular to the local surface. One constituent of the TUR is the so-called fluctuating output, or, equivalently, the time-integrated generalized current, given by a linear functional, which reads [6] Here g(x) ∈ L 2 (0, b) describes an arbitrary weight function with non-vanishing mean (i.e., b 0 dx g(x) = 0). The precision of the output functional from (2) is given by with · as the average over the noise history. In the NESS, i.e., for t 1, the precision from (3) becomes independent of the weight function g(x) [6]. The second component of the TUR product is the total entropy production in the NESS. For the KPZ equation from (1) the total entropy production is given by [6] Based on the experience from [1,3] for a Markovian dynamics, the TUR product Q is expected to fulfill In [6] we have shown analytically for the KPZ equation that for small λ i.e., the TUR is obviously fulfilled, however not saturated. In the present paper, we numerically obtain the two TUR constituents in (3) and (4), and hence Q, by direct numerical simulation of (1).

Discretization of the KPZ-Equation
Throughout this paper we will use a direct numerical integration technique to simulate the height h(x, t) of the Kardar-Parisi-Zhang equation. There are various approaches to this regarding spatial and temporal discretization (see, e.g., [23,24,26,27,28,29,30,31,32,33]). In the following we present the details and reasoning of our approach.

Spatial Discretization
We consider a one-dimensional grid with grid-points x l subject to periodic boundary conditions with lattice-spacing δ given by where b is the fixed length of the grid and L is the number of grid-points. At each grid-point we have for a fixed time t the value of the height field h l (t) ≡ h(x l , t) = h(lδ, t), with x l = lδ and l = 0, . . . , L−1. The time evolution of h l (t) is then governed by (1), i.e., where h L (t) = h 0 (t) due to the periodic boundary conditions. Furthermore, L l and N l denote the discretizations of the linear and nonlinear term at the grid-point x l , respectively, and η l (t) ≡ η(x l , t) = η(lδ, t) represents the discretized noise. Regarding the diffusive term L l in (8), we choose the standard discretization, namely the nearest-neighbor discrete Laplacian, see, e.g., [24,25,27,28,32]. The discretization of the nonlinear term N l is more subtle. During the last few decades different discretizations of the nonlinear term have been proposed for numerically integrating the KPZ equation [24,25,27,32]. In the case of one spatial dimension, they all belong to the family of so-called generalized discretizations [27], with γ ∈ R and 0 ≤ γ ≤ 1. In the following, we will highlight the cases γ = 0, γ = 1 and γ = 1/2. For γ = 0, this discretization reads which is simply the arithmetic mean of the forward and backward taken slope, respectively, of the height field at the grid-point x l [27]. The case γ = 1 yields This is the square of the central difference discretization of ∂ x h, which is a commonly used choice for numerically integrating the KPZ equation, see, e.g., [24,32].
This form was applied to the KPZ equation in, e.g., [24,25,26]. It is closely related to the discretized non-linearity proposed in [23,28,31] of the 1d-Burgers equation equivalent to (8). Following [24], we will name (13) the improved discretization (ID), for the following reasons. It has been shown analytically in [24] using (13) that the discrete Fokker-Planck equation corresponding to (8) possesses a steady state probability distribution for all λ > 0, which is equal to the linear (Edwards-Wilkinson, λ = 0) steady state distribution. It was further shown that the stationary solution of the Fokker-Planck equation is reached for a non-vanishing conserved probability current, which indicates a genuine non-equilibrium steady state in the discretized system. This implies that the case γ = 1/2 accurately mimics the NESS-behavior of the continuous case, with the exact form of the total entropy production ∆s tot from [6]. The above mentioned properties of the operator N (1/2) l distinguish the case γ = 1/2 from, e.g., γ = 1, which does not fulfill the fluctuation-dissipation relation in (1 + 1) dimensions that is essential for obtaining the discrete NESS probability distribution equivalent to the continuous case. Furthermore, the choice γ = 1/2 in (10) is the only one that displays the above behavior [27]. We note that for spatially smooth enough functions h, any discretization from (10) (0 ≤ γ ≤ 1) has an approximation error O(δ 2 ). This implies that for sufficiently small δ the differences between their respective outcomes can be made arbitrarily small. However, the solution h(x, t) of (1) is at every time t a very rough function in space (see also section 6). The various discretizations in (10) thus lead to significantly different results, e.g., with respect to the surface width in [24,25,29] and in the present paper with respect to certain integral norms of the KPZ non-linearity being essential for the KPZ-TUR (see section 6).

Temporal Discretization
Regarding the temporal discretization of (8), we choose the stochastic Heun method (see, e.g., [30,34]), as its predictor-corrector nature reflects the Stratonovich discretization used in [6]. To be specific, the predictor step applies the Euler forward scheme to (8), which yields the predictor y l (t + ∆t) according to Here, l = 0, . . . , L − 1 like above and {ξ l (t)} are stochastically independent N (0, 1)-distributed random variables (see, e.g., [30,32,34]). The prefactor in front of ξ l (t) ensures that the noise has the prescribed variance according to (1). The predictor from (14) is then used in the subsequent corrector step as The form in (15) displays the above mentioned Stratonovich time discretization. For the sake of simplicity, we start at t = 0 from a flat profile, in particular h l (0) = 0, l = 0, . . . , L − 1, and we impose periodic boundary conditions, i.e., h L (t) = h 0 (t). We slightly reformulate the expressions in (14) and (15) by introducing a set of effective input parameters { ν, ∆ 0 , λ} given by with δ from (7). Hence, the predictor-corrector Heun method reads where we set The effective spatial step-size ∆x in the simulation is now simply given by which is a common choice, see, e.g. [25,26,30,32]. From the parameter set { ν, ∆ 0 , λ}, which enters the simulation, the physical parameter set {ν, ∆ 0 , λ} can be obtained from (16). Finally, the calculation of the constituents of the TUR requires expectation values, denoted by · · · . Those are approximated by ensemble-averaging over a certain number E of independent realizations.

Regularizations
Since the KPZ equation is strictly speaking a singular SPDE (see, e.g., [23,29,35]), it has to be regularized in some way. From a physical point of view, this can be done by either introducing a smallest length-scale (e.g., in form of a lattice-spacing δ [28]) or, in Fourier-space, by defining an upper cutoff wave number [25]. In the course of the analytical derivation of a KPZ-TUR in [6], we took the second approach and introduced the cutoff wave number 2πΛ/b. This caused the physical entities like output functional, diffusion coefficient and entropy production rate to depend on this cutoff parameter (see eqs. (80), (85) and (110), respectively, in [6]) and to become singular for Λ → ∞. Here, we use the real-space direct numerical simulation, described in section 3, with lattice-spacing δ from (7) to calculate the relevant physical quantities, which will depend on δ and diverge for δ → 0. For comparison purposes, a relation between the cutoff parameter Λ and the lattice-spacing δ = b/L has to be established. To this end, consider a function f (x, t), x ∈ [0, b], the values of which are known at L grid-points x l = l δ, l = 0, . . . , L − 1. Then the discrete Fourier-transform of (∂ x f ) 2 calculated at these x l is exact for all Fourier modes |k| ≤ Λ, with This is the content of the 3/2-rule by Orszag [36]. It is used in pseudo-spectral methods, as the so-called dealiasing procedure (see, e.g., [25]). With this relation between the number of grid-points L and the wave number cutoff Λ, we will now proceed with the numerical approximation of the TUR constituents and their respective scaling forms.

Mean and Variance of the Output Functional
As we showed in [6], the KPZ-TUR is independent of the choice of g(x) in (2) and thus, for the sake of simplicity, we set g(x) ≡ 1 in the following. The output functional thus becomes the instantaneous spatially averaged height. We define i.e., via a composite Simpson's rule, with periodic boundary conditions h (17) and ∆x = 1 from (19). This implies that we approximate L Ψ (t)/b, rather than (21) itself, which simplifies the comparison of the numerically obtained results with the theoretical ones. The expected scaling of Ψ (N ) (t) is derived as follows. From [6] the corresponding (dimensionless) theoretical prediction Ψ s (t s ) is known to lowest non-vanishing order in λ eff as Here, x s ≡ x/b, t s ≡ t/T and h s ≡ h/H are scaled, dimensionless quantities with reference values b, T = b 2 /ν and H = (∆ 0 b/ν) 1/2 , respectively, and represents an effective, dimensionless coupling constant, see [6]. Hence, after rescaling, (23) can also be written as where (16) and (20) were used. The left hand side of (25) is what we approximate with Ψ (N ) (t) from (22), and thus is the expected scaling behavior in the number of grid-points L and time t for t 1. For the variance of Ψ (N ) (t), we know from [6], that for sufficiently large Λ to lowest non-vanishing order in λ eff . Hence, by following the same steps as above, we get for t 1. Using (26) and (28) the scaling form for the precision is given by for t 1.

The Total Entropy Production
The last entity missing for formulating the numerical version of the KPZ-TUR is the total entropy production ∆s tot . It is given as defined in (4) where we note that the integrand on the r.h.s. is nothing but the square of the KPZ nonlinearity (∂ x h) 2 and we thus can approximate the integrand using any of the discretizations from (10). To be specific, by means of the composite Simpson's rule, we get the approximation The prefactor of L 3 /b 3 arises from the fact that (31) uses (17) with ∆x = 1. Lastly, the time integral in (4) is approximated via where ∆t is a discrete time-step and t = N ∆t.
Proceeding similarly like in subsection 4.2, we get from the theoretical predictions in [6] the expected scaling behavior for ∆s tot (N ) and the TUR product with respect to the number of grid-points L and time t as for t 1, and 5 Numerical Results

Employed Parameters and Fit Functions
In this section we present the numerical results obtained from (17) by using three different discretizations according to (10). If not explicitly stated otherwise, we employ for the numerical simulations the ID discretization (γ = 1/2) from (13). We compare the numerics to the expected scaling forms from section 4. The numerics is performed for the following set of input parameters. For all simulations we set ν = ∆ 0 = 1 and take λ from λ = 0.01 to λ = 0.1 on a range of grid-points, which varies from L = 16 to L = 1024. In the range of L = 16 . . . 64 we use a time-step size of ∆t = 10 −4 and an ensemble size of E = 500. For L = 128 . . . 1024 a larger time-step of ∆t = 10 −2 and smaller ensemble size, E = 250, is used. This reduction is due to the strongly increasing run-time of the simulations for larger numbers of grid-points. We test the scaling predictions for the TUR constituents. At first, we check whether (26) and (28) is fulfilled. This is done by fitting the numerical data of Ψ (N ) (t) 2 and var Ψ (N ) (t) according to the fit-function f 1 , with and f 2 , given by respectively, where a L and b L are L-dependent fit-parameters. Subsequently, we compare a L and b L with c 1 (L) and c 2 (L), respectively. The scaling prediction for the precision ( 2 ) (N ) according to (30) is tested by fitting with fit-parameter d L , to the numerically obtained data for the precision and comparing c 3 (L) to d L for the respective values of L. Finally, by fitting the numerical data for ∆s tot (N ) according to with e L as the fit-parameter and subsequently comparing c 4 (L) to e L , the scaling prediction for the total entropy production (33) is checked.    (35), (36), respectively, for the fits as shown in Fig. 1.

Expectation Squared and Variance for the Spatial Mean of the Height Field
denote the absolute values of the respective relative errors in percent.
indeed valid. The approximation becomes more accurate for a growing number of grid-points L as the relative error ∆ 1 shows the clear trend of decreasing for growing L. The slight deviation in this trend observed between L = 64 and L = 256 is due to the fact that we changed ∆t from ∆t = 10 −4 for L = 64 to ∆t = 10 −2 for L = 256 as well as E = 500 for L = 64 to E = 250 for L = 256. However, as the effect is rather small, we did not see the need to adjust the parameters ∆t and E for L = 256 . . . 1024 in order to achieve a higher accuracy. We now turn to the variance of the mean height field according to (28). The predicted power-law behavior of var Ψ (N ) (t) in (28) can be observed in Figs. 1(c) and 1(d). The predicted scaling factor c 2 (L) in (28) is reproduced well by the numerical data and its respective fit-functions (36) with fit-parameter b L . In contrast to the results for Ψ (N ) (t) 2 , no clear trend in the relative error ∆ 2 = |c 2 (L) − b L |/c 2 (L) can be observed, i.e., ∆ 2 does not become smaller with growing L. An improvement in the approximation may be obtained by an increase of the ensemble-size E and a further decrease of the time-step size ∆t. This, however, would lead to a significantly longer run-time of the simulations. As in all cases the relative error is below 10%, the gain from a further improved accuracy following the above mentioned steps may be outweighed by the increasing run-time.
We also performed the same numerical simulation for λ = 0.01 (data not explicitly shown) instead of λ = 0.1 before. For Ψ (N ) (t) 2 this causes the system to take longer to reach its NESS-behavior, namely t ≈ 10 4 in comparison to t ≈ 10 2 for λ = 0.1. This is due to the weaker driving by the non-linearity weighted with λ = 0.01 opposed to λ = 0.1 in Fig. 1 with { ν, ∆ 0 } fixed. Nev-ertheless, the general trend that with an increase of the number of grid-points L a decrease in the relative error ∆ 1 is achieved could still be seen clearly. Regarding the variance of Ψ (N ) (t), we could not determine a significant difference between λ = 0.01 and λ = 0.1.    To be specific, the numerical data converges to the asymptotic behavior according to (30) after t ≈ 10 2 . This is roughly the same time it took for Ψ (N ) (t) 2 in Fig. 1. That these two times coincide is to be expected, since var Ψ (N ) (t) did not show a discernible amount of time to converge to the asymptotic scaling form (see Fig. 1).

Precision and Total Entropy Production
We will now turn to the scaling behavior of the total entropy production ∆s tot (N ) . In Figs. 2(c), 2(d), we show the plots of the numerically obtained data for ∆s tot (N ) and the according fits. We find that the scaling behavior is recovered nicely, albeit with a significantly greater relative deviation as compared to ( 2 ) (N ) (see Tab. 2).
We also calculated the precision for λ = 0.01 (data not explicitly shown), where it could again be seen that the time needed by the system to reach its NESS-behavior is at least one order of magnitude longer for λ = 0.01 than for λ = 0.1. This is due to the same reason as discussed for Ψ (N ) (t) 2 above. It becomes apparent for ∆s tot (N ) that the only effect of the reduction of λ by one order of magnitude from λ = 0.1 to λ = 0.01 is a rescaling of the scaling factors. While for the other entities like Ψ (N ) (t) 2 and var Ψ (N ) (t) some impact by the change of λ can be observed in regard to the scaling factors and the relative errors, the relative errors of the scaling factors of ∆s tot (N ) do not change significantly with λ. In particular, the only influence on the relative error of ∆s tot (N ) is achieved by an increase in the number of grid-points L. It seems, however, that the relative errors do not become smaller than roughly 10% even for large L. In section 6 we will discuss this observation in more detail and present an analytical explanation for this discrepancy.

Thermodynamic Uncertainty Relation
In the previous sections we have derived the scaling forms of all the TUR constituents and tested their scaling predictions numerically. Here, we will combine these results for the numerical thermodynamic uncertainty product Q (N ) = ∆s tot (N ) ( 2 ) (N ) . In Fig. 3, we plot the TUR product Q (N ) . It can be seen that it approaches a stationary value. Since in the stationary state the data of Q (N ) fluctuates stochastically around a certain value, we introduce Q (N ) τ , i.e., the temporal mean of Q (N ) for times t ≥ τ . This yields a quantity that can be compared to Q from (34) shown as dashed lines in Fig. 3. Note that the value of τ ≥ 10 3 is chosen heuristically based on the observations from Fig. 3. From Tab. 3, it can be seen that Q (N ) τ ranges from 4.16 to 4.58. Hence, for all calculated configurations the TUR product is significantly greater than 2 and thus the numerical calculations support the theoretical prediction from (34) and [6] well. It can be further inferred from Tab. 3 that all the Q (N ) τ 's underestimate the theoretically predicted values. This is due to the above discussed observation that, at least for large L and E, the relative errors of Ψ (N ) (t) 2 and var Ψ (N ) (t) tend to zero, whereas the relative error of ∆s tot (N ) seems to tend to approximately 10%. Therefore, the TUR product is inherently underestimated by the numerical scheme from (17) and the ID discretization of the non-linearity from (13). We have shown that for Ψ (N ) (t) 2 and var Ψ (N ) (t) the predicted scaling Table 3:  forms from (26) and (28), respectively, fit the numerically obtained data for these two quantities very well. Especially for large values of the number of grid-points L, we observed for Ψ (N ) (t) 2 a clear decrease in the relative error between the theoretical predictions and the numerical results, i.e., ∆ 1 → 0 (see Tab. 1). The relative error ∆ 2 of the variance var Ψ (N ) (t) did not depend on L and seems to be solely caused by stochastic fluctuations due to the limited ensemble size E (see Tab. 1). With the above two quantities, both components of the precision ( 2 ) (N ) from (30) were found to follow the predicted scaling forms and thus also the numerically obtained precision behaves as expected (see Tab. 2). In Tab. 2 we have seen for ∆s tot (N ) , that the scaling of the numerical data fits well with the theoretically predicted one from (33). It was observed, however, that even for large L the relative error did not get smaller than roughly 10%. We conclude that this is an inherent issue with the numerical scheme from (17) with the non-linearity according to (13). Further discussion of this point will follow in section 6. For the TUR product we observe that all simulated systems tend to a stationary value for Q (N ) = ∆s tot (N ) ( 2 ) (N ) (see Fig. 3). However, the numerical value is underestimating the theoretically expected one in all cases (see Tab. 3). Nevertheless, the numerical data shows clearly that the TUR product is well above the value of 2 and ranges for our simulations roughly between 4 and 5, which is strong support for the analytical calculations from [6].

Implications of Poor Regularity
As we have already mentioned in section 3, the solution h(x, t) to the KPZ equation from (1) is a very rough function for all times t. The spatial regularity of h(x, t) cannot be higher than that of the solution to the corresponding Edwards-Wilkinson equation, h (0) (x, t), i.e., the KPZ equation with a vanishing coupling constant, λ = 0 in (1) (see, e.g., [29,35,37]). For h (0) (x, t) it can be checked that for all t > 0 where H s denotes the Sobolev space of order s ∈ R of 1-periodic functions, is not a well-defined quantity. Using Hölders inequality for the expectation, one gets L 2 is not well defined either. These two expressions do, however, play an important role in determining the TUR constituents. Hence, some form of regularization is needed to make these expressions well-defined. In [6], this was accomplished by introducing a cutoff Λ of the Fourier-spectrum, i.e., |k| ≤ Λ.
Here, we will follow a different path. Since the solution of the Edwards-Wilkinson equation given by is expected to be a reasonable approximation to the solution of the KPZ equation (1) for λ 1, we approximate the KPZ non-linearity (∂ x h) 2 by (∂ x h (0) ) 2 . The Fourier-coefficients h (0) k are given by where µ k = −4π 2 k 2 (see [6]) and η k is the k-th Fourier-coefficient of η(x, t) = k η k e 2πikx , i.e., the Fourier-series of the KPZ noise from (1). This procedure is equivalent to solving the KPZ equation by a low order perturbation solution with respect to λ, which was performed in [6]. We then replace in 1 0 dx (∂ x h (0) ) 2 and 1 0 dx (∂ x h (0) ) 4 the non-linearity (∂ x h (0) ) 2 with any of its generalized discretizations N where we have defined as the continuum variant of (10). For simplicity, as the operator only acts on x we omit the time t in the above equation. Of course, the expressions from (43), (44) will diverge for δ → 0. The necessary regularization of these expressions is performed by introducing a smallest δ > 0. Both ways of regularization are based on the physical idea of introducing a smallest length scale [28], here in real space and in [6] in Fourier space, so their respective results are directly comparable to one another.

Expected Integral Norms of the Non-Linearity
The expectation of the L 1 -norm of N (γ) δ [h (0) ] from (43) is evaluated for t 1 and δ 1 as where the details of the calculation are given in Appendix A. Similarly, the expectation of the L 2 -norm squared from (44) reads as shown in Appendix B. The expressions in (46) and (47) being divergent for δ → 0 reflect the regularity issues from above.

Approximations of the TUR Constituents
We now establish how the expressions in (43) and (44) are related to the respective constituents of the thermodynamic uncertainty relation.

The Output Functional
Consider the dimensionless form of the KPZ equation from (1) (see, e.g., also [6]). Performing a spatial integration within the boundaries (0, 1) yields with the definition of the output functional Ψ (t) from (21) Due to the periodic boundary conditions the diffusive term vanishes. A subsequent averaging leads to In the NESS, (49) becomes We now approximate the right hand side of (50) by (43) and thus we find with (46) for the output functional in lowest non-vanishing order of λ eff and for t 1

The Total Entropy Production
From [6] we know that in the NESS ∆s tot = σ t holds with σ as the entropy production rate given by (see [6]) where we now approximate the right hand side of (52) by (44). Hence, with (47) we obtain for the entropy production rate for t 1 (53)

The Variance
We have var[ where with h 0 (t) the 0-th coefficient of the Fourier series h(x, t) = k h k (t)e 2πikx for the KPZ solution. The h k may be expanded in terms of the effective coupling constant λ eff and to lowest non-vanishing order it reduces to h k (t) ≈ h (0) k (t), where the latter corresponds to the solution of the Edwards-Wilkinson equation (λ eff = 0). Hence, to lowest non-vanishing order we get from (55) where (42) and the relation η 0 (r)η 0 (s) = δ(r − s) have been used. The second term in (54) is known from (51) and thus gives no contribution to the O(λ 0 eff )-term from (56). However, for completeness we note that the next nonvanishing term in a λ eff -expansion of (Ψ (t)) 2 is O(λ 2 eff ) and the prefactor of λ 2 eff contains the contribution 1/((γ + 1) 2 δ 2 )t 2 , which cancels the second term in (54). This is similar to the continuum case in [6]. Thus, to lowest non-vanishing order in λ eff , the variance is for all 0 ≤ γ ≤ 1 given by in dimensionless form.

The Discrete TUR Product
As the variance from (57) is to lowest order identical to the theoretical predicted one, the TUR product as a function of γ, denoted by Q Since Q (γ) δ is monotonously increasing with γ, the case of γ = 0 represents a lower bound on the TUR product. Hence, the TUR is clearly not saturated as was predicted in [6]. Compared to [6], we here follow an independent path in obtaining the TUR product in (58). Instead of using a Fourier space representation, we derive the TUR product from real space calculations. In particular, we introduce a smallest length scale δ in real space as the regularizing parameter opposed to a cutoff parameter in the Fourier spectrum in [6]. In other words, we work here with the full Fourier spectrum but have to approximate the non-linearity by (45), whereas in [6] we calculated the exact non-linearity, but only on a finite Fourier spectrum with |k| ≤ Λ. Below we connect these two approaches. We note that the divergences of Ψ (t) 2 and ∆s tot are in the present representation in 1/δ 2 for δ → 0, whereas in [6] these expressions diverge in Λ 2 for Λ → ∞. Thus, the analysis above may well be understood as an alternative way of calculating the thermodynamic uncertainty relation.

Comparison of the Approximated TUR Constituents with the Theoretical Predictions
To compare the results of the approximated TUR components in (51) and (53) to the theoretical predictions from [6], we need to express the lattice-spacing δ in terms of the Fourier-cutoff Λ used in [6]. On the interval (0, 1), the latticespacing is given by δ = 1/L, with L the number of grid-points. Using the link between L and Λ from (20) leads to where the last step holds for large enough Λ. Thus, with the theoretical predictions for Ψ (t) 2 and ∆s tot from [6] given by we can calculate the relative deviation ∆ of Ψ (t) δ , respectively. For the expectation of the output functional squared we obtain for Λ 1 with (59), (60) and (51) Analogously, we get for the relative deviation of ∆s tot (γ) δ for Λ 1 and with (59), (61) and (53) With the theoretical prediction from [6], where the last step holds for Λ 1, we calculate the relative deviation of the thermodynamic uncertainty product according to In Tab. 4 we show the relative deviations of all three quantities for some significant values of γ. As can be seen, the overall best result is obtained for Overview of the relative errors of the approximated TUR components from (51) and (53) as well as of the TUR product itself from (58). A negative sign in ∆ indicates that the respective approximated value overestimates the theoretically predicted one and vice versa. γ = 1/2. For all other choices of γ as displayed in Tab. 4, either all the relative errors are greater or, if one of the three errors is chosen to be zero, the two others turn out to be larger than the respective ones for γ = 1/2. In fact, γ = 1/2 minimizes the target function for w = 2. Choosing, e.g., w = 1, 3 results in γ = 0.48, 0.52, respectively. Hence, γ = 1/2 provides in a natural sense a much better approximation than γ = 0, 1. The main purpose of the above analysis was to confirm and explain our key numerical findings from Fig. 1, Tab. 1 and Fig. 2, Tab. 2. Namely, that for γ = 1/2 the error of Ψ (N ) (t) 2 nearly vanishes, while ∆s tot is underestimated by roughly 10% and consequently the TUR product is underestimated by roughly 10% as well (see Tab. 3). These findings are confirmed by the corresponding analytical results in Tab. 4 and (58). Furthermore, we infer from the analysis above that the deviation for ∆s tot (N ) (and thus for the TUR) cannot be reduced by changing the parameters of the numerical scheme like lattice-size δ = ∆x or the time step ∆t. It is instead caused by an intrinsic property of the non-linear operator N with γ = 0, 1/2, 1. Fig. 4 shows the data for Ψ (N ) (t) 2 and ∆s tot (N ) from (26) and (33), respectively. We quantify the significant differences between the respective graphs in Tab. 5. A comparison of the numerically found relative errors ∆ in Tab. 5 to those analytically obtained in Tab. 4 shows very good agreement, which supports the above analysis. As the variance is not dependent on the respective choice of γ (see (57)), which was also reproduced by the numerics, we refrain from explicitly showing this plot as there is no discernible difference in the three graphs. Finally, we show the TUR product Q (N ) for the three different choices of γ in Fig. 5, indicating a clear distinction between the three different discretizations and good agreement with the analytically calculated values from (58) represented by the dashed lines in the plot. Comparison of the predicted scaling factors c 1 (L) and c 4 (L) from (26) and (33) to a L and e L from (35)

Conclusion
We have performed direct numerical simulation of the KPZ equation driven by space-time white noise on a spatially finite interval in order to test the analytical results from [6] regarding the KPZ-TUR in the NESS that were based on a perturbation expansion in Fourier space. Due to the spatial roughness of the solution to the KPZ equation (see section 6), the discretization of the nonlinear term is of great importance. It may be chosen from a set of different variants introduced over the last few decades, which all belong to the so-called generalized discretization N (γ) δ , with 0 ≤ γ ≤ 1 [27]. The numerical data in section 5 was obtained with γ = 1/2, whereas in subsection 6.5 we also used γ = 0, 1, to illustrate the lower and upper bounds of the KPZ-TUR product, respectively. The choice of γ = 1/2 leads to the so-called improved discretization [24] for the KPZ equation. N (1/2) δ is distinguished by the fact that it preserves the continuum steady state probability distribution of h(x, t) [24,27,28]. This implies that also the continuum expression for the total entropy production as derived in [6] remains true in the discrete case (δ > 0). As the limit δ → 0 inherently diverges due to the surface roughness, we believe this feature of N (1/2) δ to be of importance. We have further analytically shown in section 6 and confirmed numerically in subsection 6.5 that the discretization with γ = 1/2 leads to the most accurate approximation of the results in [6], which again highlights the significance of N A central result of this paper, numerically obtained in subsection 5.4 and analytically shown in subsubsection 6.3.4, is that for all choices of γ the TUR product clearly does not saturate the lower bound Q = 2. In particular, we have found as lowest value 4 for γ = 0 and as largest value 6 for γ = 1. Our preferred choice of γ = 1/2 leads to a TUR product of 9/2 (see subsubsection 6.3.4), which is also found within the numerical data in subsection 5.4. This 10% underestimation of the theoretical prediction from [6] is independent of the lattice spacing δ and the time-step size ∆t. By using an idea presented in [29], which consists basically of testing how the discretized non-linearity of a rough SPDE acts on the solution of the corresponding linearized equation, we were able to show analytically that this deviation is an intrinsic property of the N (1/2) δ -operator. Whereas it recovers the correct scaling of Ψ (t) 2 it underestimates the scaling factor of ∆s tot by 10% (see section 6). Furthermore, the analysis in section 6 may be seen as an alternative way, compared to [6], of deriving the KPZ-TUR. We thus conclude that the value 9/2 for the TUR product obtained with γ = 1/2 is the most reliable result that can be achieved by direct numerical simulation of the KPZ equation. Regarding future work, the findings in [24,25,26] lead us to believe that a pseudo spectral simulation of the KPZ equation might yield an even closer approximation to the value 5 as found in [6] than the one obtained in this paper by direct numerical simulation.
To obtain the result in (46), we define then the expression in (45) may also be written as Using (68), (43) for 0 ≤ γ ≤ 1 reads The first two terms in (69) are calculated via where Denoting by (·) the complex conjugate, the right hand side of (70) is evaluated as follows where we have used in the second step that 1 0 dx e 2πi(k−l)x = δ k,l with δ k,l the Kronecker symbol. The third step employs the two-point correlation function of h 2 µ k (p + q) 2 δ 2 = k∈Z\{0} 1 − cos 2πk(p + q)δ (2πk(p + q)δ) 2 .
(73) With the substitution x = 2πk(p + q)δ and δ 1, we may rewrite (73) as The integral in (74) can be evaluated by either using the residue theorem or by employing an adequate CAS, and yields π/2. Hence, the expression of (70) is given for t 1 and δ 1 by 2 cos x − cos 2x − 1 where we have substituted x = 2πkδ with δ 1 and used the value of the integral in (74). Combining (75) and (76) gives (46).

B Expectation of the
To ease the calculation of (44), let us first rewrite (45) in the following way, Consequently, where we have used (71). The four point correlation function can be evaluated via Wick's theorem and with Π k,l (t, t ) ≡ e µ k t+µ l t 1 − e −(µ l +µ l )(t∧t ) µ k + µ l , µ k = −4π 2 k 2 as above [6]. With (80) and (81), the expression in (79) becomes where we have substituted x = 2πlδ for δ 1 and used (81) for t 1. Next, we will calculate where we have again used Wick's theorem and (80) with (81) as well as an index shift k − l ↔ l to obtain the prefactor of two in the second term. The first term in (84) reads for t 1 l∈Z\{0} Π l,l (t, t) C since the second sum in (85) has the same form like the one in (76). The second term in (84) may be evaluated with (81) for t 1 by substituting x = 2πlδ for δ 1 according to k∈Z l∈Z\{0,k} Π l,l (t, t)Π k−l,k−l (t, t)C Here, we used again an index shift k − l ↔ l to obtain the factor of two in front of the first sum of (87). The first term in (87) has, after inserting (81) for t 1 and substituting x = 2πlδ the same form as the second sum in (85) and thus vanishes. The second term in (87) becomes for t 1 k∈Z l∈Z\{0,k} Π l,l Π k−l,k−l C where we substituted in the second step x = πlδ. Hence, combining (83), (85), (86) and (88) leads to which is the result given in (47).