domain: spi phys

text:
J Lefieux
email: julien.lefieux@onera.fr

E Garnier
email: eric.garnier@onera.fr

N D Sandham
email: n.sandham@soton.ac.uk

DNS Study of Roughness-Induced Transition at Mach 6

In hypersonic flows, it is useful to be able to trip the transition to turbulence upstream of air intakes for example. In the present work, direct numerical simulations have been performed of the flow past an isolated roughness element at Mach 6. First, the capability of two solvers to compute laminar and transitional flow involving freestream disturbances was demonstrated. A series of simulations was then carried out without acoustic perturbation of the freestream. The Reynolds number was increased from 14,000, to 28,000 and then to 40,000. The first two cases remain laminar within the computational domain, whereas the last case undergoes a selfsustained transition to turbulence. A response to a perturbation impulse shows the presence of a varicose mode at the intermediate Reynolds number and a sinuous mode at the largest Reynolds number.

I. Nomenclature

II. Introduction

I Mach 6 flows, the weak growth factor of natural instabilities makes it difficult for the flow to go through transition to turbulence. Shock-wave interactions with boundary layers are likely to create separation upstream of engine intakes, thus decreasing the flight performance of future hypersonic vehicles. Prediction of transition in hypersonic flows is also very important because of the impact it has on heat and drag [START_REF] Wang | Response of a hypersonic boundary layer to wall blowing-suction[END_REF]. Roughness elements may be used to act as a source of disturbance for a developing laminar boundary layer to trigger the transition to turbulence upstream of the engine intakes.

An analogy to low-speed flows can be made with respect to the following mechanisms. Due to the recirculation region upstream of the roughness element, a steady horseshoe vortex which is similar to the laminar horseshoe vortex system from Baker [START_REF] Baker | The laminar horseshoe vortex[END_REF] is wrapped around the roughness [START_REF] Bartkowicz | Numerical Simulations of Roughness Induced Instability in the Purdue Mach 6 Wind Tunnel[END_REF] [START_REF] Iyer | Roughness-induced transition in high speed flows[END_REF]. The streamwise vortices and the associated low-momentum streak can lead to a convective Kelvin-Helmholtz instability in the wake, characterized by periodic shedding of hairpin-shaped vortical structures [START_REF] Acarlar | A study of hairpin vortices in a laminar boundary layer. Part 1. Hairpin vortices generated by a hemisphere protuberance[END_REF], having a dominant role in transition. In hypersonic flow conditions, early experiments focused on the effect of roughness elements on a cone [START_REF] Van Driest | Boundary layer transition at supersonic speeds -Three dimensional rougness effects (spheres)[END_REF]. The importance of using physics based models such as the e N based on linear stability theory over empirical considerations [START_REF] Reshotko | Is Re θ /M e a meaningful transition criterion?[END_REF] had already been given in [START_REF] Morkovin | Transiton in open flow systems: a reassessment[END_REF]. The criteria commonly used to predict the transition, such as the Van Driest and Blumer correlation for cones and hemispheres, Re k , or the Shuttle Orbiter correlations, didn't take into account the three-dimensional shape of the roughness. An alternative transition prediction, based on momentum deficit due to the roughness element and taking into account the shape of the roughness was mentioned by Bernardini et al. [START_REF] Bernardini | Parameterization of Boundary-Layer Transition Induced by Isolated Roughness Elements[END_REF].

Advances in numerical methods and computing capabilities now make it possible to study hypersonic transitional flows with a good accuracy. The factors affecting whether or not a roughness element will be effective in tripping a high-speed flow to turbulence, include the Reynolds number [START_REF] Subbareddy | Direct numerical simulation of high-speed transition due to an isolated roughness element[END_REF], the roughness size, and its shape [START_REF] Van Den Eynde | Numerical Simulations of Transition due to Isolated Roughness Elements at Mach 6[END_REF], but also the wall temperature [START_REF] Redford | Compressibility Effects on Boundary-Layer Transition Induced by an Isolated Roughness Element[END_REF] and Mach number [START_REF] Iyer | Roughness-induced transition in high speed flows[END_REF]. The way in which roughness induces transition is now better understood. As described in [START_REF] Schneider | Effects of Roughness on Hypersonic Boundary-Layer Transition[END_REF], a wake with counter-rotating streamwise vortices and a possibly unstable shear layer are generated. The likelihood that transition in the wake of an isolated roughness element will occur depends on "... a race between unsteady fluctuation growth and the rapid relaxation of the basic state toward a spanwise-uniform Blasius flow". The streamwise vorticity behind small roughness elements is likely to grow via instability mechanisms such as stationary crossflow [START_REF] Saric | Stability and transtion of three-dimensional boundary layers[END_REF], Görtler vortices [START_REF] Saric | Görtler vortices[END_REF] or transient growth [START_REF] Reshotko | Role of Transient Growth in Roughness-Induced Transition[END_REF]. Transition starting from a receptivity process is also possible, due to an interaction between the roughness and the acoustic waves or other freestream disturbances leading to instability waves.

Global stability analyses have been carried out for subsonic flows but remain relatively scarce for hypersonic flows. Subbareddy et al. [START_REF] Subbareddy | Direct numerical simulation of high-speed transition due to an isolated roughness element[END_REF] performed a dynamic mode decomposition and mentioned a self-sustaining mechanism that causes the flow to transition at Mach 6 and a sufficiently large Reynolds number. Two possible feedback mechanisms were found: one due to the unsteadiness of the strong bow shock inducing pressure variations and providing energy for the periodic vortex motion; another one caused by oscillations in the near wake that trigger motion of the stagnation point which thus affects the bow shock. A global instability whose symmetric core is located in the shear layer separating the outer flow region and the boundary layer immediately downstream of the roughness element was observed at low speeds [START_REF] Citro | Global stability and sensitivity analysis of boundary-layer flows past a hemispherical roughness element[END_REF]. A global stability and sensitivity analysis of a three-dimensional Falkner-Skan-Cooke (FSC) boundary layer with a cylindrical surface roughness was performed in [START_REF] Brynjell-Rahkola | Stability and sensitivity of a cross-flowdominated Falkner-Skan-Cooke boundary layer with discrete surface roughness[END_REF] where unstable global modes responsible for transition were observed. A global stability analysis made in [START_REF] Loiseau | Investigation of the roughness-induced transition: global stability analyses and direct numerical simulations[END_REF] showed that the global instability of the near-wake region and the associated recirculation region were responsible for the sinuous instability, whereas the whole three-dimensional shear-layer surrounding the central low-speed streak induced the varicose instability.

The aim of the present work is first to validate elsA and OpenSBLI solvers on a laminar and a transitional cases. Then a numerical global stability study based on DNS at Mach 6 will be presented.

III. Numerical methods

A. Governing equations

In this work, the full compressible Navier-Stokes equations are solved for a Newtonian fluid with viscosity µ.

∂ ρ ∂t + ∂ ρu j ∂ x j = 0 ∂ ρu i ∂t + ∂ ρu i u j ∂ x j + ∂p ∂ x i = ∂τ i j ∂ x j ∂ ρE ∂t + ∂(ρE + p)u i ∂ x i = - ∂q i ∂ x i + ∂u i τ i j ∂ x j ,
where ρ is the density, u i is the velocity component in the i-coordinate direction (i=1,2,3); E = p/ρ/(γ -1) + u i u i /2 is the total energy per unit mass. The symmetric viscous stress tensor τ i j is defined as

τ i j = µ Re ∂u j ∂ x i + ∂u i ∂ x j - 2 3
∂u k ∂ x k δ i j ,
where δ i j is the Kronecker delta function defined as δ i j = 1 for i = j and δ i j = 0 for i j. The equation of state is considered for the calculations of the fluid properties. Fourier's law of heat conduction is used for the heat flux vector q j .

q j = - µ (γ -1)M 2 ∞ Pr Re ∂T ∂ x j .
The molecular viscosity is computed using Sutherland's law 1 gathers the numerical parameters which have been used in the validation of elsA and OpenSBLI, as well as those used in the reference test case with SBLI (Legacy) [START_REF] Van Den Eynde | Numerical Simulations of Transition due to Isolated Roughness Elements at Mach 6[END_REF].

