Serge Meynet

Alexis Barge

Vincent Moureau

Guillaume Balarac

Ghislain Lartigue

Abdellah Hadjadj

P R Prandtl

ROUGHNESS-RESOLVED LES OF ADDITIVE MANUFACTURING-LIKE CHANNEL FLOWS

Additive manufacturing (C)HX (Compact) heat exchanger ES Effective slope f Fanning friction factor H

INTRODUCTION

In aeronautical engines, heat exchangers have a fundamental role especially in controlling the temperature of the liquid fuel and oil streams. The mass of heat exchangers is also critical to reduce the total weight of the engine. The design of compact and efficient heat exchangers (CHX) is therefore of paramount importance for the future aeronautical engines. Additive manufacturing (AM) offers the opportunity to explore new designs and boost the optimization of CHX. [START_REF] Carozza | Heat exchangers in the aviation engineering[END_REF]. For instance, as shown by Saltzman et al. [START_REF] Saltzman | Design and evaluation of an additively manufactured aircraft heat exchanger[END_REF], they obtained an improvement around 10 percent compared to heat exchangers built with traditional manufacturing.

Nonetheless, surface roughness generated with AM can be substantially larger than with traditional manufacturing techniques. Surface state of metal AM test samples analysed in several studies have indeed underlined this key point [START_REF] Saltzman | Design and evaluation of an additively manufactured aircraft heat exchanger[END_REF][START_REF] Cabanettes | Topography of as built surfaces generated in metal additive manufacturing: A multi scale analysis from form to roughness[END_REF][START_REF] Ventola | Rough surfaces with enhanced heat transfer for electronics cooling by direct metal laser sintering[END_REF][START_REF] Simonelli | Effect of the build orientation on the mechanical properties and fracture modes of slm ti-6al-4v[END_REF][START_REF] Snyder | Build direction effects on additively manufactured channels[END_REF][START_REF] Stimpson | Scaling roughness effects on pressure loss and heat transfer of additively manufactured channels[END_REF]. Moreover the wall roughness is not isotropic in space and varies according to the main direction of printing [START_REF] Ventola | Rough surfaces with enhanced heat transfer for electronics cooling by direct metal laser sintering[END_REF][START_REF] Simonelli | Effect of the build orientation on the mechanical properties and fracture modes of slm ti-6al-4v[END_REF]. Some patterns called welding tracks can appear on the surface [START_REF] Cabanettes | Topography of as built surfaces generated in metal additive manufacturing: A multi scale analysis from form to roughness[END_REF]. This typical roughness has a significant impact on the performances in terms of pressure drop and heat transfer capacity [START_REF] Saltzman | Design and evaluation of an additively manufactured aircraft heat exchanger[END_REF][START_REF] Snyder | Build direction effects on additively manufactured channels[END_REF][START_REF] Stimpson | Scaling roughness effects on pressure loss and heat transfer of additively manufactured channels[END_REF]. In addition, an important gap between performances obtained from simulations and/or empirical correlations and those measured a posteriori exists.

Many papers have been devoted to the study of turbulent flow over different roughness types and surfaces experimentally [START_REF] Jiménez | Turbulent flows over rough walls[END_REF] as well as numerically [START_REF] Piomelli | Recent advances in the numerical simulation of roughwall boundary layers[END_REF][START_REF] Napoli | The effect of the slope of irregularly distributed roughness elements on turbulent wall-bounded flow[END_REF][START_REF] Busse | Direct numerical simulation of turbulent flow over a rough surface based on a surface scan[END_REF]. The presence of surface roughness causes an overall thickening of the boundary layer in addition to the modification of the streamwise structures. Secondary flow motions can also be observed due to spanwise inhomogeneities on the surface and this induces a better flow mixing. In addition, the transition between laminar to turbulent regimes appears earlier than over smooth surfaces. Concerning the mean velocity profile, a downward shift compared to the log-law profile is observed [START_REF] Clauser | Turbulent boundary layers in adverse pressure gradients[END_REF][START_REF] Hama | Boundary layer characteristics for smooth and rough surfaces[END_REF]. The latter is commonly called the roughness function (∆U + ). Roughness has also an impact on the Nusselt number since the increase of surface exchange enhances heat transfer [START_REF] Stimpson | Roughness effects on flow and heat transfer for additively manufactured channels[END_REF].

Despite a better understanding of the flow behavior over inhomogeneous rough surfaces, some attempts have been performed to get a universal correlation for evaluating the equivalent sand-grain roughness height k s [START_REF] Flack | Review of hydraulic roughness scales in the fully rough regime[END_REF][START_REF] Forooghi | Toward a universal roughness correlation[END_REF][START_REF] Stimpson | Scaling roughness effects on pressure loss and heat transfer of additively manufactured channels[END_REF]. This parameter k s can be expressed as a function of the statistical roughness parameters such as the arithmetic average height (S a ), the root-mean-square height (S q ), the skewness (Sk) and the kurtosis (Ku). The two latter correspond respectively to the third and the fourth moments of the height distribution of the surface. Studies have underlined the predominance of the effective slope (ES), defined as the average slope of the surface height along the streamwise direction, in predicting roughness effects [START_REF] Napoli | The effect of the slope of irregularly distributed roughness elements on turbulent wall-bounded flow[END_REF].

S a = 1 L x L y Lx Ly | h(x, y) | dxdy (1)
S q = 1 L x L y Lx Ly h(x, y) 2 dxdy (2) Sk = 1 S q 3 L x L y Lx Ly h(x, y) 3 dxdy (3) Ku = 1 S q 4 L x L y Lx Ly h(x, y) 4 dxdy (4) ES = 1 L x L 0 | dh dx | dx (5)
where h(x, y) is the surface height, the mean plane is supposed to be at h = 0, L x and L y are respectively the streamwise and spanwise lengths.

A well-known threshold equal to ES ≈ 0.35 delimits two regimes: i) the sparse or waviness regime below this limit and ii) the roughness or dense one [START_REF] Macdonald | Turbulent flow over transitionally rough surfaces with varying roughness densities[END_REF]. Nevertheless, such universal correlation that would cover the two regimes has not been found up to now. Furthermore, few articles propose numerical simulation of AM surfaces.

In the aim of building new rough-wall models dedicated to AM surfaces, the creation of a roughness-resolved Large-Eddy Simulation (RRLES) database of representative channel flows emerges as a worthwhile process. The present paper describes a methodology to perform parametric RRLES from a limited set of surface parameters to control the surface height distribution and the streamwise and spanwise spatial auto-correlation functions. The first section of the article describes the numerical setup, configurations and the used methodology for performing RRLES. The post-processing and obtained results are presented and discussed in the second section.

NUMERICAL SETUP

In this section, chosen configurations and methodology are addressed. In order to ensure periodic channels, recycling conditions are applied on velocity and temperature. This process is particularly exposed here.

Configurations and Meshes

First of all, rough surface (RSG) and body-fitted unstructured mesh (RRMG) generators have been developed [START_REF] Meynet | Automatic surface and volume mesh generation for roughness-resolved les of additivemanufacturing heat exchangers[END_REF]. The former creates triangulated surfaces stored in the STL format that are used by the latter to compute level set functions. The level set functions are then employed to refine an existing 3D unstructured mesh and to create a conforming body-fitted mesh with controlled cell size, quality and size gradation. The resulting meshes are suitable for performing RRLES with higher-order finite-volume schemes. The RRMG has been developed within the YALES2 code. YALES2 is an unstructured low-Mach number code developed at CORIA for the Direct Numerical Simulations and Large-Eddy Simulations in complex geometries [START_REF] Moureau | Design of a massively parallel cfd code for complex geometries[END_REF]. It heavily relies on the parallel volume and surface mesh adaptation.

The RSG/RRMG are utilized to generate rough planar channels but may also be applied to more complex shapes such as cylindrical tubes or plate/fin configurations. RRLES are conducted for periodic channels of size 8H x 3H x 2H in the streamwise (x), spanwise (y) and crosswise (h) directions with H the half height of the channel which is equal to 1.0 mm. This domain size was confirmed to be sufficient to have a negligible impact of the periodic recycling on the wall friction [START_REF] Kuwata | Direct numerical simulation of turbulence over systematically varied irregular rough surfaces[END_REF]. In addition, the half height was selected in order to be close to hydraulic diameters and channel heights encountered for some AM experiments as the L-2x-In case in Stimpson et al. [START_REF] Stimpson | Roughness effects on flow and heat transfer for additively manufactured channels[END_REF].

Concerning the chosen configurations, three different cases representing two printing directions plus a streamwise and spanwise isotropic case have been selected for analysis. A smooth channel is also considered as a reference. Roughness parameters which are targeted in this article are the root-mean-square height (S q ), the height distribution skewness (Sk) and kurtosis (Ku) and the effective slope (ES). While the surface height distributions are the same in the three cases, the effective slope is different due to the different space autocorrelation functions.