µ(T) = T 3 2 1 + S * /T * r T + S * /T *
The governing equations are discretized onto a structured grid. OpenSBLI has capability for full three-dimensional curvilinear coordinates. The curvilinear system transforms the equations from a curved to an orthogonal base. A fourth-order central difference method is used to compute the spatial derivatives at the interior points of the domain. The stable boundary treatment of Carpenter et al. [START_REF] Carpenter | A Stable and Conservative Interface Treatment of Arbitrary Spatial Accuracy[END_REF] is used at the boundaries of the domain to ensure overall fourth-order accuracy. All metric terms for the coordinate transformation are also computed with a fourth order accuracy to avoid additional numerical error. The 3 r d order explicit Runge-Kutta low-storage time-stepping scheme described by Jacobs et al. [START_REF] Jacobs | OpenSBLI: A framework for the automated derivation and parallel execution of finite difference solvers on a range of computer architectures[END_REF] was used in the present simulations. High-order Weighted Essentially non-Oscillatory (WENO) and Targeted Essentially non-Oscillatory (TENO) schemes were used with OpenSBLI in the present work. A validation for the WENO-Z implementation/similarity solution is available in [START_REF] Lusher | Shock-wave/boundary-layer interactions in the automatic source-code generation framework OpenSBLI[END_REF]. Instead of selecting only one stencil over the smoothest region of the flow, WENO schemes take a weighted linear combination of all stencils. The weighting is proportional to the local smoothness over each. TENO (Targeted ENO) schemes [START_REF] Fu | A family of high-order targeted ENO schemes for compressible-fluid simulations[END_REF] have lower numerical dissipation than standard ENO and WENO methods. Instead of weighting all candidate stencils to build a flux reconstruction, they discard stencils deemed to contain discontinuities while retaining as many smooth stencils as available. In that way it benefits from improved weighting formulations and a staggered stencil arrangement to decrease the number of stencils crossed by a shock.

elsA is based on a finite volume discretization method. It has capability for implicit and explicit time integration schemes. In the first part of the validation (laminar case), the implicit Gear method of 2 nd order was used, with 20 Gear sub-iterations in order to increase the time step and reduce the computational time. In the second part of the validation (transitional case), the Runge-Kutta 4 th order explicit scheme was used to be compatible with the compact scheme. The second order AUSM+PW [START_REF] Kim | Methods for the accurate computations of hypersonic flows. I. AUSMPW+ scheme[END_REF] scheme was improved from the previous AUSM schemes by introducing pressure-based weight functions. The 6 th order compact scheme from Lele [START_REF] Lele | Compressibility effects on turbulence[END_REF] was also used in the transitional case, coupled with a 10 th order low-pass filter [START_REF] Visbal | On the Use of Higher-Order Finite-Difference Schemes on Curvilinear and Deforming Meshes[END_REF].

C. Global stability method

The 3D compressible Navier-Stokes equation may be written in a generic way as

∂q ∂t = RHS(q),
where q is the conservative variable vector ρ, ρu, ρv, ρw, ρE T . Let us consider q b as a steady flow field obtained with an implicit time-integration scheme, which will be referred to as the base flow. We use this base flow as the initial condition of the simulation. The rate of change of q to satisfy the equilibrium can then be written

∂q b ∂t = RHS(q b ).
If the flow is decomposed into its base-flow component q b and a perturbation q (q = q b + q ), then we obtain

∂q ∂t = RHS(q) -RHS(q b ).
RHS(q b ) acts like a forcing term in the governing equation. Thus, one can track the time and spatial evolution of an initial disturbance on the base flow, while the flow goes back towards the undisturbed base flow. Therefore, the most unstable (or least damped) mode can be detected. To check the validity of this method, simulations were performed with no initial disturbance and the base flow could be maintained for four flow through times. This method was found to be useful in cancelling out any discretisation differences between the implicit and explicit methods in elsA. Then, white noise disturbances were introduced in the density close to the bump and upstream of it, with a maximum amplitude four orders of magnitude smaller than the freestream.

IV. Comparison of OpenSBLI to elsA on a 3D laminar case

A. Computational set up

Computational domain & roughness element geometry

The same principle as in [START_REF] Van Den Eynde | Numerical Simulations of Transition due to Isolated Roughness Elements at Mach 6[END_REF] was used for the computational domain and for the grid generation. The numerical domain and setup can be seen on figure 1. The dimensions of the domain L x , L y and L z are respectively 150, 20 and 50 inflow displacement thicknesses. It has periodic boundary conditions at the spanwise boundaries and an isothermal no-slip condition at the wall. OpenSBLI uses Neumann boundary conditions at the outflow and top boundaries. At the inlet, the pressure is extrapolated from the interior of the domain. In elsA, the outflow and top boundaries have a supersonic outlet boundary condition. To make the compact scheme stable, a characteristic boundary condition was used for the outflow, along with a sponge layer. The inflow boundary condition in elsA is decomposed into a subsonic and a supersonic part. The subsonic part is an injection boundary condition where the total pressure and total enthalpy are given in each cell. The supersonic part imposes the conservative variables in each cell, so as to match the full similarity solution.

An isolated roughness element is placed at a distance of x r = 53 inflow displacement thicknesses from the inlet. The grid was stretched in the wall-normal direction to ensure a sufficient resolution of the boundary layer. The stretching function between the uniform computational grid (0 < η < 1) and the uniform physical grid (y 0 < y < L y ) is the following:

y = y 0 + (L y -y 0 ) sinh(b y η) sinh(b y )
, where b y = 3.19 is chosen as the stretching factor. The geometry of the bump can be described by the following equation:

y 0 ( r) = -k 2 tanh S W k tanh S H (2 r -W) + tanh S H (-2 r -W) ,
where r = √ x2 + z2 and S = cot S * . The variable S * can be seen as a smoothness factor, which controls the slope of the roughness sides and is equal to S * = 1.35 for the smooth bump case, and S * = Π 4 for the flat-top roughness element case.

Polynomial fit for boundary layer

The initialization of the flow was done based on the full compressible laminar similarity solution in [START_REF] White | Viscous Fluid Flow[END_REF]. The laminar boundary-layer profiles for u, v and T were generated numerically in Python during the code generation process in OpenSBLI. The profiles obtained at the inlet were copied over the whole domain for the initialization. This similarity profile is normalized such that the laminar displacement thickness at the domain inflow δ * 0 is unity and grows according to

δ * ( x * ) δ * 0 = ∆ √ 2Re x * Re δ * 0 , with Re x * = 1 2 Re δ * 0 ∆ 2 + Re δ * 0 x * δ * 0 ,
where x * is the dimensional streamwise coordinate in a reference frame starting from the leading edge; and ∆ is a scaling factor, depending on Mach number (Mach 6 in the current work) and wall temperature (T w = 7.02 T * r ) which is equal to ∆ = 9.071.

B. Grid convergence study

In the following simulations, a 5th order WENO-Z scheme and the AUSM+PW scheme were used in OpenSBLI and elsA respectively.

A grid convergence study was performed with both elsA and OpenSBLI. Three different grids were used and compared: Fine (900x225x350), Medium (600x150x234), and Coarse (300x100x200). The Fine mesh is close to the one used in [START_REF] Van Den Eynde | Numerical Simulations of Transition due to Isolated Roughness Elements at Mach 6[END_REF] which used a grid as fine as in the work of De Tullio and Sandham [START_REF] De Tullio | Influence of boundary-layer disturbances on the instability of a roughness wake in a high-speed boundary layer[END_REF], based on an extensive grid convergence study.

figure 2 shows a pseudo-Schlieren plot of the density gradients to better understand the topology of the flow. Note the edge of the boundary layer relative to the bump height (δ/h = 1.6). One can also notice the separation and reattachment shocks occurring near the bump. The separation/recirculation regions upstream and downstream of the bump are shown in red, using contours of small negative streamwise velocity (u = -10 -5 ). The recirculation regions are better seen

C. Results of the comparison