In addition of printing direction, the idea is to generate surfaces which mimic roughness height distribution encountered in AM. An overall range of roughness parameters values can be found in the literature even though there are differences between upward-and downward-facing surfaces [START_REF] Saltzman | Design and evaluation of an additively manufactured aircraft heat exchanger[END_REF][START_REF] Cabanettes | Topography of as built surfaces generated in metal additive manufacturing: A multi scale analysis from form to roughness[END_REF][START_REF] Ventola | Rough surfaces with enhanced heat transfer for electronics cooling by direct metal laser sintering[END_REF][START_REF] Simonelli | Effect of the build orientation on the mechanical properties and fracture modes of slm ti-6al-4v[END_REF][START_REF] Snyder | Build direction effects on additively manufactured channels[END_REF][START_REF] Stimpson | Scaling roughness effects on pressure loss and heat transfer of additively manufactured channels[END_REF]. For this purpose, values of the latter for S q , Sk and Ku are set respectively to 20 µm, around 0.2 and 4.0. Final values are exposed in Tab. 1 and probability density functions of height are plotted in Fig. 2 With the RRMG, conforming body-fitted meshes are obtained from these geometries. A balance on the element count has to be reached in order to resolve sufficiently the flow while minimizing the cost. In addition, the meshing process has to keep the topology intact. The mesh generation process starts from a cartesian grid which is tessellated and adapted numerous times.

The initial numbers of elements in each direction (N x , N y , N h ) for the cartesian grid are fixed in order to initially get isotropic cells of 20 µm. The cell size gradient is equal to 0.1 and the cell size on the rough surface is set to ∆y w = 7 µm. The final number of mesh elements for RRLES is about 130 million cells.

All characteristics and performances of meshing process are summarized in tables Tab. 2 and Tab. 3. The mesh for the smooth case is a stretched cartesian grid with the dimensionless wall grid spacing of ∆ + x = 13.7, ∆ + y = 5.5, ∆ + h = [0.5, 7.3] in accordance with operating conditions.

For each configuration, the mesh size is fixed for all the tested Reynolds numbers. The mesh resolution for the present study is discussed in the appendix. It must be noted that the dimensionless roughness S a /∆y w is only larger than two, which may seem insufficient for the wall resolution. However, with bodyfitted grids, this latter has to be sufficient for capturing the wall roughness (curvature and position as illustrated in Fig. 3) and the wall-normal velocity gradient.

Parameters

Values L x , L y , L z (mm) 8.0; 3.0; 2.0 Initial N x , N y , N z 400 ; 150 ; 100 Cell size on STL 7 µm Max cell size gradient 0.1

Statistics and Performances

For collecting statistics, each RRLES is split into two steps. The initialization step is performed during a given transient time and then statistics are accumulated. For each step, the durations are summarized in Tab. 4. For the second step, the dimensionless time t + defined by t + = tuτ H is around 20 in average. Interestingly, due to inhomogeneities of the surface, the time to establish a fully-developed turbulent channel flow is reduced compared to smooth cases.

RRLES CPU costs are presented in Tab. 5. For clarity, the CPU hours are averaged over all Reynolds numbers cases for each configuration.

Re

Numerics and Models Working hypotheses

In each case, incompressible flow simulations are performed. The chosen target bulk Reynolds number range is the fully developed turbulent flow. This is why, RRLES are performed at Re = 5 000, Re = 8 000, Re = 17 000 and Re = 25 000. The fluid kinematic viscosity is set to ν = 1.517 10 -5 m 2 /s and the maximal CFL number used is equal to 0.8. The WALE subgrid scale model is retained as it is widely used for LES of boundary layer flows [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF]. A fourth-order central finite-volume scheme is used and the four-step fourth-order scheme TFV4A is applied for velocity and scalar transport prediction [START_REF] Kraushaar | Application of the compressible and low-Mach number approaches to Large-Eddy Simulation of turbulent flows in aeroengines[END_REF].

Methodology, boundary conditions

The automatic generation of fully periodic channels is very challenging and ensuring periodic meshes is complex for unstructured grid. The periodicity is though necessary to reach statistical convergence of the flow. Instead of imposing strict periodic boundary conditions in the streamwise and spanwise directions, a Lagrangian recycling method has been developed. This method does not require to add a body force as source term in Navier-Stokes momentum equation to compensate the wall friction. This approach can also be used as a precursor in spatially developing boundary layers by increasing L out . The idea behind this recycling method is simply to use a time-varying inlet boundary condition with a velocity coming from a distant plane in the channel. This 3step recycling process is based on passive Lagrangian particles as illustrated in Fig. 4. A flow rate is imposed at the inlet and for optimizing performances, the recycling is done every N time steps, N > 1. A linear interpolation in time is performed between two recycled planes. No rescaling is applied on velocity profile as such profile is unknown a priori. Thus, contrary to what is proposed in Xiao et al. for instance, no target or corrections on velocity field are applied as input of this method. In this article, we select N = 20 as it gives the better compromise between performances and velocity interpolation errors. The complete process is described here:

1. Lagrangian particles are created at the inlet of the grid (blue spheres in Fig. 4) and are translated to the recycling plane. 2. The target field at the recycling plane, for instance velocity, is interpolated for translated particles (red spheres in Fig. 4) 3. This is the relocation step: particles are transferred at the inlet and the field at this location is updated.

The computational domain is then defined with a recycling zone and a buffer area to avoid any influence of the outlet boundary treatment onto the recycled velocity. Thus the location of the recycling plane is primordial and should be set at a given distance from the outlet within the domain. For two parallel planes, this distance was found to be at least equal to the length between both planes and in our cases, the latter corresponds to 2H. Indeed, perturbations of the velocity field due to the outlet boundary condition may be recycled and injected at the inlet if L out is below this threshold.

In order to decrease interpolation margin errors between the inlet and the recycling plane, both locations should be the equivalent at the wall. The generated rough surface and conformal mesh both respect this periodicity condition only in terms of wall topology. A no-slip wall boundary condition is applied on rough planes and the other sides in the spanwise direction are considered as slip walls. This methodology is also applied in the smooth configuration. Validation of the numerical methods may be found in the appendix.

Finally, from STL generation to post-processing, all calculation steps are integrated into a workflow tool. This allows to manage automatically series of runs on a distant super-computer.

SIMULATION ANALYSIS

This section details the different tools used for the subsequent analysis of the RRLES.

Non-dimensional velocity and temperature

For scaling velocity and temperature profiles, a calculated effective distance introduced by Kuwata & Kawaguchi [START_REF] Kuwata | Direct numerical simulation of turbulence over systematically varied irregular rough surfaces[END_REF] is used. Indeed, due to irregularities of the surface height, this kind of distance is not straightforward to determine as in a smooth-wall case. The effective distance is defined in Eq. 6 with h w the minimal height of the surface. The variable φ corresponds to the x-y plane porosity which is the ratio between the x-y plane surface occupied by the fluid and the total x-y plane area.

h e = h hw φ dh (6)
Computation of the friction velocity u τ (Eq. 7) is based on the difference between average pressure at the inlet and at the recycling plane as exposed in Fig. 6. The shear Reynolds number is then calculated as Re τ = uτ H ν . and the quantity h + e is defined as h + e = heuτ ν .

u τ = H ⟨p i ⟩ -⟨p r ⟩ ρL (7)
Concerning the friction factor, the Fanning definition

f = 2 uτ U b 2
is used with U b the bulk velocity. All these quantities are monitored at each iteration in the LES simulation.

For the analysis of heat transfer, a passive filtered scalar Z is used. This latter can be defined as Z = T -Tp T∞-Tp with T p the temperature imposed at a wall and T ∞ the bulk temperature. This scalar is considered passive and this hypothesis is valid if the temperature difference has no significant impact on the density, which is assumed here. This is why, the temperature can be replaced by this dimensionless scalar. The equation for this scalar is the following:

∂ Z ∂t + ∇ • (ū Z) = ∇ • (D z ∇ Z) (8)
Thus, Z = 1 is imposed at the upper wall and Z = 0 everywhere including the bottom wall. The laminar Prandtl number P r of this scalar is set to P r = 0.71 and the turbulent Prandtl number is equal to unity. The diffusivity D z includes the molecular and turbulent diffusivities.

Budget equations

In the present RRLES, different turbulence budget equations have been computed, dumped and stored. This is the case for the mean kinetic equation (Eq. 9). Quantities denoted with the bracket ⟨⟩ symbol are time-averaged and u ′ corresponds to the fluctuating component of velocity within the Reynolds decomposition. The quantity ν t refers to the turbulent viscosity which is modeled through the WALE model.

Mean Kinetic Equation (MKE)

1 2 ∂⟨u i ⟩⟨u i ⟩⟨u j ⟩ ∂x j 1 + ∂⟨u j ⟩⟨p⟩ ∂x j 2 + ∂ ⟨u i ⟩ τ SGS ij ∂x j 3 -2 ∂ ⟨u i ⟩ ⟨νS ij ⟩ ∂x j 8 + ∂ ⟨u i ⟩ ⟨u ′ i u ′ j ⟩ ∂x j 4 = ⟨S ij ⟩ τ SGS ij 5 + ⟨S ij ⟩ ⟨u ′ i u ′ j ⟩ 6 -2 ⟨S ij ⟩ ⟨νS ij ⟩ 7 (9)
with

τ SGS ij = -2ν t S ij and S ij = 1 2 ∂ui ∂xj + ∂uj ∂xi .

RESULTS AND DISCUSSION

In this section, the wall-unit velocity and temperature profiles as well as the momentum and kinetic energy balance equations described in the previous section are analyzed. This analysis should provide a better understanding of the impact of the wall roughness on the flow and especially the effective slope parameter.

Impact on Turbulence