For elsA, the results obtained with AUSM+PW are presented here. Jameson and Roe schemes were also used and gave similar results. The Fine mesh was used for the comparison of elsA and OpenSBLI. Very good agreement was reached between them. As can be seen on figure 4, the change in velocity in all three directions is predicted the same in both codes. The only noticeable discrepancies can be found on v for instance at x = 53 or x = 58. The maximum deviation found at the reversed flow peak location (most negative value of v) at x = 53 by elsA with respect to OpenSBLI is 2.5%. ρ and T profiles also show good agreement between the two codes.

This first comparison validates the computation of viscous and metric terms in OpenSBLI, as well as the ability of 2 nd order schemes in elsA to predict hypersonic laminar flows accurately.

V. Comparison of elsA and OpenSBLI to SBLI on a transitional case A. Computational set-up

Validation of elsA and OpenSBLI with reference to the work done in [START_REF] Van Den Eynde | Numerical Simulations of Transition due to Isolated Roughness Elements at Mach 6[END_REF] with the code SBLI was carried out. The comparison is made on a transitional case, which is more discriminating for the choice of the numerical schemes used for the simulations and more sensitive to boundary conditions. The same grid as in [START_REF] Van Den Eynde | Numerical Simulations of Transition due to Isolated Roughness Elements at Mach 6[END_REF] was used for this case. The numerical domain and setup can be seen on figure 1. The grid and domain sizes are shown in Table 2. The domain was shortened in the streamwise direction from 300 to 190 inflow displacement thicknesses as the transition, which we are interested in, is supposed to occur before the end of the domain. Once again the initialization of the flow was done based on the full compressible laminar similarity solution. The Reynolds number based on the inflow displacement thickness is now Re δ * 0 = 14, 000. The freestream reference temperature is T * r = 273.15K, and the static pressure is p = 9451.319Pa. A normalization of the parameters was performed on ρ, u and T, leaving an adimensionalized static pressure of p ad = 0.01984. Figure 5 shows what the In order to keep the Reynolds number moderate and still be able to trigger instabilities, acoustic disturbances were introduced in the freestream in the work of [START_REF] Van Den Eynde | Numerical Simulations of Transition due to Isolated Roughness Elements at Mach 6[END_REF], with a 3-dimensional source placed at x f = 12.0 and y f = 6.0. The same acoustic forcing was used in the simulations with elsA and OpenSBLI solvers. This forcing consists of a source term s, added to the right hand side of the continuity equation and described by

s(x, t) = a exp(-r 2 ) M m=0 N n=1 cos(β m z + φ m ) sin(ω n t + φ n ),
where a = 3x10 -4 is the amplitude, r = (xx f ) 2 + (yy f ) 2 gives the distribution, β m = 2πm/L z is the spanwise wave number and ω n = 2πnF 0 is the frequency with F 0 = 0.02. The acoustic source has thus M = 20 spanwise modes and N = 20 temporal modes to mimic the broadband nature of acoustic fluctuations in turbulent boundary layers. The simulation integration time of statistical quantities is governed by the lower and upper bound of the frequency, while the spanwise wave number range sets the required grid resolution. φ n and φ m are random phases which prevent large spikes in the imposed forcing. These random phases used in the older simulations with SBLI were unknown meaning that different phases with elsA and OpenSBLI are likely to be present.

B. Results

Figure 6 shows the instantaneous contours of density gradient magnitude in the centreplane, computed with elsA. The contour range was adjusted to highlight the interaction between the density fluctuations caused by the acoustic forcing, the roughness element and the boundary layer, leading to transition at the end of the domain. Again the red zones show the recirculation regions which are larger than in the laminar case, due to the shape of the roughness element which is now a flat-top cylindrical roughness element. First, elsA and OpenSBLI can be compared qualitatively, by looking at figure 7 which shows instantaneous Q criterion iso-surfaces (Q = 0.003) coloured by the streamwise vorticity at different time instances. Good agreement is obtained in terms of the flow structures predicted by both solvers. They display longitudinal streaks, rolling up in the wake of the roughness element, forming hairpin-like vortices that eventually break down to turbulence.

A quantitative comparison of the three solvers is made by looking at the skin friction coefficient, averaged over one period of the forcing along the domain centreline, as shown on figure 8(a). Good agreement is obtained in the laminar part of the flow and the recirculation regions are predicted the same by both codes. The disturbance energy e is integrated over 99% of the boundary-layer thickness at the roughness centreline and shown for each code as a function of x in figure 8(b). Also the average exponential growth rate σ is computed as the average growth rate σ in the range of 100 < x < 140, and shown in the legend of figure 8(b). The receptivity process is linear, and the boundary-layer disturbances grow exponentially before reaching saturation, indicating the occurrence of nonlinear interactions. The black solid line in figure 8(b) corresponding to SBLI results [START_REF] Van Den Eynde | Numerical Simulations of Transition due to Isolated Roughness Elements at Mach 6[END_REF] exhibits an offset due to different forcing. The work of [START_REF] Van Den Eynde | Numerical Simulations of Transition due to Isolated Roughness Elements at Mach 6[END_REF] showed minor impact of the freestream disturbance amplitude on the growth rate of linear instabilities. This is also observed in the present validation, where the exponential growth rates of all four cases agree very well. The breakdown to turbulence, however, takes place earlier for the SBLI test case due to the increased integrated forcing energy compared to the computations by elsA and OpenSBLI.

VI. Global stability study

A. Numerical methodology

Numerical domain and grid

The same roughness element as in section V was used for this study. The domain was shortened in the streamwise direction to L x = 90 and in the spanwise direction to L z = 40. The roughness element is now placed at a distance of 40 inlet displacement thicknesses from the inflow. The Reynolds number was increased from 14,000 to 28,000 and then to 40,000.

A grid convergence study was made at Re = 40, 000 with OpenSBLI, with a Coarse, Medium and Fine grids having respectively (N x ,N y ,N z ) = (800,100,200), (1200,150,300) and (1800,225,450) grid points in each direction. This showed consistency of the results, in terms of whether or not transition occurs, although refining the grid resulted in earlier transition, as shown by the increase of skin-friction coefficient on figure 9, suggesting that even finer grids are needed. Under these considerations, the number of grid points in each direction was chosen as (N x ,N y ,N z ) = (1800,300,875)

Calculation of the laminar base flow

The numerical solution of the laminar base flow was obtained by using elsA. The convective fluxes were computed using a second order Roe scheme, along with a minmod limiter. Time integration was achieved using an implicit (backward Euler) scheme. Once again the flow is initialized with the compressible similarity solution for the boundary layer. The system is integrated in time until a decrease of sixth order of magnitude in the average residual is achieved. This method was found to be successful in obtaining a steady base flow, even though the Re = 40, 000 undergoes a self-sustained transition to turbulence. Figure 10 shows iso-surfaces of Q criterion (Q = 0.005) coloured by the streamwise velocity for both base flows.

B. Results

Reynolds number comparison

Unsteady simulations were first carried out with OpenSBLI, using a 5 th order TENO scheme for the spatial discretization. Figure 11(a) compares the three Reynolds number in terms of vortical structures induced by the roughness. It is clear that case Re = 14, 000 (i) stays laminar, whereas case Re = 28, 000 (ii) is on the edge to transition. Case Re = 40, 000 (iii), however, is experiencing a self-sustained transition to turbulence. As can be seen in figure 11(b) which is based on a time-averaged flow field, increasing the Reynolds leads to a larger upstream recirculation region (shown by the red line) but a smaller downstream recirculation region in the symmetry plane. The highest Reynolds number cases also exhibits more vortices within the horseshoe vortex system. The time-averaged pressure contours are also shown on figure 11(b).