A first qualitative analysis can be done by visualizing the vortices which are generated over the rough surface. The same Q-criterion iso-contours are plotted for the three different cases in Fig. 7 at a Reynolds number of 8 000. At a first glance, it is clearly noticed that the parallel case promotes elongated vortices in the flow direction while the two other cases also feature small vortices trapped between valleys. These trapped vortices are more coherent in the front case due to the orthogonal wavy pattern of the wall. Even if the roughness height distribution is the same, important differences are visually perceptible.

Mean streamwise velocity profiles are plotted in Fig. 8 for the three cases at Reynolds numbers of 8 000 and 17 000 as well as the smooth wall configuration. The impact of ES is clearly visible as the velocity profile is shifted downward with increasing ES. Differences are also noticeable between isotropic and anisotropic surfaces. Indeed, if the ES in isotropic case is lower than the front configuration, the mean velocity near the wall has a steeper increase for the isotropic case. These results indicate that the orientation of the wall roughness has a nonlinear impact on the mean velocity in the channel. Moreover, increasing the Reynolds number tends to shift downward the profile in the log-law region as expected. Concerning, the fluctuating velocity, streamwise and spanwise components are plotted for Reynolds numbers Re = 8 000 Re = 17 000 in Fig. 9 and 10, respectively. For the streamwise fluctuating velocity u rms , and considering the two anisotropic cases, the ES has an impact on the maximum value and this latter is shifted to the flow stream when changing the roughness direction. This peak is also less pronounced for higher ES value. As for the mean velocity, the isotropic case has a slightly different behavior and while its ES value is in between the two anisotropic cases, the u rms profile is not in between the two of the anisotropic cases confirming the non-linear behavior of the rough wall with ES.

For spanwise velocity fluctuations, the maximum of the profile is less affected by the geometry but its location follows the same trend as the streamwise fluctuations. Increasing the Reynolds number to 17 000 tends to shift the peaks to the channel center, to reduce the sharpness of the peaks for both velocity fluctuations and to slightly increase the maximum value of peaks for spanwise fluctuations.

The reduction of the anisotropy of the velocity fluctuations observed in each case in comparison with the smooth configuration is consistent with previous studies [START_REF] Shafi | Anisotropy of the reynolds stresses in a turbulent boundary layer on a rough wall[END_REF][START_REF] Piomelli | Recent advances in the numerical simulation of roughwall boundary layers[END_REF]. It has been shown that this is caused by the redistribution of the mean roughness wake kinetic energy into turbulence, generating additional normal and spanwise stresses [START_REF] Yuan | Roughness effects on the reynolds stress budgets in near-wall turbulence[END_REF]. The streamwise turbulence is also mainly converted into wake energy as the streaks are destroyed by the roughness elements. Increasing the Reynolds number contributes to intensify this phenomenon. The presence of the wake is slightly visible on the mean velocity profiles when they become negative in the near-wall region for the isotropic and front cases. Mean kinetic energy balance plotted in Fig. 11 shows that this region is dominated by the work of the pressure drag against the flow direction. Those observations are no more valid for the parallel case since the surface topology does not create significant re-circulation zones, letting the viscous drag to be dominant. Thus, the decreasing of the turbulence anisotropy and the pressure loss are less pronounced for this case.

For each Reynolds number and all roughness types, the roughness velocity function ∆U + is studied. The latter is evaluated at h + e = 100 and compared to the boundary-layer log law. Results are summarized in Tab. 6. The first remark is that ∆U + is amplified with higher Reynolds number. Values are mainly higher in the front case except at Re = 25 000 where the isotropic case has a larger value. These results point out that the additive manufacturing direction has a strong impact on the mean velocity profile although roughness height distributions are the same. Concerning the effect of the effective slope on pressure loss, the friction factor is investigated. The Colebrook correlation [START_REF] Colebrook | Turbulent flow in pipes, with particular reference to the transition region between the smooth and rough pipe laws[END_REF] (Eq. 10) is taken as a reference in order to compare friction factor values obtained in the present RRLES database to the Moody diagram. The definition used in this correlation is Darcy's one and the relationship between Darcy and Fanning friction factors is f Darcy = 4f . Therefore, in Fig. 12, friction factors from the Colebrook equation are divided by a factor of 4 for getting comparative results. Two relative roughness S a /D h values (S a /D h = 10 -2 and S a /D h = 5.10 -3 ) for Colebrook formula are plotted as the relative roughness of the present cases is intermediate:

S a /D h = 7.8 10 -3 . 1 f Darcy = -2 log S a 3.7D h + 2.51 Re f Darcy (10)
As expected, for each Reynolds number, friction factor is higher with increasing ES. For the front and isotropic cases, values are above the Colebrook correlation ones and tend to slightly increase with the Reynolds number. However, the trend for parallel case is reversed: the curve is below Colebrook expectations and decreases. There is a clear distinction in the friction factor behavior when varying the printing direction. Therefore this correlation is inadequate for describing the tendency and evaluating a priori the friction factor. It also highlights again the need for better correlations for this type of roughness. fluctuations observed close to the wall have wave lengths similar to those of the wall roughness. This is particularly visible in the front and isotropic cases. In the parallel case, the chosen cut cannot represent the wall roughness as the roughness features are aligned with the cut plane. As a result, less short wave lengths are visible in this cut close to the wall.

In order to quantify if temperature transport by turbulence is enhanced by wall roughness or not, mean temperature profiles are plotted for two different Reynolds numbers in Fig. 14. At Re = 8 000 (Fig. 14a), the temperature transport in the boundary layer is enhanced by wall roughness as the mean temperature gradient increases towards a constant value in the whole channel height. The three rough cases have a slightly different behavior close to the wall. It seems that the isotropic case is the configuration that is the closest to the smooth case. For the two anisotropic cases, the presence of coherent structures in the roughness valleys seems to slightly increase the temperature wall gradient and thus the wall heat transfer.

For Re = 17 000, the conclusions are different and more intuitive. Indeed, the parallel roughness and smooth cases have almost the same temperature profiles as it would be expected due to the alignment of the roughness elements with the flow. The front and isotropic cases have also the same temperature profile with a clear shift toward the constant gradient imaginary profile. In theses cases, the wall roughness promotes temperature transport by turbulence. For this higher Reynolds case, the dependence of the wall heat transfer with the effective slope is better recovered. Again, the different non-linear behaviors of the temperature with both the effective slope and the Reynolds number illustrate the challenge of finding a Nusselt number correlation for rough walls.

CONCLUSIONS

This paper presents roughness-resolved Large-Eddy Simulations that are representative of the flow obtained in additive-manufactured heat exchangers. The aim of these simulations is to provide a rich database that will ease the development and assessment of rough-wall models. The most challenging task in the building of this database is the generation of representative rough surfaces and conformal unstructured meshes. Three different configurations of parallel rough plane channels with the same roughness height distribution but different effective slopes have been chosen and modeled. The impact of the effective slope parameter, which is directly linked to the alignment of the wall roughness with the flow, has a strong impact on the flow topology, velocity and temperature profiles, as expected. In these cases, the existing empirical correlations find their limits and new correlations are needed.

acknowledgment

The STREAM project has received funding from the Clean Sky 2 Joint Undertaking (JU) under grant agreement No 865378. The JU receives support from the European Union's Horizon 2020 research and innovation program and the Clean Sky 2 JU members other than the Union. This article reflects only the authors' view, the Joint Undertaking is not responsible for any use that may be made of the information it contains. This work was granted access to the HPC resources of CINES/TGCC/IDRIS under the allocation 2020-A0092B06880 and 2021-A0112B06880 made by GENCI and of CRIANN under the allocation 2012006.

Appendix: Numerical methods assessment

We validate our numerical methods by reproducing the classical Re τ = 180 pressure-gradient driven periodic-channel flow DNS test case from [START_REF] Kim | Turbulence statistics in fully developed channel flow at low reynolds number[END_REF] (refered to as KMM hereafter). The geometry is L x = 4πH, L y = 2πH, L h = 2H with the same RRLES direction denomination and the dimensionless mesh resolution is ∆ + x = 8.6, ∆ + y = 8, ∆ + h = [0.38, 5.3]. A first test (T1) is done with the periodic condition methodology and a second one (T2) uses the recycling boundary condition. For the latter, we set the recycling plane at a distance 2H above the outlet in the longitudinal direction and the input flow rate is chosen to impose the same bulk velocity measured in [START_REF] Kim | Turbulence statistics in fully developed channel flow at low reynolds number[END_REF]. Results are summarized in Fig. . We discuss now the RRLES mesh size used in this paper. Fig. .16 compares the results obtained from the Re = 17 000 isotropic configuration with the mesh described in Tab. 2 (denoted M2) and those obtained with a coarser mesh (M1) whose cell size is 10 µm. It is seen that no significant differences may be observed between the two simulations, both on velocity first-and second-order moments. This suggests that results will be probably weakly improved with a mesh finer than M2, for the analyzed moments. Thus, we consider that the dimensionless wall normal resolution y + obtained with the mesh M2 for this case is a satisfactory reference. Due to the surface alternating of peaks and valleys, making the wall normal mean velocity gradient to vary in space, it is more relevant to observe the statistical distribution of the wall resolution than its mean value alone to evaluate the mesh quality. Fig. .17 shows the probability density function of y + obtained for the considered Reynolds numbers. Obviously, the distributions for the two lowest Reynolds numbers are satisfactory. On the other hand, the distribution for Re = 25 000 is larger than our limit. Nevertheless, it is similar to the one obtained for Re = 17 000 with M1. According to previous results, we can consider that M2 is adequate for the simulations performed in this paper. Analysis for the front and parallel surfaces (not shown here) leads to the same conclusions.

Figure 2 :

 2 Figure 2: PROBABILITY DENSITY FUNCTION OF HEIGHT FOR THE THREE CASES

  A) FRONT CASE (ES = 0.29) B) PARALLEL CASE (ES = 0.06) C) ISOTROPIC CASE (ES = 0.24)

Figure 7 :

 7 Figure 7: Q-CRITERION ISO-CONTOUR COLORIZED BY VELOCITY AT RE = 8 000

Figure 8 :

 8 Figure 8: VELOCITY PROFILES (RE = 8 000 CONTINUOUS LINE, RE = 17 000 DASHED LINE)

Figure 11 :

 11 Figure 11: MEAN KINETIC BUDGET TERMS FOR RE = 8000 NORMALIZED BY u 4 τ /ν FROM SMOOTH CASE. 1: -, 2: -, 4: -, 6: -, -7: -, 8: -. U I 0 , U F 0 REFER TO MEAN VELOCITY CROSS-ZERO SCALES FOR ISOTROPIC AND FRONT CASES.

  CASE (ES = 0.24)

Figure 13 :

 13 Figure 13: INSTANTANEOUS TEMPERATURE FIELD AT RE = 8 000

  15. Good agreement between tests and references shows that our numerical schemes and recycling boundary condition are appropriate to perform infinite channel flow simulations.

Figure . 15 :

 15 Figure .15: VELOCITY MEAN AND FLUCTUATION PROFILES IN THE CHANNEL FLOW TEST CASE.

  A) PERPENDICULAR PRINTING DIRECTION TO THE FLOW (ES = 0.29) B) PARALLEL PRINTING DIRECTION TO THE FLOW (ES = 0.06) C) ISOTROPIC ROUGHNESS CASE (ES = 0.24)

Figure 1 :

 1 Figure 1: GEOMETRY OF CONSIDERED CONFIGURATIONS

21

 21

Figure 3 :Figure 4 :

 34 Figure 3: ZOOM-IN OF THE MESH ON THE UPPER SURFACE FOR ISOTROPIC CASE (ES = 0.24)

Figure 5 :Figure 9 :

 59 Figure 5: COMPUTATIONAL DOMAIN SPLIT INTO TWO AREAS: THE RECYCLING AND THE BUFFER ZONES

Figure 10 :

 10 Figure 10: FLUCTUATION OF SPANWISE VELOCITY (RE = 8 000 CONTINUOUS LINE, RE = 17 000 DASHED LINE)

Figure . 17 :

 17 Figure .17: DIMENSIONLESS MEAN WALL RESOLUTION PDF FROM RRLES ISOTROPIC CASES.

Table 1 :

 1 for the three geometries. ROUGHNESS PARAMETERS OF CHOSEN GEOMETRIES

		ES	S a (µm) S q (µm) S h (µm) Sk	Ku
	Front	0.29 15.6	20.0	187	0.21 3.92
	Parallel 0.06 15.6	20.0	181	0.21 4.01
	Isotropic 0.24 15.6	20.0	216	0.21 3.90

Table 2 :

 2 MESHING CHARACTERISTICS

		ES	Cells	CPUs	CPU hours /CPU RAM
	Front	0.29 130.1M 1024 4096h 294Mb
	Parallel 0.06 129.7M 1024 4198h 292Mb
	Isotropic 0.24 122.9M 1024 5307h 279Mb

Table 3 :

 3 MESHING PERFORMANCES

Table 4 :

 4 TIME ACCUMULATION FOR STATISTICS (FTT: FLOW-THROUGH TIME)

		Init run FTT Stats run FTT
		(ms)	Init	(ms)	Stats
	5000	0.7	8.8	4.0	50.6
	8000	0.5	10.7	2.4	51.2
	17 000	0.3	12.7	1.2	51.0
	25 000	0.02	12.7	0.8	50.7
		ES	Cells	CPUs	CPU hours
	Front	0.29 130.1M	560	32 800h
	Parallel 0.06 129.7M	560	18 500h
	Isotropic 0.24 122.9M	560	32 500h

Table 5 :

 5 RRLES PERFORMANCES

Table 6 :

 6 ROUGHNESS FUNCTION VALUES EVOLUTION

	∆U +	Re 5 000 8 000 17 000 25 000 Re Re Re
	Front	5.90	7.84	11.1	13.9
	Parallel	0.41	0.69	1.34	2.89
	Isotropic 3.10	5.48	10.2	15.0

A) RE = 8 000 B) RE = 17 000 The roughness influence on heat transfer is also studied with this RRLES database. Indeed, modifications to the turbulent boundary layer induced by the roughness inhomogeneities can modify the temperature boundary layer and potentially the mean temperature field. As explained above, a passive scalar Z is transported and averaged (Eq. 8). Instantaneous Z field for the three rough cases are illustrated in Fig. 13. This figure shows that the temperature


Dian Wu
email: dian.wu@epfl.ch

Riccardo Rossi
email: riccardo.rossi@epfl.ch

Filippo Vicentini
email: filippo.vicentini@epfl.ch

Giuseppe Carleo
email: giuseppe.carleo@epfl.ch

From Tensor Network Quantum States to Tensorial Recurrent Neural Networks

We show that any matrix product state (MPS) can be exactly represented by a recurrent neural network (RNN) with a linear memory update. We generalize this RNN architecture to 2D lattices using a multilinear memory update. It supports perfect sampling and wave function evaluation in polynomial time, and can represent an area law of entanglement entropy. Numerical evidence shows that it can encode the wave function using a bond dimension lower by orders of magnitude when compared to MPS, with an accuracy that can be systematically improved by increasing the bond dimension.

I. INTRODUCTION

Tensor networks (TN) have been extensively used to represent the states of quantum many-body physical systems [1][2][3]. Matrix product states (MPS) are possibly the simplest family of TN, and are suitable to capture the ground state of 1D gapped Hamiltonians [4,5]. They can be contracted in polynomial time to compute physical quantities exactly, and optimized by density matrix renormalization group (DMRG) [6] when used as variational ansätze. More powerful TN architectures that cannot be efficiently contracted in general have been proposed later, notably projected entangled pair states (PEPS) [7]. They are of interest mainly because they can describe area-law states in dimensions D > 1 [8]. In practice, they are usually approximately contracted and optimized with low-rank truncation [9][10][11]. They can also be optimized in the variational Monte Carlo (VMC) framework [12][13][14][15] using stochastic sampling rather than exact contraction.

In recent years, a growing trend has been to use TN in machine learning tasks such as supervised [START_REF] Stoudenmire | Advances in Neural Information Processing Systems[END_REF][START_REF] Cheng | [END_REF] and unsupervised [18,19] learning. It is of particular interest to use them as generative models [20,21], where it is possible to draw perfect samples with tractable probability densities. Considering that many generative models are based on recurrent neural networks (RNN), there have been attempts to use tensor operations in RNN to improve their expressivity [START_REF] Sutskever | Proceedings of the 28th International Conference on Machine Learning[END_REF][START_REF] Irsoy | 3rd International Conference on Learning Representations[END_REF]. This kind of architectures has shown leading performance when applied to quantum physical systems [START_REF] Hibat-Allah | [END_REF].

Soon after the original proposal of neural quantum states (NQS) [25], several attempts have been made to link them to ground states of well-known Hamiltonians and TN quantum states. While an arbitrary Hamiltonian might require an NQS with a depth growing exponentially with the system size [26], several Hamiltonians with desirable topological properties have compact representations in NQS [27][28][29][30].

Considering the relation between neural networks (NN) and TN, the first works focused on the restricted Boltzmann machines (RBM), which are one of the simplest classes of NN. It is impossible to efficiently map an RBM onto a TN, as they correspond to string-bond states with an arbitrary nonlocal geometry [28]. This result was later refined to show that an RBM may correspond to an MPS with an exponentially large bond dimension, and only short-range RBM can be mapped onto efficiently computable entangled plaquette states [31]. Similar results have been obtained that deep Boltzmann machines with proper constraints can be mapped onto TN that are efficiently computable through transfer matrix methods [32].

Among other classes of NN, it has been shown that simple cases of convolutional neural networks (CNN) and RNN, which are also instances of arithmetic circuits or sum-product networks [START_REF] Poon | IEEE International Conference on Computer Vision Workshops[END_REF], can be efficiently mapped onto TN. Such efficient mapping is no longer possible with modifications to those simple cases, including overlapped convolution kernels, stacked recurrent layers [START_REF] Levine | [END_REF], or nonlocal score function [START_REF] Khrulkov | 6th International Conference on Learning Representations[END_REF], which suggests that the reuse of information in general NN architectures is an extensive source of expressivity. In the opposite direction, a mapping from arbitrary TN to feedforward NN can be constructed using NN layers that contract a set of tensor indices at a time [START_REF] Sharir | [END_REF].