The Re = 40, 000 case undergoes a self-sustained transition to turbulence, in the sense that no artificial disturbances are needed for instabilities to develop and to sustain by themselves. This is likely to be due to the presence of a global instability with its origin near the roughness element. The fact that the bow shock wrapping around the roughness-element is a potential source of vortex structure unsteadiness was mentioned in [START_REF] Subbareddy | Direct numerical simulation of high-speed transition due to an isolated roughness element[END_REF] and cannot be ruled out. Figure 12 shows an instantaneous of the Q criterion coloured by streamwise vorticity. One can notice the presence of counter-rotating streamwise vortices away from the symmetry plane due to the upstream separation-induced vortices wrapping around the roughness. As described in [START_REF] Chaudhry | Direct Numerical Simulation of Roughness-Induced Transition in the VKI Mach 6 Tunnel[END_REF], the shear layers bounding these streaks roll up and result in hairpin-shaped structures. Close to the symmetry plane and downstream of the roughness, large streamwise streaks are also present, inducing the formation of hairpin-shaped structures. Farther downstream, these trains of hairpin vortices eventually lead to the formation of secondary vortices, resulting in the spanwise spreading of turbulent structures, in a similar way to the results of [START_REF] Iyer | Boundary layer transition in high-speed flows due to roughness[END_REF].

2. Response to a perturbation impulse at Re = 28, 000 and Re = 40, 000

The response to a perturbation impulse was carried out at Re = 28, 000 and Re = 40, 000, following the method described in section III.C. The initial disturbance consists of an increment of 1 × 10 -4 in the density, in a region of 3 × 3 grid points around the location (x, y, z) = (0.33, 0.5, 0.75), corresponding to a location within the boundary layer, upstream of the roughness element. This effectively introduces a disturbance with a sharp-edged distribution, which will excite a range of frequencies at low amplitude. The simulation restarts from the base flow with this initial disturbance and the response of the flow to this perturbation is monitored for four flow-through domains. First of all, the power spectral density (PSD) computed from the wall pressure perturbation response at (x, z) = (88, 0.75) is shown on figure 13. The PSD was adimentionalized by the variance of the signal for better comparison. The non-dimensional frequencies F = 0.14 and F = 0.24 are the most amplified frequencies for Re = 28, 000 and Re = 40, 000 respectively. F = 0.14 was often mentioned in the literature [START_REF] De Tullio | Influence of boundary-layer disturbances on the instability of a roughness wake in a high-speed boundary layer[END_REF] and was related to a varicose instability (or Mack mode). F = 0.24 is also mentioned by [START_REF] De Tullio | Laminar-turbulent transition induced by a discrete roughness element in a supersonic boundary layer[END_REF] as being related to a sinuous mode.

Figure 14 shows the time histories of the wall pressure perturbation p for different streamwise location. On the vertical axis, x indicates the streamwise distance. The maximum amplitude of the pressure perturbation obtained at the end of the domain (x = 90) is ±2.5 × 10 -6 and ±5 × 10 -8 for Re = 28, 000 and Re = 40, 000 respectively. In the Re = 28, 000 case (figure 14(a)), a first wave-packet induced by the initial disturbance can be observed (0 < t < 100). It convects downstream towards the end of the domain. At t = 150, a second wavepacket has reached the end of the domain and undergoes an exponential growth, before converging towards a self-sustained oscillatory state (red line in figure 14(a)). In the Re = 40, 000 case (figure 14(b)), the initial disturbance once again convects downstream, but the exponential growth at x = 90 occurs more rapidly. Even in the absence of disturbances, oscillations self-sustain in the wake of the roughness element, indicating the presence of a global instability. Figure 15 shows the remaining modes at t = 180 by displaying the density perturbation. The structures left at this time show a varicose (symmetrical) mode at Re = 28, 000, whereas a sinuous instability is present for the highest Reynolds number. The shape of the structures is qualitatively similar to the modes found by De Tullio [START_REF] De Tullio | Laminar-turbulent transition induced by a discrete roughness element in a supersonic boundary layer[END_REF] at a lower supersonic Mach number (M=2.5).

VII. Conclusion

The flow over an isolated roughness element has been studied using direct numerical simulations at Mach 6. Comparison of elsA and OpenSBLI was made on a laminar case (Re = 8, 200) over a smooth bump and was found to be successful. This validates the computation of metric and viscous terms in OpenSBLI as well as the ability of 2 nd order schemes in elsA to predict hypersonic laminar flows accurately. Validation of elsA and OpenSBLI with respect to SBLI on a transitional case (Re = 14, 000) where the transition is triggered by acoustic perturbations was also satisfying. The second order AUSM+PW scheme as well as the 6 th order compact scheme in elsA are found to give close results to the 5 th order TENO scheme in OpenSBLI. A global stability study has been carried out using a DNS approach. The flow for three Reynolds numbers was compared without addition of acoustic disturbances. The lowest Reynolds number case stays laminar. At the intermediate Reynolds number, the flow in on the edge to transition but remains laminar within the simulated domain. The largest Reynolds number undergoes a self-sustained transition to turbulence. This is due to the presence of a global instability. The response to an initial perturbation impulse at both the intermediate and large Reynolds numbers shows evidence of a varicose instability (at F = 0.14) and a sinuous instability (at F = 0.24) respectively. The origin of the self-sustained oscillatory mechanism has not yet been determined although the instability appears to be linked to the downstream recirculation region and to the shear-layer induced by the roughness element.

  u, v, w = streamwise, wall-normal, and spanwise velocities W

Fig. 1

 1 Fig. 1 Computational domain used in the current work

Fig. 2

 2 Fig. 2 Pseudo-Schlieren plot showing density gradients in the roughness centreplane. Red region showing the recirculation regions by using contours of small negative streamwise velocity (u = -10 -5 )

  N x x N y x N z 2415 x 205 x 468 1600 x 205 x 468 ∆x [0.15 (upstream), 0.05 (roughness), 0.135 (downstream)] ∆z [0.04 (roughness), 0.15] Table 2 Grid and domain sizes used for the transitional case. Values of ∆x are the grid spacing upstream, near and downstream of the roughness (Lengths are in mm) .

Fig. 3 Fig. 4 Fig. 5 Fig. 6

 3456 Fig. 3 OpenSBLI grid convergence study in the symmetry plane

Fig. 7 Fig. 8

 78 Fig. 7 Instantaneous Q criterion iso-surfaces (Q = 0.003) for OpenSBLI (left) and elsA (right), coloured by the streamwise vorticity

Fig. 9 Fig. 10

 910 Fig. 9 Grid convergence study for the Re = 40, 000 case -Skin friction coefficient averaged over span and time

  (a) Isosurface of Q criterion Q = 0.05 (b) Time-average pressure in the symmetry plane. The red line shows separated regions.

Fig. 11

 11 Fig. 11 Reynolds number comparison. From top to bottom, Re = 14, 000 (i), Re = 28, 000 (ii) and Re = 40, 000 (iii)

Fig. 12 Fig. 13 PowerFig. 14 6 Fig. 15

 121314615 Fig. 12 Visualization of instantaneous Q criterion coloured with streamwise vorticity w x at Re = 40, 000

th order Runge-Kutta Table 1 Numerical parameters for each code used in this study Table

			,
			r
	where S * = 110.4K is the Sutherland reference temperature and T * r = 273.15K
	B. Numerical method
		elsA	OpenSBLI	SBLI (Legacy)
		[ONERA]	[University of Southampton]	[University of Southampton]
		Finite volume	Finite difference	Finite difference
	Convective fluxes AUSM+PW (2 nd order)	5 th order WENO	4 th order central difference
		Lele (6 th order)	5 th order TENO	+ Total Variation Diminishing (TVD)
	Viscous fluxes	5 points centred stencil	4 th order central difference	4 th order central difference
	Time integration	2 nd order Gear	3 r d order Runge-Kutta	3 r d order Runge-Kutta
		4

Acknowledgments