In this work, we propose a new family of variational ansätze that shares characteristics with NN and TN. We start by establishing an exact mapping from MPS to an RNN with a linear update rule for the memory and a quadratic output layer, which we refer to as 1D MPS-RNN. We then propose a generalization of this architecture to 2D lattices using a linear combination of memory components at previous neighboring sites in each memory update, which we call 2D MPS-RNN. After a theoretical discussion of the limitation on entanglement entropy due to the linear memory update, we propose a multilinear update with inspiration from the tensorial nature of the PEPS architecture. This ansatz can represent states with area law of entanglement entropy, which is particularly interesting in quantum physics. Unlike PEPS, our ansatz is still a generative model, which supports perfect sampling and wave function evaluation in polynomial time. Finally, we conduct numerical experiments to compare

M (0) γ (0)
M (1) γ (1) M (2) γ (2) ∏ √ h (-1) h (0) h (1) h (2)

σ 0 σ 1 σ 2 p(σ 0 ) p(σ 1 |σ 0 ) p(σ 2 |σ 0 , σ 1 )
ψ h (2) = h (1) M (2) σ 2 p(σ 2 |σ 0 , σ 1 ) ∝ h (2) † γ (2) h (2) (a) (b) the performances of MPS and our ansätze on the antiferromagnetic Heisenberg model (AFHM) on square lattice and triangular lattice. The results show that our ansätze produce lower variational energy with bond dimensions smaller by orders of magnitude compared to MPS, even in the presence of frustration, and the tradeoff between computation cost and accuracy can be systematically controlled by the bond dimension, which is the only hyperparameter of our architectures.

II. RNN QUANTUM STATES

We can write any wave function ψ of a quantum spin-1/2 system consisting of V sites in the form

ψ(σ) = i p(σ i |σ <i ) e iφ(σ) , (1)
where σ = (σ 0 , . . . , σ V -1 ) is the spin configuration in the z-basis, and p(σ i |σ <i ) is the conditional probability of measuring the spin i given the measurements of previous spins σ <i = (σ 0 , . . . , σ i-1 ). An RNN quantum state [START_REF] Hibat-Allah | [END_REF] uses a memory vector h (i) at each step to summarize the information about the previous spins σ ≤i , which is updated by a function

h (i) = f (i) update (σ i , h (i-1)
). Using this information, we compute each conditional probability from an output function p(σ i |σ <i ) = f (i) out (h (i) ). Therefore, the RNN quantum state is defined by the update function f (i) update , the output function f (i) out , and the phase function φ, as illustrated in Fig. 1 (a). The purpose of introducing h (i) instead of working directly with the spin values is the compression of the exponential amount of information in the wave function while keeping a rich expressivity. The trade-off between tractability and expressivity is usually controlled by the complexity of f (i) update and f (i) out , and the size of h (i) . A notable feature of RNN quantum states is the perfect sampling from the probability distribution of the spins, without the need of Markov chains and therefore avoiding the problem of autocorrelation, which is particularly advantageous when stochastically optimizing the wave function by VMC.

III. MAPPING MPS TO RNN

In the following, we present an exact mapping from an MPS of bond dimension χ to a specific RNN architecture of memory dimension χ. Accordingly, in this manuscript, we use the term "bond dimension" when referring to the dimension of the memory of our RNN architectures to underline the fact that they are conceptually similar. It is perhaps not completely surprising that a mapping from MPS to RNN exists as MPS has been shown to allow perfect sampling [38] and as it achieves an efficient compression of the information of the wave function for many quantum systems.

The MPS ansatz for a wave function ψ is defined by

ψ(σ) = χ-1 s0,...,s V =0 V -1 i=0 M (i) σi;si+1,si , (2)
where

M (i)
σi is a complex χ×χ matrix that depends on the spin σ i at the site i. Note that we keep the indices s 0 and s V at both ends to simplify the discussion. To rewrite the MPS as an RNN, we identify the intermediate result of the tensor contraction as the memory h (i) , which satisfies a local update rule

h (i) = M (i) σi h (i-1) , (3)
where h (i) ∈ C χ is a vector for each site i, and implicitly depends on the previous spins σ ≤i = (σ 0 , . . . , σ i ).

The boundary condition is h (-1) = (1, . . . , 1). This is useful as the conditional probability of the MPS is then proportional to a positive semi-definite quadratic form of the memory

p(σ i |σ <i ) ∝ [h (i) ] † γ (i) h (i) , (4)
where the explicit form of γ (i) as a contraction of the

M (i)
σ matrices is given for completeness in the Supplemental Material S1. The phase of ψ can be obtained from the memory at the last site:

φ(σ) = arg s h (V -1) s . (5)
We now have all the elements to introduce the ansatz that we refer to as the vanilla MPS-RNN: Eq. (3) defines For each site of the square lattice (x, y), the RNN takes as input the memory vectors h (x-1,y) and h (x,y-1) from previous neighboring sites and the spin value σx,y at the current site, and outputs the updated memory vector h (x,y) to be used in future neighboring sites and the conditional probability p(σx,y|σ <(x,y) ).

∏ √ ψ h (x-1,y) h (x,y-1) h (x,y) h (x,y) σ x,y p(σ x,y |σ <(x,y) ) (a) (b)
the memory update in terms of the variational parameters

M (i)
σ , the conditional probability is obtained from Eq. ( 4), and the phase is obtained from Eq. ( 5). See Fig. 1 (b) for TN diagrammatic illustrations of the memory update and the conditional probability. We elevate each γ (i) to be a free variational parameter, independent of M (i) σ . Therefore, while any MPS can be exactly mapped to a vanilla MPS-RNN with the same bond dimension, the opposite is not true, and the freedom of γ (i) brings additional expressivity. As M (i) σ and γ (i) depend on the spatial position i, we can say that we are encoding the spatial dimension into the "time" dimension of the RNN.

We slightly generalize the vanilla MPS-RNN architecture to improve its numerical performance by adding the normalization of the memory at each step and including a vectorial term in the memory update, as described in the Supplemental Material S2. We call the resulting ansatz 1D MPS-RNN, which is the one we use in numerical experiments.

IV. 2D MPS-RNN

From the known success of MPS and its exact mapping to 1D MPS-RNN, we conclude that the latter can represent the ground states of gapped 1D systems. Therefore, the nonlinearity of conventional RNN is not necessary to efficiently represent short-range entanglement in one dimension, as 1D MPS-RNN only uses a linear memory update. By taking inspiration from the PEPS architecture, we consider a minimal generalization of the architecture to 2D systems, where we seek to efficiently approximate their ground states using linear or multilinear memory updates while keeping the computation time and memory scaling polynomially with the system size.

In the following, for simplicity, we limit our discussion to a square lattice of size V = L × L. As the RNN se- quentially outputs the conditional probability of a given spin at each step, we need to define a 1D ordering for the sites on the 2D lattice. We use the "snake" ordering as commonly used in MPS [39] and RNN [START_REF] Hibat-Allah | [END_REF] for 2D inputs, which is illustrated in Fig. 2. A minimal generalization of Eq. ( 3) to 2D consists in using a linear combination of the memory components at two previous neighboring sites, and we call the resulting architecture 2D MPS-RNN. We provide the TN diagram for the memory update of 2D MPS-RNN in Fig. 3, which is a generalization of Fig. 1 (b). The equations are detailed in the Supplemental Material S3.

h (2,2) = T (2,2) h (1,2) h (2,1) σ 2,2 + M (2,2) y h (2,1) σ 2,2 + M (2,2) x h (1,2) σ 2,2 + v (2,2) σ 2,2
A notable difference between 2D MPS-RNN and TN is that in the former each memory vector is used in two future neighboring sites. In a TN, on the contrary, each tensor only appears once in the contraction. This reuse of information can lead to a more efficient state compression using fewer parameters [START_REF] Levine | [END_REF]. We remark that while the shallow RAC TN in [40] is only a subset of MPS, our vanilla MPS-RNN is a superset of MPS.

The direct information flow between vertically neighboring sites is an apparent advantage of 2D MPS-RNN over 1D architectures. Thanks to this, the memory at a given site no longer needs to be carried along O(L) memory updates to influence the memory of a vertical neighbor of the 2D lattice. However, because of the linearity of the memory update proposed above, a 2D MPS-RNN with bond dimension χ can be exactly simulated by a 1D MPS-RNN with bond dimension L χ, as proved in the Supplemental Material S4. Therefore, 2D MPS-RNN can reduce the bond dimension χ at most linearly in L compared to 1D MPS-RNN rather than exponentially. In particular, its entanglement entropy cannot have an area law because it can scale at most logarithmically with the length L of a horizontal cut at a fixed χ.

V. TENSOR-RNN

We are therefore led to consider a nonlinear generalization of the memory update Eq. ( 3), and a minimal choice is a multilinear function of the memories at pre- vious neighboring sites. The resulting memory update is sketched in Fig. 3 with the help of a TN diagram and detailed in the Supplemental Material S3. We call this architecture tensor-RNN. It can represent states with area law of entanglement entropy, and an example of this is constructed in the Supplemental Material S5.

Compared to other methods employing 2D RNN as NQS, such as [START_REF] Hibat-Allah | [END_REF][START_REF] Hibat-Allah | [END_REF], our ansätze are developed by generalizing TN and therefore retain some advantages of TN, namely (1) their expressivity is controlled by a single hyperparameter, the bond dimension, which avoids the laborious architecture search of RNN; (2) their linear or multilinear architectures are simple enough to enable theoretical analysis tools inspired from TN, such as entanglement entropy; and (3) they provide an effecive way of initialization when used in VMC, which we discuss below.

The ansätze we have discussed form a hierarchy: MPS 1D MPS-RNN 2D MPS-RNN tensor-RNN, where an architecture on the right can exactly simulate an architecture on the left using the same bond dimension. We refer to this procedure as hierarchical initialization. During the variational optimization, we first use the DMRG algorithm to optimize a MPS, then use the optimized parameters to initialize a 1D MPS-RNN. After the gradientbased optimization of the 1D MPS-RNN ansatz, we use it to initialize a 2D MPS-RNN, which after optimization can be used to initialize a tensor-RNN. This procedure provides reasonable starting points for each optimization and avoids being stuck early in a local minimum with high energy.

VI. NUMERICAL EXPERIMENTS

To numerically evaluate the performances of our proposed ansätze, we start with the AFHM, defined by the Hamiltonian Ĥ = i,j Ŝ(i) • Ŝ(j) , where i, j denotes a pair of nearest neighbors in the lattice. We perform numerical experiments on a 10×10 square lattice with open boundary conditions (OBC). We apply the Marshall sign rule [43] to the Hamiltonian, which makes the ground-state wave function positive. We use the standard VMC method to optimize our ansätze, with automatic differentiation (AD) [44] to compute the gradients, perfect sampling [45], Adam optimizer [START_REF] Kingma | 3rd International Conference on Learning Representations[END_REF], and the hierarchical initialization described above. The variational energies are compared to the quantum Monte Carlo (QMC) result [41]. More information on the numerical implementation can be found in the Supplemental Material S7. Fig. 4 (a) shows that we can systematically obtain lower variational energies by increasing the bond dimension χ. We remark that the memory updates of 2D MPS-RNN and tensor-RNN can reduce the bond dimension by 1-2 orders of magnitude, which confirms our expectation for 2D systems. The significant compression achieved by tensor-RNN over 2D MPS-RNN shows that the tensorial part of the memory update is crucial to capturing 2D quantum entanglement, as expected from the TN analogy.

Finally, we apply our ansätze on a frustrated system, namely the AFHM on a 10 × 10 triangular lattice with OBC. The ground state wave function can no longer be made positive by the Marshall sign rule, and the hierarchical initialization is particularly useful in this case, as it alleviates the difficulty of learning the sign structure of the wave function for frustrated systems [START_REF] Westerhout | [END_REF]48]. Fig. 4 (c) shows that the reduction of bond dimension in 2D MPS-RNN and tensor-RNN is unaffected by the frustration. Moreover, a comparison of square and triangular lattices provides evidence that the multilinear memory update of tensor-RNN is important to capture the sign structure.

VII. CONCLUSIONS

We have shown that an RNN can exactly represent an MPS of bond dimension χ with a linear memory update, which we call 1D MPS-RNN. Using analogies with TN and analytical arguments, we have proposed two minimal generalizations of the RNN architecture better suited to capturing 2D quantum entanglement while still using linear or multilinear memory updates. The tensorial version, tensor-RNN, can represent states with the area law of entanglement entropy while keeping the convenient properties of efficient, perfect sampling and wave function evaluation. We have provided numerical evidence for the square AFHM that our ansätze can achieve performance comparable to an MPS of bond dimension higher by orders of magnitude. These results have also been confirmed in the presence of frustration in the case of the triangular AFHM, where we have found that the hierarchical initialization from MPS to tensor-RNN is particularly useful to learn the sign structure of the wave function, which is generally a challenge for NQS.

Several possible future directions can be envisaged. In this work, we have not used one of the essential features of neural networks: nonlinear activation. Future work is needed to understand to what extent nonlinearities help increase the expressivity of representing quantum states and how they affect the optimization of tensorial RNNs. It has also been shown that deep architectures in modern neural networks can more efficiently produce entanglement [START_REF] Sharir | [END_REF], and our tensorial architectures can be generalized in a multi-layer fashion. Another promising research direction for the architectures concerns physical symmetries, and we may incorporate the quantum number conservation and the SU(2) symmetry of MPS [49] into our ansätze. With those techniques, we expect tensorial RNNs to find broader applications.

Our code is available at https://github.com/ wdphy16/mps-rnn Supplemental Material for "From Tensor Network Quantum States to Tensorial Recurrent Neural Networks"

S1. CONDITIONAL PROBABILITY OF MPS

M (0) M (1) M (2) M (3) M (4) M (5) M (0) † M (1) † M (2) † M (3) † M (4) † M (5) † s 3 t 3 ↑ ↓ σ 2 h (2)
γ (2) Figure S1. Tensor contraction of a MPS with 6 sites to obtain the unnormalized conditional probability p(σ2|σ0 = ↑, σ1 = ↓).

The construction of conditional probabilities in the vanilla MPS-RNN comes from the fact that in MPS we can efficiently trace out the spins σ >i and specify the spins σ <i to compute the unnormalized conditional probability p(σ i |σ <i ), as shown in Fig. S1, and the normalization is given by σ∈{↑,↓} p(σ|σ <i ) = 1.

(S1)

When mapping an MPS to a vanilla MPS-RNN, each coefficient matrix γ (i) can be computed from the traced out matrices {M (j)

σ |j > i}: γ (i-1)
ti,si = si+1,...,s V ti+1,...,t V σi,...,σ V -1

V -1 j=i [M (j) σj ;tj+1,tj ] * M (j) σj ;sj+1,sj . (S2)

S2. 1D MPS-RNN

We append a few generalizations to the vanilla MPS-RNN defined by Eqs. ( 3), (4), and (5) to ease its practical use in numerical experiments with gradient-based optimization. As each γ (i) matrix is positive semidefinite, we diagonalize it by γ (i) = [U (i) ] † diag(η (i) ) U (i) , where U (i) is a unitary matrix and η (i) is a real and positive vector. We further absorb U (i) by h (i) ← U (i+1) h (i) and M

(i) σ ← U (i+1) M (i) σ [U (i)
] † , so that Eq. ( 4) is simplified without loss of generality to

p(σ i |σ <i ) ∝ s η (i) s |h (i) s | 2 . (S3)
During the variational optimization, we regard the parameters η (i) as independent of

M (i)
σ , and we use log η

(i)
s as free parameters to ensure the positivity of η (i) s . These free parameters give our ansatz richer expressivity than the original MPS.

Based on Eq. ( 3), we add a vector term v (i) σ ∈ C χ , which plays the role of a bias in linear layers of NN, and normalize h (i) at each step to improve numerical stability:

h(i) (σ i , σ <i ) = M (i) σi h (i-1) (σ i-1 , σ <i-1 ) + v (i) σi , (S4)
h (i) s (σ i , σ <i ) = h(i) s (σ i , σ <i ) t,σ h(i) t (σ, σ <i ) 2 . ( S5
)
Note that here we explicitly show the dependence of h (i) on the spins σ i and σ <i . When training the ansatz, not including the vector term v

(i)
σ can cause gradient explosion because there is no bounded nonlinear activation between RNN steps as in conventional RNN architectures to suppress high-order powers of parameters, and the v

(i)
σ term helps regularize the gradients. For a trained ansatz, we can always absorb v

(i) σ into an extra row of M (i)
σ , and the normalization does not affect the conditional probability. Eq. (S4) is also a special case of Eq. (S8) for tensor-RNN with T (x,y) σ = 0, M (x,y) y;σ = 0. We also obtain the phase from a generalization of Eq. ( 5):

φ(σ) = i φ (i) , φ (i) = arg [w (i) σi ] T h (i) + c (i) σi , (S6)
where w

(i) σ ∈ C χ and c (i)
σ ∈ C are new variational parameters. There is a local phase at each RNN step, which directly contributes to the global phase and receives gradient information from the loss. This representation of the phase is particularly useful when optimizing on frustrated systems, although the physical correspondence of the local phases remains an open question. We call the above architecture 1D MPS-RNN, which is defined by Eqs. (S4), (S3), and (S6).

In the simplest form of commonly used RNN architectures in machine learning, the recurrence equation is defined by the following instead of Eq. (S4):

h (i) (σ i , σ <i ) = f (M h (i-1) + σ i w + v), ( S7
)
where M is a matrix of parameters, w and v are vectors of parameters, and they are independent of the site index i or the input σ i . The function f is a nonlinear activation function, usually the sigmoid function. More advanced recurrence equations, such as long short-term memory (LSTM) [50] and gated recurrent unit (GRU) [START_REF] Cho | Proceedings of SSST@EMNLP 2014, Eighth Workshop on Syntax, Semantics and Structure in Statistical Translation[END_REF], have been proposed to better represent long-term correlations. In comparison, we use site-and input-dependent parameters such as

M (i)
σ and h

(i)
σ to create the exact mapping between MPS and RNN, and we do not use the nonlinear activation.

S3. MEMORY UPDATE FOR 2D MPS-RNN AND TENSOR-RNN