This work was funded by the "Ministère des armées -Direction Générale de l'Armement (DGA/DSTL)" as part of the UK-FR PhD program. It was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (http://www.csd3.cam.ac.uk/), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council.


Adrien Kassel
email: adrien.kassel@ens-lyon.fr

Thierry L Évy

Thierry Lévy
email: thierry.levy@sorbonne-universite.fr

A colourful

Keywords: 2010 Mathematics Subject Classification. 05C30, 05C22, 15A15 matrix-tree theorem, graph, forests, cycles, Laplacian, determinant, Q-determinant, holonomy, ordered products, simplicial complexes, pseudoforests, circular and bicircular matroids

come

Introduction

It is well known in combinatorics that determinants of matrices can be interpreted as weighted sums over paths in graphs. One famous and prominent example is the matrix-tree theorem and its variants which interpret minors of a discrete Laplacian as weighted counts of spanning forests. A beautiful proof of its most classical version was given by Doron Zeilberger in [Zei85, Section 4], a short paper which emphasizes this combinatorial approach to matrix algebra, a point of view which seems to have been newer and less popular at the time of its publication.

In the course of our research in probability and stochastic processes on graphs with gauge symmetry [Kas15, KL16, KL19], we were led to reconsider this classical theorem, adding in new parameters in the guise of matrices attached to edges of a directed graph. While doing so we revisited Zeilberger's proof, in its memorable colourful re-enactment performed for us by David Wilson some years ago, and found a concise and pictorial way to prove in a unified way several such theorems, notably the classical one of Kirchhoff [START_REF] Kirchhoff | Über die Auflösung der Gleichungen, auf welche man bei der Untersuchung der linearen Vertheilung galvanischer Ströme geführt wird[END_REF], the signed version of Zaslavsky [START_REF] Zaslavsky | Signed graphs[END_REF]Theorem 8A.4], the abelian-group-weighted version of Chaiken [Cha82, Eq. (15) and the following unnumbered equation], a special case of which is the complex-weighted one of Forman [START_REF] Forman | Determinants of Laplacians on graphs[END_REF] Theorem 1], and the quaternion case of Kenyon [Ken11,Theorem 9]. Let us emphasize that this proof is one of those which does not use the Cauchy-Binet theorem, a step which always blurs somewhat the understanding of cancellations in the signed combinatorial sums involved.

We believe this proof to be crisp and clear enough to be shared with the combinatorics community, and hope it can publicise Zaslavsky's, Chaiken-Forman's, and Kenyon's apparently lesser known theorems. We also believe this approach may inspire the discovery of new formulas and we indeed provide at least a few new identities, which although not definitive, show some insight into the combinatorics of graphs with parameters and the power of Zeilberger's approach.

The note is organised as follows. We first recall definitions about determinants, then proceed to give our revisited and unified proof. Last, we prove new results and discuss in an open-ended way some natural perspectives.

A unified shorthand notation for determinant and Q-determinant

Let R be a ring and S a commutative ring. Let τ : R → S be a morphism of Z-modules, that is, an additive map. Let us assume that τ is central, in the sense that for any r, r ∈ R, we have τ (rr ) = τ (r r). For all integers n 1 and all n × n matrices M ∈ M n (R), we define the τ -determinant of M as the element of S given by the formula

(1) det τ (M ) = σ∈Sn ε(σ) c cycle of σ c=(i 1 ...ir) τ (M i 1 i 2 M i 2 i 3 . . . M iri 1 ) ,
where S n denotes the symmetric group on n elements. We do not claim that the function det τ has particularly good properties on M n (R) in general. In particular, although it is Z-multilinear in the columns and the rows of M , it is not alternating in general. However, its interest for us is that it unifies on the one hand the usual notion of determinant of a matrix with complex entries, or with entries in a ring of polynomials with complex coefficients, or indeed with entries in any commutative ring, by taking S = R and τ = id R ; and on the other hand the less usual notion of Q-determinant [Moo22] (Moore's original work can only be found easily as a conference abstract, so we refer to [Meh04, Section 5.1] for a more detailed textbook treatment) of a matrix with quaternionic entries, by taking R = H, S = R and τ = , the real part map. It also applies to the situation where R is a polynomial ring with quaternionic coefficients, and S the polynomial ring with the same indeterminates and real coefficients.

A matrix-tree theorem

Let n 2 be an integer and let V = {1, . . . , n} be the set of vertices of a complete graph. Let E be the set of edges of this graph, that is, the set of ordered couples of distinct vertices: E = {(i, j) ∈ V 2 : i = j}. We think of the edge (i, j) as going from i to j. 1 Let H be a ring, which we do not assume to be commutative. Examples that we have in mind are R, C, or H. Let R be the polynomial ring H[a ij : (i, j) ∈ E] over n(n -1) indeterminates. Let {h ij , (i, j) ∈ E} be n(n -1) elements of H. As the notation suggests, we think of a ij and h ij as being attached to the edge (i, j).

We form the matrix ∆ ∈ M n (R) by setting, for all (i, j) ∈ E,

(2)

∆ ij = -h ij a ij
and, for all i ∈ V,

(3) ∆ ii = j∈V\{i} a ij .
The use of a colour will play an instrumental role in the computation, by helping us to keep track of the origin of the coefficients. 2 Let us choose a commutative ring K and a central morphism of Z-modules τ : H → K. Here, the examples we have in mind for K are R or C, and τ the identity map or the real part. Let us define S = K[a ij : (i, j) ∈ E]. The morphism τ extends to a morphism R → S that we still denote by τ .

1 Since our edges will be weighted, all results can be easily re-interpreted to take into account graphs with multiple (weighted) edges. We can also allow self-loops by adding the diagonal coefficients aii(1 -hii).

2 For those reading this paper in black and white, note that the right-hand side of Eq. (2) is red; moreover, in the rest of the text, replace the word red by marked. However, for those who can afford it, we recommend the coloured version of the paper.

Let us fix an integer m ∈ {1, . . . , n}. We consider the m × m principal submatrix ∆ [m] of ∆ obtained by erasing the rows and columns of indices greater than m. The version of the matrix-tree theorem that we are going to state computes the τ -determinant of ∆ [m] .

Let us define the boundary of our graph, or its well, as the subset W = {m + 1, . . . , n} ⊂ V. Let us also define U = {1, . . . , m} = V \ W, the set of inner vertices.

We say that a subset F ⊂ E is a cycle-and-well-rooted spanning forest (with well W) if it contains exactly m edges, one coming out of i for each i ∈ U (see Figure 1). 3 We denote by F → m the set of all cycle-and-well-rooted spanning forests with well W. The arrow is meant to remind us that each element of F → m is a configuration of directed edges. Consider F ∈ F → m . We set a F = (i,j)∈F a ij . Moreover, using the decomposition into connected components of the undirected graph underlying (V, F), we can partition the directed graph (V, F) into disjoint pieces, each of which is of one of two kinds: either a tree rooted at a vertex of W with all edges consistently pointing towards the root, or a cycle-rooted tree (also called a unicycle) contained in U. To F, we attach the collection C (F) of the cycles of its unicycles. It is a (possibly empty) family of pairwise disjoint simple cycles in U. For each such cycle c, visiting the vertices i 1 , . . . , i r in this cyclic order, we can form the ill-defined element h c = h i 1 i 2 . . . h iri 1 of H, and the well-defined element τ (h c ) = τ (h i 1 i 2 . . . h iri 1 ) of K.

Theorem 3.1. In the ring S = K[a ij , (i, j) ∈ E], we have the equality

det τ ∆ [m] = F∈F → m a F c∈C (F) 1 -τ (h c ) .
Proof. Let us write the definition (1) of the τ -determinant of ∆ [m] and separate, in the contribution of each permutation, the fixed points from the cycles of length at least 2. We have

det τ ∆ [m] = σ∈Sm ε(σ) c=(i 1 ...ir) τ (∆ i 1 i 2 . . . ∆ iri 1 ) = σ∈Sm ε(σ) i: σ(i)=i τ (∆ ii ) c=(i 1 ...ir), r 2 τ (∆ i 1 i 2 . . . ∆ iri 1 ) . (4)
For each permutation σ ∈ S m , we will interpret the term of this sum indexed by σ as a sum of a certain number of contributions, each being attached to a certain configuration of black and red 4 edges on our graph. There will be, for each permutation σ, a number of edge configurations equal to (n -1) f , where f is the number of fixed points of σ. Each edge configuration will be a cycle-and-well-rooted spanning forest, in which some of the cycles are made of red edges.

To be more precise, let us consider a particular permutation σ. For each cycle c = (i 1 . . . i r ) of σ with length at least 2, we have

τ (∆ i 1 i 2 . . . ∆ iri 1 ) = (-1) r τ (h i 1 i 2 . . . h iri 1 )a i 1 i 2 . . . a iri 1 .
We think of this term as associated with a cycle of red edges on our graph, namely the cycle formed by the edges (i 1 , i 2 ), . . . , (i r , i 1 ).

Then, for each fixed point i of σ, we have

τ (∆ ii ) = τ (a i1 + . . . + a ii + . . . + a in ) = a i1 + . . . + a ii + . . . + a in .
We think of these n -1 terms as being respectively associated with the n -1 possible ways of drawing one edge in our graph coming out of the vertex i. This edge, as the notation suggests, is a black edge.

3 These combinatorial objects or their variants are also called mappings, pseudoforests, or bases of the bicircular matroid in the literature. 4 Or marked.

Expanding the term indexed by σ in (4) over the sum produced by each fixed point of σ makes it appear as the sum of (n -1) f contributions, where f is the number of fixed points of σ, each contribution being attached to a configuration of black and red edges on our graph. Each edge configuration contains exactly one edge coming out of each vertex of U = {1, . . . , m} and is thus, in the terminology introduced above, a cycle-and-well-rooted spanning forest. Moreover, in each cycle-and-well-rooted spanning forest that appears, some of the cycles are made of red edges. Here the phrase some of the cycles is a shorthand for an arbitrary, possibly empty and possibly full, subset of the set of cycles.

Let us call cycle-coloured cycle-and-well-rooted spanning forest a configuration of black and red edges on the graph of the kind that we just described, and illustrated in Figure 1 We claim that each cycle-coloured cycle-and-well-rooted spanning forest appears once and only once in the whole sum (4). Indeed, consider such an edge configuration. Form the permutation σ of {1, . . . , m} of which the cycles of length at least 2 are exactly the red cycles of the edge configuration. The fixed points of σ are thus the vertices of U out of which comes a black edge. Then this particular edge configuration appears exactly once in the term of (4) indexed by that permutation σ, and in no other.

Finally, the contribution of a cycle-coloured cycle-and-well-rooted spanning forest to (4) is the product of

• the indeterminates a ij attached to the edges of the configuration, regardless of their colour, • a term -τ (h c ) for each red cycle c of the configuration. The minus sign which appears in the contribution of a red cycle of length r is the product of a sign (-1) r coming from the fact that the off-diagonal entries of ∆ carry a minus sign, and the contribution (-1) r-1 of this cycle to the signature of σ.

At this point, we have replaced the sum (4) by a larger sum indexed by (black and red) cycle-coloured cycle-and-well-rooted spanning forests. We are now going to reorganise the terms of this larger sum and make it into a sum over (black 5 ) cycle-and-well-rooted spanning forests. Indeed, let us consider a cycle-and-well-rooted spanning forest F ∈ F → m . We will assign to F the sum of all the contributions coming from the various ways of colouring none, some, or all of the cycles of F in red. There are 2 b contributions, where b is the number of cycles of F, and their sum is equal to the product of

• the indeterminates a ij attached to the edges of F, that is, a F , • a term 1 -τ (h c ) for each cycle c of F, accounting for the fact that this cycle can be left black or coloured in red. This is exactly the announced formula.

In many applications, the indeterminates {a ij , (i, j) ∈ E} are specialised in a symmetric way, that is, in such a way that the equality a ij = a ji is satisfied for each edge (i, j) ∈ E. Theorem 3.1 can be nicely adapted to this situation. Indeed, in this symmetric case, the monomial a F attached to a cycle-and-well-rooted spanning forest F is not affected by reversing the orientation of one or several cycles of F. To acknowledge this invariance, let us declare equivalent two elements of F → m which differ only by the orientation of some of their cycles, and let us denote by F m the quotient set.

Let [F] be the class in F m of a cycle-and-well-rooted spanning forest F. The obvious statement that the elements of [F] are deduced from F by changing the orientation of some of the cycles of F comes with the nice twist that the two orientations of a cycle of length 2 are in fact one and the same. For this reason, if c is a cycle of F of length 2, then c = c -1 , hence τ (h c ) = τ (h c -1 ) is well defined and depends only on [F]. On the other hand, if c is a cycle of F of length at least 3, then it is the pair {τ (h c ), τ (h c -1 )} that depends only on [F].

Summing the expression of det τ ∆ [m] provided by Theorem 3.1 over equivalence classes yields the following result. We denote the length of a cycle c by (c).

Corollary 3.2. In the quotient ring S/(a ij -a ji : (i, j) ∈ E), we have the equality

det τ ∆ [m] = [F]∈Fm a F c∈C (F) (c)=2 1 -τ (h c ) c∈C (F) (c) 3 2 -τ (h c ) -τ (h c -1 ) .

Special cases

4.1. The classical matrix-tree theorem. Consider H = K = R and let τ be the identity map. Take h ij = 1 for each edge (i, j). Then in the sum in Theorem 3.1, non-zero contributions arise only from cycle-and-well-rooted spanning forests without any cycle, that is, from spanning forests rooted in the well. We recover exactly the classical matrix-tree theorem.

4.2. Forman's theorem. Consider H = K = C and let τ be the identity map. Take h ij to be any complex number for each edge (i, j). Then we recover Forman's matrix-tree theorem with holonomies6 along the edges that are complex numbers ([For93, Theorem 1]). If the variables a ij are symmetric, and if the holonomies are taken to be complex numbers such that h ji = h ij for each edge (i, j), then we can apply Corollary 3.2 to find that cycles of length 2 contribute by a factor 1 -(h c ) and cycles of length at least 3 by a factor 2 -2 (h c ). Assuming moreover that the complex numbers h ij have modulus 1 forbids the presence of cycles of length 2 in any of the cycle-and-well rooted forests in the sum, since in this case, h c = 1 for any such cycle.

4.3. Chaiken's theorem. Consider H = K = Z[Γ],
where Γ is an abelian group, and τ the identity map. Then Theorem 3.1 is the formula in the middle of p.326 of [START_REF] Chaiken | A combinatorial proof of the all minors matrix tree theorem[END_REF] in the case of principal minors. By further specialising Γ to be the multiplicative group {±1}, we recover Zaslavsky's theorem about signed graphs ([Zas82, Theorem 8A.4]), also a special case of Forman's theorem.

4.4. Kenyon's theorem. Consider H = H, K = R and let τ be the real part map. Choose the quaternions h ij to be of norm 1 for each edge (i, j), and such that h ji = h ij . Then, up to the identification of the unit sphere of H with SU(2), Corollary 3.2 reduces to Kenyon's theorem with holonomies in SU(2) ([Ken11, Theorem 9], which Kenyon stated in the case of an empty well only, although his proof extends readily to the case of a non-empty well).

by the setup where the graph is embedded in a manifold and where hij is the parallel transport of a connection along edge ij.

Graphs with holonomy in a space of matrices

In our research about stochastic processes on graphs with gauge symmetry [Kas15, KL16, KL19], our main motivation was to understand the determinant of ∆ [m] when the elements h ij are matrices, for instance unitary matrices. This is an ill-posed question in the sense that several expressions, with various merits and demerits, can be given for this determinant.

For example, if we are interested in N × N complex matrices, we can apply Theorem 3.1 with H = M N (C), K = C and τ (•) = 1 N Tr(•). One difficulty with the formula that we obtain in this way is that it seems difficult to make sense of the definition det τ ∆ [m] . We propose an alternative approach, in which we imitate and adapt the proof of Theorem 3.1.

Let us describe the situation again. We have a ring H, not assumed to be commutative, a commutative ring K and a central morphism τ : H → K. We have the set of vertices V = {1, . . . , n}, in which a certain subset W = {m + 1, . . . , n} is called the well. To each edge (i, j) ∈ E is attached an indeterminate a ij and an N × N matrix h ij ∈ M N (H). We form the matrix ∆ following (2) and (3), with a ij understood as a ij I N in (3). It is an n × n matrix of N × N matrices, which we see in the most natural way as an N n × N n matrix with entries in H.

We will compute the τ -determinant of the N m × N m matrix ∆ [N m] obtained from ∆ by erasing the rows and columns of indices greater than N m.

For this, we introduce, above our usual graph G = (V, E), a new graph which we denote by G N . The set of vertices of this new graph is V × {1, . . . , N }. Thus, a vertex of G N is a pair (i, k) where i ∈ {1, . . . , n} is a vertex of G and k ∈ {1, . . . , N } can be used to index a row or a column of an N × N matrix. The edges of G N are all pairs ((i, k), (j, l)) such that (i, j) is an edge of G and (k, l) ∈ {1, . . . , N } 2 . It is an important feature of the graph G N that it contains no 'vertical' edge (i.e. an edge of the form ((i, k), (j, l)) To each edge ((i, k), (j, l)) of G N , we associate on the one hand the indeterminate a ij and on the other hand the element h ((i,k),(j,l)) = (h ij ) kl of H, that is, the (k, l)-th entry of the matrix h ij .

with i = j). 7 G G 5 G G 5
We define the well of the graph G N as the set W × {1, . . . , N } of all vertices located 'above' a vertex of W. Our formula involves a particular subset of the set F → m,N of cycle-and-wellrooted spanning forests in G N . To describe this set, let us define a horizontal edge as an edge ((i, k), (j, l)) such that k = l. Then, we say that an element F of F → m,N is horizontal if every edge of F that is not in a cycle is horizontal. We denote by H F → m,N the corresponding subset of F → m,N . Moreover, we say that a cycle of F is horizontal if it consists only of horizontal edges. We say that a cycle is skew if it is not horizontal. We denote by H (F) the set of horizontal 7 We are not aware of a standard name for this simple construction, but as pointed out to us by the referee, a more combinatorially idiomatic characterisation is that G N is the lexicographic product of G with the edgeless graph on N vertices. cycles of F, and by S (F) the set of its skew cycles. For every cycle c of F, we define h c as the product of the complex numbers h ((i,k),(j,l)) associated with the edges which constitute c.

Our result is the following.

Theorem 5.1. In the ring S = K[a ij , (i, j) ∈ E], we have the equality

det τ ∆ [N m] = F∈H F → m,N a F c∈H (F) 1 -τ (h c ) c∈S (F) -τ (h c ) .
Proof. The proof is very similar to the proof of Theorem 3.1. The main difference is that we see the N m indices of the rows and of the columns of ∆ [N m] as corresponding to the N m vertices of the new graph G N : the vertex (i, k) corresponds to the (N (i -1) + k)-th row, or column. We expand the determinant of ∆ [N m] as a sum indexed by permutations, which we really think of as permutations of the set {1, . . . , m} × {1, . . . , N }. We separate, for each permutation, the fixed points from the non-trivial cycles.

The contribution of each non-trivial cycle is understood as being associated with a cycle of red edges on the graph G N . This is possible because the diagonal N × N blocks of ∆ are diagonal n × n matrices (recall the third paragraph of the current Section 5), so that the contribution of any permutation containing a cycle which would produce a 'vertical' edge vanishes. To be clear, the contribution of a permutation σ such that there exists i ∈ {1, . . . , m} and k, l distinct elements of {1, . . . , N } such that σ(i, k) = (i, l) is zero, because ∆ (i,k),(i,l) = 0.

The fixed points contribute by producing a sum over all ways of completing the configuration of red cycles that we just obtained, by adding a horizontal black edge coming out of every vertex of {1, . . . , m} × {1, . . . , N } that is not already visited by a red cycle.

Let us emphasize that there is something arbitrary in our decision of assigning the contributions of fixed points to horizontal edges (we come back to this point in Section 7.1). Let us at least spell out what this exactly means: if (i, k) is a fixed point of σ, then in expanding the sum a i1 + . . . + a ii + . . . + a in , we will attach, for each j = i, the weight a ij to the edge ((i, k), (j, k)).

The outcome of this procedure is an expression of the determinant of ∆ [N m] as a sum over cycle-and-well-rooted spanning forests in G N , with well {m + 1, . . . , n} × {1, . . . , N }. Thanks to our decision of assigning the contributions of fixed points to horizontal edges, all the forests which appear are horizontal, in the sense that every edge that is not in a cycle is horizontal. Moreover, each of these forests has a certain number of its cycles made of red edges. To complete our description, we must say that any cycle that is not horizontal must be made of red edges, because all black edges are horizontal.

Just as in the proof of Theorem 3.1, we can check that every horizontal cycle-and-well-rooted spanning forest in G N with all its skew cycles and an arbitrary subset of the set of its horizontal cycles made of red edges contributes exactly once to the sum, with a contribution equal to the product of

• the indeterminates a ij attached to the edges of the configuration, regardless of their colour, • a term -τ (h c ) for each red cycle c. Reorganising the sum by summing first over sets of configurations which differ only by the colour of their edges, and taking into account the fact that skew cycles are always red, whereas horizontal cycles can be either black or red, we obtain the announced formula.

In the case where the indeterminates {a ij , (i, j) ∈ E} are symmetric, the discussion of the end of Section 3 applies and, by summing over equivalence classes of cycle-and-well-rooted spanning forests and using self-explanatory notation, we obtain the following result.

Corollary 5.2. In the quotient ring S/(a ij -a ji : (i, j) ∈ E), we have the equality

det τ ∆ [N m] = [F]∈H F m,N a F c∈H (F) (c)=2 1 -τ (h c ) c∈H (F) (c) 3 2 -τ (h c ) -τ (h c -1 ) c∈S (F) (c)=2 -τ (h c ) c∈S (F) (c) 3 -τ (h c ) -τ (h c -1 ) .
Let us stress that, just as in Theorem 5.1, cancellations may occur in the right-hand side of the equation of Corollary 5.2. When a ij = a ji , we can easily factor ∆ as a product of rectangular matrices. Provided we are in a case where it holds for the τ -det8 , the Cauchy-Binet formula implies that the coefficients of the monomials a F such that the degree of a ij exceeds N are zero. However, we do not know a direct proof of this fact. 6. Applications 6.1. The case N = 1. We recover Theorem 3.1 and Corollary 3.2. 6.2. The unitary case. When H = R, C or H and the matrices h ij are unitary, we have for any cycle c ∈ G N that τ (h c ) and τ (h c -1 ) are complex conjugates so that we can slightly simplify the equation in Corollary 5.2 using real parts, writing 2 (τ (h c )) in place of τ (h c ) + τ (h c -1 ). Moreover, in that case, the Cauchy-Binet formula alluded to right above shows that the coefficient of each monomial a F is nonnegative. 6.3. The case of block-diagonal matrices. If the matrices h ij are block-diagonal in a consistent way, then the determinant factorises, a fact that is clear from the point of view of the left-hand side of the formula. 6.4. Simplicial complexes. Analogs of the classical matrix-tree theorem have been proven for simplicial complexes, starting with Kalai [START_REF] Kalai | Enumeration of Q-acyclic simplicial complexes[END_REF] and followed by works of Adin [START_REF] Adin | Counting colorful multi-dimensional trees[END_REF], Lyons [START_REF] Lyons | Random complexes and 2 -Betti numbers[END_REF] and Duval-Klivans-Martin [START_REF] Duval | Simplicial matrix-tree theorems[END_REF] among others.

Recall that a simplicial complex C on the vertex set C 0 = {1, . . . , v} is a collection of nonempty subsets of C 0 , called cells, which contains the singletons {1}, . . . , {v} and which contains every non-empty subset of every set that it contains.

For concreteness we consider here the simplicial complex C of all non-empty subsets of C 0 of cardinality less than or equal to d + 1 for some fixed d ∈ {1, . . . , v -1}. For any k ∈ {0, . . . , d}, it contains v k+1 subsets of cardinality k + 1, which we call cells of dimension k, or k-cells. A graph is a simplicial complex of dimension 1, and its Laplacian is an operator acting on 0chains. In higher dimension, the analogs of the matrix-tree theorem concern a natural operator on (d-1)-chains, the so-called 'top-down' Laplacian (coming from the structure of chain complex with boundary operator ∂ and coboundary operator δ and defined as ∆ = ∂δ). Here we consider a version of this operator twisted by elements of a ring H. For the purposes of this short paper, and in tune with Zeilberger's opening lines of [START_REF] Zeilberger | A combinatorial approach to matrix algebra[END_REF], we do not detail the definition of this operator from the algebraic-topological point of view of the chain complex, but rather directly define it combinatorially as a matrix indexed by the set of (d -1)-cells. This matrix depends on a choice of orientation of the cells, but in a way that does not affect its principal minors.

Let C d-1 be the set of all (d -1)-cells of C. We call two elements σ, τ of C d-1 adjacent, and we write σ ∼ τ , if σ ∪ τ is a d-cell. We denote by E the set of ordered pairs of adjacent (d -1)-cells in C. To each (σ, τ ) ∈ E we assign an indeterminate a στ and an element h στ of H.

We fix an arbitrary choice of orientation for each (d -1)-cell. This choice determines for each pair σ ∼ τ a sign ε στ , the exact definition of which need not worry us here. The important fact is that for any chain c = (σ 1 . . . σ k ) of adjacent (d -1)-cells with σ k ∼ σ 1 , the product of the signs ε σ j σ j+1 when j runs cyclically through {1, . . . , k}, denoted by ε c , is independent of the choice of orientation of the cells σ j . 9 We can now write the non-zero entries of the 'matrix' 10 ∆ by setting, for all (σ, τ ) ∈ E,

(5) ∆ στ = -ε στ h στ a στ and for all σ

∈ C d-1 , (6) ∆ σσ = 1 d τ ∈C d-1 \{σ} τ ∼σ a στ .
This operator ∆ can indeed be written in the form ∂δ whenever the weights satisfy some symmetry conditions, namely: for σ ∼ τ , the weight a στ is only a function of the d-cell ρ = σ ∪ τ and moreover there exist weights h σρ and h ρτ such that h στ = h σρ h ρτ .

The similarity between Eqs. ( 5) & (6) and Eqs. ( 2) & (3) suggests to look at ∆ as a Laplacian matrix on the weighted graph G = (C d-1 , E) with weights ãστ = a στ /d and hστ = d ε στ h στ , when σ ∼ τ ; and ãστ = 0 and hστ = 1, otherwise.

We are now exactly in the setup of Section 5 and we can thus use our Theorem 5.1 to compute the principal minors of ∆. For simplicity, we only state the N = 1 case, which is thus a corollary of Theorem 3.1. Fix m ∈ {1, . . . , v d }. Theorem 6.1. In the ring S = K[a ij , (i, j) ∈ E], we have the equality

det τ ∆ [m] = 1 d m F∈F → m a F c∈C (F) (1 -d (c) ε c τ (h c )) .
As in the discussion at the end of Section 5, an application of the Cauchy-Binet formula (when it holds, see Footnote 8) in the case where the weights are symmetric enough that ∆ = ∂δ (which means in particular that for all σ ∼ τ , the weight a στ is a function x ρ of the d-cell ρ = σ ∪ τ ) shows the cancellation of many terms in the sum: all terms for which the monomials a F contain x ρ with degree exceeding 1. 11 Moreover the operator ∆ then has a non-zero kernel (the image of the coboundary operator δ on (d-2)-cells) so all minors cancel for m large enough. In particular, choosing h = 1 and m = v d , the determinant is zero, something which is not clearly apparent from the right-hand side of the formula.

The case h = 1 has been studied considerably, and in that case, the coefficient of each monomial x F in the maximal non-zero principal minors of ∆ is the square of the torsion of a certain (d -1)-subcomplex called a simplicial spanning tree. 12 Our formula expresses this torsion in terms of a weighted count of cycles. This can be compared to the formulas of [START_REF] Bernardi | Directed rooted forests in higher dimension[END_REF] expressing the torsion in terms of related combinatorial sums involving cancellations.

Are there direct proofs of these identities obtained by comparing two polynomial expansions? In particular, can one compare our Theorem 6.1 with Theorem 31 and Remark 32 of [START_REF] Bernardi | Directed rooted forests in higher dimension[END_REF] 9 As pointed out to us by the referee, the above combinatorial setup can be extended from simplicial complexes to finite CW complexes by allowing the signs εστ to become arbitrary integers. See [START_REF] Duval | Simplicial and cellular trees[END_REF] for a survey on spanning trees in the cellular case.

10

To be precise, we did not specify an ordering of the indices so it is only a matrix up to permutation, a fact which bears no consequence on the computation of its principal minors.

11 In that special case, the equation of Theorem 6.1 holds with the sum restricted to a certain subset of F → m . 12 By torsion, we mean the cardinality of a certain torsion subgroup of the homology of that subcomplex.

and go from one statement to the other? Our statement is yet another evidence that homological quantities may be expressed in geometrico-combinatorial terms as shown by Bernardi and Klivans. We leave these open questions to the interested reader. 7. Loose ends 7.1. Horizontal edges. As noted in the course of the proof of Theorem 5.1, there was something arbitrary about interpreting the black edges as horizontal.

By averaging, in the proof of Theorem 5.1, over all possible inclinations that we can assign to black edges (not merely the horizontal one), or equivalently by applying Theorem 3.1 to the graph G N , we obtain a variant of Theorem 5.1. The main difference is that we now sum over all cycle-and-well-rooted spanning forests of G N . There are many more terms in the sum (which thus cancel out somehow), but the formula is slightly slicker. We leave the easy details to the interested readers.

Theorem 7.1. With the notation of Section 5, we have, in the ring S = K[a ij , (i, j) ∈ E], the equality

det τ ∆ [N m] = 1 N N m F∈F → m,N a F c∈C (F) 1 -N (c) τ (h c ) .
We note the close similarity between the identities in Theorem 7.1 and Theorem 6.1 about simplicial complexes. This suggests a maybe not so surprising combinatorial analogy between the study of vector bundles of rank 1 over simplicial complexes of dimension d with v = dn vertices, and the study of vector bundles of rank N = d over graphs with n vertices.

This similarity also suggests that there could exist a formula for simplicial complexes of the more concise form of Theorem 5.1. 7.2. Gauge invariance. One shortcoming of our Theorem 5.1 is that the terms in the sum of the right-hand side are not manifestly gauge-invariant, nor non-negative when holonomies are unitary. We propose solutions to this problem in an upcoming work. 7.3. Non-principal minors. One question we have not adressed in this note is the computation of non-principal minors of discrete Laplacians. In the classical case, it is known that these minors also enumerate spanning trees.

Chaiken [START_REF] Chaiken | A combinatorial proof of the all minors matrix tree theorem[END_REF]p. 326] computed a formula for non-principal minors of so-called 'gaingraphs/voltage-graphs' which in our language are graphs with holonomy taking values in an abelian group.

It is possible to work out formulas in our case too, and we leave this endeavour to the curious reader.

  below.

Figure 1 .

 1 Figure 1. A cycle-coloured cycle-and-well-rooted spanning forest. The vertices belonging to the well are white.

Figure 2 .

 2 Figure 2. The new graph has N 2 directed edges over each directed edge of the old graph.

We call holonomies the elements hij and their products along paths in the graph. This terminology is inspired

This is the case when H is commutative and τ = id or when H = H, h is unitary, and τ = .

Acknowledgments. We thank David Wilson for teaching us Zeilberger's proof in an enlightening way using coloured chalks on a blackboard of the Newton Institute in Cambridge during the Random Geometry programme in 2015. Unknowingly, it planted the seed for the above proofs. We also thank Omid Amini and Grégory Miermont for stimulating and inspiring discussions which encouraged us to write this note. We are also grateful to Omid Amini for his help in checking Theorem 6.1. Finally, we thank the anonymous referee for a careful reading and helpful feedback.

A.K. thanks CNRS for support through a PEPS grant in 2018 and the Newton Institute for hospitality in 2015 during the programme Random Geometry. T.L. would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme Scaling limits, rough paths, quantum field theory when work on this paper was pursued. This work was supported by EPSRC grant numbers EP/K032208/1 and EP/R014604/1.