As illustrated in Fig. 3 σx,y is a complex tensor of dimensions χ × χ × χ for each site (x, y) and depends on the spin σ x,y . In the matrix terms, the subscripts x and y (in upright font) denote two matrices at the same site and acting on different input memory vectors. The sign ± takeswhen y is even and + when y is odd, as defined by the snake ordering and shown in Fig. 2. The boundary conditions are h (-1,0) = 1 and h (-1,y) = h (L,y) = h (x,-1) = 0 for all x and for y > 0. We call this architecture tensor-RNN, which is defined by Eqs. (S8), (S3), and (S6). The 2D MPS-RNN is a special case of tensor-RNN with T (x,y) σ = 0. To discuss the computational complexity of tensor-RNN, we focus on one RNN step and use simplified notations in Eq. (S8):

h(1) s = t,u T s,t,u h (2) t h (3) u + t M (x) s,t h (2) t + t M (y) s,t h (3) t + v s . (S9)
Evaluating the memory update takes O(χ 3 ) time in each RNN step, which is dominated by the tensor term t,u T s,t,u h

(2)

t h (3)
u for all s, and the evaluations of Eqs. (S3) and (S6) take a lower order of time compared to Eq. (S8). It takes O(V χ 3 ) time when running through V RNN steps to evaluate the wave function, where V = L 2 is the number of lattice sites. This is exponentially lower than the time to exactly contract a PEPS, and polynomially lower than typical methods of approximate contraction with SVD truncation. For example, the method in [42] uses O(χ 6 ) time. Without the tensor term, the time complexity becomes O(V χ 2 ) for 1D and 2D MPS-RNN.

The static memory (used to store the parameters) is O(V χ 2 ) for 1D and 2D MPS-RNN, and O(V χ 3 ) for tensor-RNN. The dynamic memory (additionally used when evaluating the wave function) is O(χ) for 1D MPS-RNN as we only need to store the memory of the current site, and O(Lχ) for 2D MPS-RNN and tensor-RNN as we also need to store the memories in the last row. This is also exponentially lower than exactly contracting a PEPS, and in practice we have found that the bottleneck of scaling up our ansätze is usually time rather than memory.

When computing the gradients of the variational energy w.r.t. the parameters of tensor-RNN, following the derivative rule of tensor products, we have

∂ h(1) s ∂T s ,t,u = δ s,s h (2) t h (3) u , ∂ h(1) s ∂M (x) s ,t = δ s,s h (2) t , ∂ h(1) s ∂M (y) s ,t = δ s,s h (3) t , ∂ h(1) s ∂v s = δ s,s . (S10)
Therefore, computing the gradients also takes O(χ 3 ) time in each RNN step, which is dominated by the tensor term

∂ h(1) s ∂T s ,t,u
for all s, t, u, and it becomes O(V χ 3 ) when running through V RNN steps. As the forward evaluation of the wave function also takes O(V χ 3 ) time, we conclude that the time complexity of an optimization step is O(V χ 3 ). Without the tensor term, it becomes O(V χ 2 ) for 1D and 2D MPS-RNN. The dynamic memory for an optimization step is O(V χ) for all our ansätze, as we need to store the memories at all sites to compute the gradients using backpropagation.

S4. 2D MPS-RNN SIMULATED BY 1D MPS-RNN

In the following, we show that a 2D MPS-RNN with bond dimension χ can be exactly simulated by a 1D MPS-RNN with bond dimension Lχ. The 1D MPS-RNN uses the same snake ordering on the 2D lattice, but there is only one matrix term at each site. The idea is that the memory components in each row of the 2D MPS-RNN are linear combinations of those in the previous row.

Without loss of generality, we assume y is even and y > 0. Before updating along the y-th row in the 1D MPS-RNN, the memory h (0,y-1) is a concatenation of the memories h (i,y-1) of the 2D MPS-RNN in the previous row:

h (0,y-1) iχ+s = h (i,y-1) s (S11)
for i ∈ {0, . . . , L -1} and s ∈ {0, . . . , χ -1}. As we update along the y-th row, the components in h (x,y) are gradually replaced by the memories of the 2D MPS-RNN in the current row:

h (x,y) iχ+s = h (i,y) s if i ≤ x, h (i,y-1) s if i > x.
(S12)

The corresponding transfer matrix is defined by

M (x,y) σ;iχ+s,jχ+t =          M (x,y) y;σ;s,t if i = j = x, M (x,y) x;σ;s,t if i = x -1 and j = x, δ s,t if i = j = x, 0 otherwise. (S13)

S5. TENSOR-RNN CAN REPRESENT A QUANTUM STATE WITH AN AREA LAW OF ENTANGLEMENT ENTROPY

In the following, we construct a tensor-RNN that has an exact area law of entanglement entropy for any region that follows the snake ordering and has at least one row. The idea is that we can always have L independent memory vectors on the boundary of the region, and the wave function Eq. (S23) is a nontrivial multilinear function of them.

We start by introducing a new type of variational ansatz, which we call late normalizing (LN) tensor-RNN. We use the same update rule as Eq. (S8) with the snake ordering and an even L:

h (x,y) s = t,u T (x,y) σx,y;s,t,u h (x±1,y) t h (x,y-1) u + t M (x,y) x;σx,y;s,t h (x±1,y) t + t M (x,y) y;σx,y;s,t h (x,y-1) t + v (x,y)
σx,y;s , (S14)

h (0,2) h (1,2) h (1,2) h (2,1) h (3,1)
Figure S2. Sketch of the LN tensor-RNN architecture. Red sites are in a region following the snake ordering, for which we compute the reduced density matrix ρ ≤ (1,2) . The dark red sites are on the boundary B(1, 2), and the red edges are the memories h (x,y) on the boundary. Green sites are in the traced-out region.

except that we do not apply the normalization Eq. (S5) or compute the conditional probability Eq. (S3) at each step, but only take the 0-th component of the memory vector at the last step to be the wave function:

ψLN (σ) = h (0,L-1) 0 , (S15)
and normalize it afterwards:

ψ LN (σ) = ψLN (σ) σ ψLN (σ ) 2 . (S16)
We remark that a general LN tensor-RNN cannot be represented by a tensor-RNN or a TN. Despite that, it is useful for our theoretical analysis, as it allows us to obtain analytical expressions for the entanglement entropy without worrying about the normalization. In the following, we construct a particular LN tensor-RNN such that ψLN (σ) is always a multilinear function of the memories on the boundary, as illustrated in Fig. S2.

We consider a LN tensor-RNN with χ = 2. In the last row y = L -1 of the lattice, the parameters are defined as

T (x,L-1) σ;s,t,u = δ s,0 δ t,0 δ u,σ if x < L -1, 0 if x = L -1, M (x,L-1) x;σ = 0, M (x,L-1) y;σ;s,t = 0 if x < L -1, δ s,0 δ t,σ if x = L -1, v (x,L-1) σ = 0, (
S17) where in δ u,σ we represent σ by {0, 1}. We can check that ψLN (σ) is a multilinear function of the memories on the boundary:

ψLN (σ) = h (0,L-1) 0 = h (x+1,L-1) 0 x x =0 h (x ,L-2) σ x ,L-1 for any x < L -1 (S18) = L-1 x =0 h (x ,L-2) σ x ,L-1 . (S19)
In other rows y < L -1, the parameters are defined as

T (x,y) σ = 0, M (x,y) x;σ = 0, M (x,y) y;σ;s,t = 0 if y = 0, δ s,t if y > 0, v (x,y) σ;s = δ s,σ if y = 0, 0 if y > 0. (S20)
For any x and 0 < y < L -1, assuming an even y without loss of generality, we can check that ψLN (σ) =

x-1

x =0 h (x ,y) σ x ,L-1 L-1 x =x h (x ,y-1) σ x ,L-1 , (S21) h (x,y) s = δ s,σx,0 . (S22)
Therefore, we always have L memory vectors on the boundary, each has length χ = 2 and only depends on the spin at the first row and the same column σ x,0 . We note that although some h (x,y) can appear twice on the boundary, like h (1,2) in Fig. S2, they still appear only once in ψLN (σ) because of the particular choice of parameters. Eventually, ψLN (σ) can be simplified to

ψLN (σ) = L-1 x=0 δ σx,0,σ x,L-1 . (S23)
We denote by ρ ≤(x,y) the reduced density matrix of the system obtained by tracing out σ >(x,y) , where σ >(x,y) is the set of all spins σ x ,y whose coordinates (x , y ) come after (x, y) in the snake ordering of the lattice. Using the wave function in Eq. (S23), the unnormalized reduced density matrix is ρ≤(x,y) σ ≤(x,y) , σ ≤(x,y) =

σ >(x,y) ψLN σ ≤(x,y) , σ >(x,y) ψ * LN σ ≤(x,y) , σ >(x,y) = L-1 x=0 δ σx,0,σ x,0 (S24)
for 0 < y < L -1, which shows the area law of the entanglement entropy

S = L log 2. (S25)
Finally, we show that this LN tensor-RNN can be represented by a tensor-RNN. As long as ψLN (σ) is a multilinear function of the memories on the boundary, ρ ≤(x,y) can be written as a multi-quadratic function of them, analogously to Eq. ( 4):

ρ ≤(x,y) σ ≤(x,y) , σ ≤(x,y) = s B(x,
H (x,y) s B(x,y) σ ≤(x,y) = (x ,y )∈B(x,y) h (x ,y ) s x ,y σ ≤(x,y) (S27)
is the direct product of all memory vectors on the boundary, and γ >(x,y) is a tensor. If we can also write γ >(x,y) as a direct product over sites: The tensor term in Eq. (S8) introduces an O(χ 3 )-growing time and memory consumption when evaluating the tensor-RNN, which becomes a challenge when χ gets large. To avoid that, we apply a Tucker decomposition [START_REF] Tucker | [END_REF] of the tensor:

T (x,y) σ;s,t,u = s ,t ,u K (x,y) σ;s ,t ,u U (x,y) o;σ;s,s U (x,y) x;σ;t,t U (x,y) y;σ;u,u , (S31)
where K (x,y) σ is a smaller tensor of dimensions χ × χ × χ , and we heuristically choose χ = χ 2/3 (rounded up to an integer) so that the additional time and memory consumption introduced by the tensor term is limited to O(χ 2 ). During the variational optimization, K x;σ , and U (x,y) y;σ are all regarded free parameters. We call this architecture compressed tensor-RNN. In the numerical results, we have found that compressed tensor-RNN intuitively produces a variational energy between 2D MPS-RNN and tensor-RNN, whether with a fixed χ or a fixed number of parameters.

S7. DETAILS OF NUMERICAL EXPERIMENTS

The RNN architectures are implemented using NetKet [53,54]. In the numerical experiments on the square AFHM, as we know the grounds states are real and positive, we discard the parameters w (i) σ and c (i) σ and fix the phase to 0. For the triangular AFHM, the ground state is no longer positive and we indeed need those parameters for the phase. The triangular lattice is defined by adding a diagonal edge to each plaquette in the square lattice.

During the hierarchical initialization, for the M (x,y) x;σ and M (x,y) y;σ in the 2D MPS-RNN that are undefined in the corresponding 1D MPS-RNN, and T (x,y) σ in the tensor-RNN, we initialize them using a element-wise Gaussian distribution with standard deviation 10 -7 to avoid degeneracy of features.

During the training, we use Adam optimizer with batch size 1, 024 for all reported results, and we have found that the stochastic reconfiguration (SR) and the plain stochastic gradient descent (SGD) generally produce higher variational energies within the same number of training steps or wall-clock time. Gradient clipping with global norm 1 is applied to alleviate the energy spikes. After the training, each reported variational energy is estimated from 10 6 samples, and its statistical error is too small to be visible in the plots.

The results are particularly sensitive to the learning rate, because of the instability discussed in Section S2. For bond dimensions χ ≤ 10, we use learning rate 10 -2 for 10, 000 steps, 10 -3 for another 10, 000 steps, and finally 10 -4 for 20, 000 steps. We manually reduce the learning rate by a factor of 10 if we find the variational energy diverges in the middle of training. For bond dimensions χ > 10, we start from learning rate 10 -3 . The training of a tensor-RNN with χ = 64 on the 10 × 10 triangular AFHM takes about 2 days on an Nvidia V100 GPU.

We have computed all the DMRG results in Fig. 4, except for χ ≥ 1, 024 on the 20 × 20 lattice. The DMRG algorithm is implemented using ITensors.jl [55], with the snake ordering, 20 sweeps, eigenvalue cut-off threshold 10 -12 , and noises decaying exponentially from 10 -3 to 10 -12 in the first 10 sweeps. 4 as a function of the number of parameters. Our ansätze only need 2 to 3 orders of magnitude fewer parameters than MPS to obtain the same energy error. At a fixed number of parameters, there is a clear difference in the energy error between the 1D and the 2D ansätze, which shows that the 2D ansätze can more efficiently compress the approximated ground state. σx,y h (x±1,y) ⊗ h (x,y-1) (areas of red triangles), the matrix terms M (x,y)

S8. ENERGY ERROR AS FUNCTION OF THE NUMBER OF PARAMETERS

S9. CORRELATIONS OF SPINS

x;σx,y h (x±1,y) and M (x,y) y;σx,y h (x,y-1) (areas of rectangles), and the vector term v (x,y) σx,y (areas of blue circles) in Eq. (S8) at each site (x, y), in our ansätze with χ = 10 trained on the 10 × 10 square AFHM. The green rectangles are the matrix terms along the snake ordering, and the cyan ones are not. The magnitude of each term a is defined by the norm s η (x,y) s |as| 2 , then averaged over the configurations σ ∼ |ψ(σ)| 2 . At each site, we normalize the four magnitudes to see their relative contributions to the memory at that site. In 1D MPS-RNN, the matrix term becomes larger near the ends of each row, because it needs to carry all information in the current row to the next row. In 2D MPS-RNN that is not the case, thanks to the vertical connections. In tensor-RNN, the tensor terms in the last row are particularly large, because they are essential for expressing the wave function as a multilinear function of the memories and representing the area law of entanglement entropy, as shown by the example in Section S5.

Figure 1 .

 1 Figure 1. (a) Computational graph of the vanilla MPS-RNN quantum state. The phase factor φ(σ) is omitted in the figure. (b) TN diagrams of the memory h (i) and the conditional probability p(σi|σ<i) for the vanilla MPS-RNN.

Figure 2 .

 2 Figure 2. (a) Computational graph of 2D MPS-RNN and tensor-RNN architectures for square lattice with snake ordering. (b) Inputs and outputs at each site of 2D MPS-RNN and tensor-RNN architectures.For each site of the square lattice (x, y), the RNN takes as input the memory vectors h (x-1,y) and h (x,y-1) from previous neighboring sites and the spin value σx,y at the current site, and outputs the updated memory vector h (x,y) to be used in future neighboring sites and the conditional probability p(σx,y|σ <(x,y) ).

Figure 3 .

 3 Figure 3. TN diagram of the memory update for tensor-RNN. The memory update for 2D MPS-RNN is obtained for T (2,2) = 0, and the memory update for 1D MPS-RNN can be derived by further setting M (2,2) y = 0.

Figure 4 .

 4 Figure 4. Energy relative error of MPS (optimized by DMRG) and our ansätze (optimized by VMC) as a function of the bond dimension χ. The Hamiltonian is AFHM on (a) 10 × 10 square lattice and (b) 10 × 10 triangular lattice, and the energy is compared to (a) QMC [41] and (b) DMRG with χ = 4096. The PEPS result in (a) is taken from [42] with χ = 10.

γ>

  (x,y) s B(x,y) ,s B(x,y) = (x ,y )∈B(x,y) γ (x ,y ) s x ,y ,s x ,y , (S28) and it does not depend on the spins σ ≤(x,y) , σ ≤(x,y) , then we can represent the LN tensor-RNN by a tensor-RNN. With our particular choice of parameters, we simply have γ (x,y) sx,y,s x,y = δ sx,y,s x,y . Along with the above definitions of T (x,y) σ , M (x,y) x;σ , M (x,y) y;σ , v (x,y) σ, we define the remaining parameters of the tensor-RNN asη (x,y) = 1 for y < L -1, η (x,L-1) s = δ s,0 ,(S29)w (x,y) σ;s = δ s,0 if x = 0 and y = L -1, 0 otherwise, c (x,y) σ = 0 if x = 0 and y = L -

  (x,y) σ , U (x,y) o;σ , U (x,y)

Figure S3 .

 S3 Figure S3. Energy relative error in Fig.4as a function of the number of parameters. Our ansätze only need 2 to 3 orders of magnitude fewer parameters than MPS to obtain the same energy error. At a fixed number of parameters, there is a clear difference in the energy error between the 1D and the 2D ansätze, which shows that the 2D ansätze can more efficiently compress the approximated ground state.

C

  Figure S4. (a) Qualitative evaluation of the correlations of all spins with the spin σ5,5 (left column) and σ1,5 (right column), in a 1D MPS-RNN (top row) and a 2D MPS-RNN (bottom row) with χ = 10 trained on the 10 × 10 square AFHM. The 2D MPS-RNN produces more isotropic correlations, especially near the ends of rows, thanks to the vertical connections. (b) Quantitative evaluation of the correlations C of the spins σ d,10 with σ0,10 (horizontal) and the spins σ 10,d with σ10,0 (vertical), in a 1D MPS-RNN and a 2D MPS-RNN with χ = 16 trained on the 20 × 20 transverse field Ising model at the critical point. The 2D MPS-RNN captures the algebraic decay of the correlations in both directions, while the 1D one fails to do so.

  Figure S6. Energy in a typical optimization procedure with hierarchical initialization (HI) through 1D MPS-RNN, 2D MPS-RNN, compressed tensor-RNN (see Sec. S6), and tensor-RNN, compared to a tensor-RNN with random initialization instead of HI. The Hamiltonian is the 10 × 10 triangular AFHM, and all ansätze have bond dimension χ = 16. The HI procedure starts at the energy of DMRG, and we can clearly see that each stage of HI lowers the energy as new variational parameters are introduced, while the wave function without HI gets stuck at a high energy.

  y) ,s B(x,y) is the boundary of the region, s B(x,y) = {s x ,y |(x , y ) ∈ B(x, y)} are the memory indices on the boundary,

	H (x,y) s B(x,y) σ ≤(x,y) γ	>(x,y) s B(x,y) ,s B(x,y)	σ ≤(x,y) , σ ≤(x,y) H s B(x,y) (x,y)	σ ≤(x,y)	*	,
					(S26)
	where B(x, y)

ACKNOWLEDGMENTS

We thank Jannes Nys and Or Sharir for their valuable remarks. Support from the Swiss National Science Foundation is acknowledged under Grant No. 200021 200336.