M Meldi

J C Vassilicos

Analysis of Lundgren's matched asymptotic expansion approach to the Kármán-Howarth equation using the EDQNM turbulence closure

In this paper we investigate whether the features of the non-equilibrium cascade, which have been identified in recent studies using high-fidelity tools, can be captured in the case of the classical dissipation scaling by turbulence closures based on the statistical description of freely decaying isotropic turbulence. Numerical results obtained using the EDQNM model over a very large range of Reynolds numbers (from Re λ = 50 up to Re λ = 10 6 ) are analyzed to perform an extensive investigation of the scaling region identified as inertial range in Kolmogorov's theory. It is observed that EDQNM results are in agreement with the results of Lundgren's matched asymptotic expansion approach to the Kármán-Howarth equation. Both predict that the Kolmogorov inertial range equilibrium is never obtained irrespective of Reynolds number. Equilibrium is reached in the vicinity of the Taylor length λ (which depends on viscosity) as Reynolds number tends to infinity and there is a gradual departure from equilibrium as the length scale moves away from λ, in particular towards scales larger than λ all the way to the integral length-scale.

I. INTRODUCTION

A new turbulence dissipation scaling has been discovered over the past 10 years ( [START_REF] Vassilicos | Dissipation in turbulent flows[END_REF][START_REF] Cafiero | Non-equilibrium turbulence scalings and self-similarity in turbulent planar jets[END_REF] and references therein) which characterises the non-equilibrium cascade in the relative nearfield of various turbulent flows: grid-generated turbulence, turbulent jets, turbulent wakes and periodic turbulence (both forced and decaying). This relative near-field is in fact quite extensive in practical terms and the turbulence dissipation scaling in that relative near-field incorporates an explicit dependence on inlet/initial conditions. This scaling can be written as ε ∼ U 0 L 0 K/L 2 where ε is the energy dissipation rate, K is the turbulent kinetic energy and L is the integral length scale. On the other hand, U 0 and L 0 are a velocity scale and a length scale associated with initial/inlet conditions. Further downstream one expects no such explicit dependence, and a turbulence dissipation scaling which is the classical Taylor-Kolmogorov scaling ε ∼ K 3/2 /L. [START_REF] Vassilicos | Dissipation in turbulent flows[END_REF][START_REF] Taylor | Statistical Theory of Turbulence[END_REF][START_REF] Kolmogorov | On the degeneration of isotropic turbulence in an incompressible viscous fluid[END_REF][START_REF] Antonia | Approach to the 4/5 law in homogeneous isotropic turbulence[END_REF]. Such a transition from one dissipation scaling to another has indeed been observed in some turbulent flows [START_REF] Vassilicos | Dissipation in turbulent flows[END_REF][START_REF] Goto | Unsteady turbulence cascades[END_REF]. However, the recent work of Goto & Vassilicos [START_REF] Goto | Unsteady turbulence cascades[END_REF] suggests that the far-field dissipation scaling is not necessarily the result of an equilibrium turbulence cascade even if it takes the exact same form as it would take in an equilibrium cascade situation.

Obligado & Vassilicos [START_REF] Obligado | The non-equilibrium part of the inertial range in decaying homogeneous turbulence[END_REF] explained why there is no equilibrium cascade at length-scales large enough to justify a Taylor-Kolmogorov dissipation scaling in the case of freely decaying homogeneous isotropic turbulence (HIT) with classical Taylor-Kolmogorov dissipation scaling. The starting point is Von Karman & Howarth [START_REF] Von Karman | On the statistical theory of isotropic turbulence[END_REF] equation in the form that it takes when expressed for structure functions as in Landau & Lifschitz [START_REF] Landau | Fluid mechanics: Vol 6 of course of theoretical physics[END_REF]:

- 2 3 ε - 1 2 ∂ ∂t S 2 = 1 6r 4 ∂r 4 S 3 ∂r - ν r 4 ∂ ∂r (r 4 ∂S 2 ∂r ) ( 1
)
where ε is the turbulence dissipation rate, ν is the kinematic viscosity of the fluid, r is the distance between the pairs of points defining the second and third order longitudinal structure functions S 2 and S 3 respectively. Using Lundgren's [START_REF] Thomas S Lundgren | Kolmogorov two-thirds law by matched asymptotic expansion[END_REF] analysis of the Kármán-Howarth equation 1 based on matched asymptotic expansions and wind tunnel data for grid-generated approximate HIT conditions, they demonstrated that spectral equilibrium can only be achieved in the vicinity of the Taylor length-scale when the Reynolds number tends to infinity. For any Reynolds number, including in the limit of infinite Reynolds number, the cascade remains out of equilibrium at any other length-scale which is not a fixed multiple of λ. The larger the length-scale compared to λ the wider the departure from equilibrium. Without equilibrium at length-scales large enough to be independent of viscosity it is not possible to justify the Taylor-Kolmogorov dissipation scaling by a balance between the turbulence dissipation ε and the large-scale interscale energy flux (or rate of energy injection into the cascade at the large scales). Goto & Vassilicos [START_REF] Goto | Unsteady turbulence cascades[END_REF] justified the far-field Taylor-Kolmogorov dissipation scaling on the basis of their concept of balanced non-equilibrium to which we return at the end of this paper.

The features of the non-equilibrium cascade have been up to now studied using highfidelity tools such as experiments and direct numerical simulations. In this paper we investigate whether features of the non-equilibrium cascade can be captured by turbulence closures based on the statistical description of freely decaying isotropic turbulence. To this purpose, the EDQNM model [START_REF] Orszag | Analytical theories of turbulence[END_REF] has been identified. This tool can be used to investigate the statistical features of homogeneous turbulence and in particular the time evolution of freely decaying isotropic turbulence in a far field state. [START_REF] Bos | Spectral imbalance and the normalized dissipation rate of turbulence[END_REF][START_REF] Lesieur | Turbulence in Fluids[END_REF][START_REF] Meldi | Further insights into self-similarity and self-preservation in freely decaying isotropic turbulence[END_REF][START_REF] Sagaut | Homogenous Turbulence Dynamics[END_REF]. One of the main favorable features of this model is that it can be used to investigate freely decaying HIT at Reynolds numbers which are not reachable by high-fidelity tools. The present analysis provides answers to two main questions. The first one is whether the EDQNM is consistent with the absence of interscale equilibrium at all Reynolds numbers in the far field except at scales in the vicinity of the Taylor length scale. The second one is if there is agreement between EDQNM and Lundgren's [START_REF] Thomas S Lundgren | Kolmogorov two-thirds law by matched asymptotic expansion[END_REF] matched asymptotic expansion analysis of the Karman-Howarth equation.

The paper is organized as follows. In section 2 we describe the EDQNM model we use and in section 3 we present our results and the answers to the questions we posed in this introduction. We conclude in section 4.

II. THE EDQNM MODEL

The EDQNM model [START_REF] Orszag | Analytical theories of turbulence[END_REF][START_REF] Lesieur | Turbulence in Fluids[END_REF][START_REF] Sagaut | Homogenous Turbulence Dynamics[END_REF] is a turbulence closure in spectral space. It resolves the numerical discretization of the classical Lin equation:

∂E(k, t) ∂t + 2νk 2 E(k, t) = T (k, t) + F (k, t) (2)
where F (k, t) is a forcing term in the spectral space to be explicitely provided. Present calculations are performed for freely decaying HIT, so F (k, t) = 0. T (k, t) is the nonlinear energy transfer. This term is related to the interscale energy flux via the expression

Π(k, t) = +∞ k T (k, t)dk.
Integrating the Lin equation one obtains an equation where Π(k, t) appears rather than T (k, t) and which is the spectral space equivalent of the Kármán-Howarth equation 1 (see [START_REF] Thomas S Lundgren | Kolmogorov two-thirds law by matched asymptotic expansion[END_REF] and references therein). The EDQNM model is based on two fundamental assumptions applied to the evolution equation of the three-point third-order velocity correlation, which is used to calculate T (k, t):

1. The quasi-normal (QN) approximation is used to express fourth-order moments by a sum of products of second-order moments. The fourth-order cumulants (i.e. deviation from a Gaussian pdf of the velocity derivatives) are represented via a linear damping term which is governed by the eddy damping rate η E .

2.

A Markovianization procedure is used assuming that the relaxation time of the thirdorder correlations is small when compared with the relaxation time of the second-order correlations.

Using these two hypotheses, a closed expression for T (k, t) is obtained:

T (k, t) = p+q=k Θ kpq (xy + z 3 )E(q, t)[E(p, t)pk 2 -E(k, t)p 3 ] dpdq pq (3)
where [k, p, q] represents the spectral system of coordinates and xy + z 3 is a geometric factor determined by the shape of the triangles respecting the condition p + q = k in the spectral space. The term Θ -1 kpq = η E (k, t) + η E (p, t) + η E (q, t) is a spectral time scale which is inversely proportional to the eddy damping rate η E . This term is usually derived by the following model relation [START_REF] Pouquet | Evolution of high reynolds number two-dimensional turbulence[END_REF]:

η E (k, t) = A k 0 p 2 E(p, t) dp + νk 2 . ( 4
)
The free coefficient A ∈ [0.3, 0.5] is usually set to optimize the value for the constant

C K governing E in the scaling region i.e. E(k, t) = C K ε 2/3 k -5/3 [17]. A different model
for the determination of the eddy damping term is the EDQNM-LMFA proposed by Bos & Bertoglio [START_REF] Bos | A single-time two-point closure based on fluid particle displacements[END_REF]. In this model, the eddy damping rate η E is calculated using a model equation for the velocity derivative cross-correlation

F CC , so that η E (k, t) = E(k, t)/F CC (k, t) + νk 2 .
The turbulent kinetic energy K(t), the energy dissipation rate ε(t) and all the other physical quantities are derived via manipulation or integration of E(k, t) and T (k, t). The calculations are performed using an adaptive spectral mesh strategy [START_REF] Meldi | An adaptive numerical method for solving EDQNM equations for the analysis of long-time decay of isotropic turbulence[END_REF] which preserves both the large-scale and the small-scale resolution. This operation is performed by updating the mesh elements so that the conditions k L (t)/k min (t) = s L = const and k η (t)/k max (t) = s η = const are verified. Here, k L (t) = L -1 (t) is the wave number associated with the integral length scale L, k η (t) = η -1 (t) is the wave number associated with the Kolmogorov scale η, k max is the largest resolved mode and k min is the smallest resolved mode. This strategy allows for the complete control of confinement effects. For every calculation s η = 10 -1 is imposed, which implies that the smallest resolved scale has a size of ≈ 0.1η. More details about the value of the constant s L are provided below.

A database of EDQNM calculations has been established in order to quantify the sensitivity of the physical quantities investigated to variations in the set up of the problem as detailed in the following list.

• Initial conditions. Different proposals have been employed to initialize the energy spectrum at time t = 0 of the free decay simulations. The first one is inspired from the functional form found in Pope [START_REF] Pope | Turbulent Flows[END_REF] and Meyers & Meneveau [START_REF] Meyers | A functional form for the energy spectrum parametrizing bottleneck and intermittency effects[END_REF]:

E I (k, t = 0) = C K ε 2/3 k -5/3 f L (kL)f η (kη) (5)
with

f L (kL) = kL [(kL) 1.5 + c L ] 1/1.5 5/3+σ , f η (kη) = exp(-β([(kη) 4 + c 4 η ] 1/4 -c η )) (6)
where C K ≈ 1.6 is the Kolmogorov constant. The dimensionless coefficients in equations 5 -6 have been set to c η = 0.4, β = 5.3; c L has been chosen such as to obtain

L(0) = 1.
The parameter σ, which controls the shape of the energy spectrum at the large scales, is discussed separately analyzed in the next bullet point. A second functional form employed to initialise the energy spectrum at t = 0 is the classical exponential relation

E(k, t = 0) = k σ exp((-σ/2)k 2 ) ( 7
)
which is controled by the parameter σ only.

• Large scale parameter σ. This parameter is strictly concerned with the initialization procedure described above, but the effects of σ on HIT decay do not disappear if sufficient resolution is provided at the large scales [START_REF] Meldi | On non-self-similar regimes in homogeneous isotropic turbulence decay[END_REF][START_REF] Meldi | Turbulence in a box: Quantification of large-scale resolution effects in isotropic turbulence free decay[END_REF]. The values chosen for this parameter are σ = 2, 4 which correspond to the well known cases of Saffman turbulence and Loitsiansky turbulence, respectively.

• EDQNM model. Runs have been perfomed using both the classical version of EDQNM and the EDQNM-LMFA model.

• Resolution at the large scales. The last aspect of the sensitivity analysis deals with confinement i.e. lack of resolution at the large scales. Most of the simulations in the database are run using s L = 500, which ensures that confinement effects are completely excluded. One simulation is instead run for s L = 10, in order to investigate the effects of saturation over the physical quantities investigated. A summary of the features of the calculations included in the database is reported in table I. Every calculation of the database is performed using an initial Reynolds number of Re λ (t = 0) = 10 6 . A transient regime is initially observed, which is governed by the features of the functional form prescribed for t = 0. During this transient, the Reynolds number increases up to Re λ ≈ 2 × 10 6 at t ≈ t 0 , where t 0 = K(t = 0)/ε(t = 0) is the initial turnover time. After this first phase, the statistics progressively lose memory of those initial condition prescribed in the scaling range and in the small scale region. A classical power law decay is then observed [START_REF] Meldi | On non-self-similar regimes in homogeneous isotropic turbulence decay[END_REF][START_REF] Comte-Bellot | The use of a contraction to improve the isotropy of gridgenerated turbulence[END_REF], which is governed by the parameter σ prescribed in the functional form of E. This parameter controls the slope of the energy spectrum at the very large scales.

Data in the form of E(k, t) are sampled in the range t/t 0 ∈ [10 3 , 10 C ], which corresponds to a free decay from Re λ (t/t 0 ≈ 10 3 ) = 10 6 to Re λ (t/t 0 ≈ 10 C ) = 50, C ∈ [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large Reynolds number[END_REF]45]. The time range chosen for sampling varies with initial conditions and with value of σ in particular.

Results are now discussed for simulation 1 of the database. For this simulation the effects of the functional form prescribed at t = 0, which are measured via the rate of convergence towards the expected power-law decay exponent, can be estimated to be ≈ 5% for Re λ = 10 6 and < 0.1% for Re λ = 10 5 . The energy spectrum E, the compensated energy spectrum E/(ε 2/3 k -5/3 ), the derivative of the non-linear energy transfer T and the interscale energy flux Π are shown for Re λ = 10 6 , 10 5 , 10 4 , 10 3 , 10 2 in figure 1. The EDQNM results we obtain comply with previous EDQNM results [START_REF] Bos | Spectral imbalance and the normalized dissipation rate of turbulence[END_REF][START_REF] Meldi | Further insights into self-similarity and self-preservation in freely decaying isotropic turbulence[END_REF] and, plotted as in this figure and without further analysis, may appear to support the picture drawn by the K41 theory [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large Reynolds number[END_REF][START_REF] Kolmogorov | Dissipation of energy in the locally isotropic turbulence[END_REF] for Re λ → +∞. The compensated energy spectrum as well as the non-linear energy transfer collapse remarkably well in the small scale region when plotted against kη.

In addition, one can see the emergence of a scaling range for sufficiently high Re λ (larger than 10 4 according to present data). This range, which can be described by the relation

E(k) = C K ε 2/3 k -5/3 ≈ 1.6ε 2/3 k -5/3
to a close approximation, appears to be characterized by a close-to-zero derivative of the non linear energy transfer and by Π/ε ≈ 1, which one may infer from figure 1(d).

The EDQNM predictions for E and T are used to obtain expressions for second-order and third-order structure functions. These quantities can be calculated exactly via the following integral / differential relations [START_REF] Tchoufag | Spectral approach to finite Reynolds number effects on Kolmogorov's 4/5 law in isotropic turbulence[END_REF][START_REF] Bos | Reynolds number effect on the velocity increment skewness in isotropic turbulence[END_REF]:

S 2 (r, t) = +∞ 0 4 E(k, t) 1 3
f (r, t) = - 1 ε ∂S 2 (r, t) ∂t (10)
F (r, t) = - 3 εr 5 r 0 r ′ 4 ∂S 2 ∂t dr ′ (11)
where S 2 is the second-order moment of the longitudinal velocity increment, S 3 is the third-order moment and f and F are non-stationarity functions [START_REF] Obligado | The non-equilibrium part of the inertial range in decaying homogeneous turbulence[END_REF]. These non-stationarity functions measure the departure from equilibrium scale by scale and appear naturally from the form of the Kármán-Howarth equation (1) given by Landau & Lifschitz [START_REF] Landau | Fluid mechanics: Vol 6 of course of theoretical physics[END_REF] (see also [START_REF] Thomas S Lundgren | Kolmogorov two-thirds law by matched asymptotic expansion[END_REF][START_REF] Danaila | A generalization of Yaglom's equation which accounts for the large-scale forcing in heated decaying turbulence[END_REF]). The -4/5 law S 3 ≈ -4 5 εr might be derived from this equation in a range of scales r where F ≈ 0 (implying equilibrium) and viscous diffusion is negligible. The function F is obtained by a normalised integration of the non-stationarity function f which is directly interpretable in terms of equilibrium as it is small at a given scale r and a given time t if questions which we now address in terms of EDQNM.

III. RESULTS

Results presented in this section are again taken from simulation 1 of the database.

Data are first compared with experimental results in Obligado & Vassilicos [START_REF] Obligado | The non-equilibrium part of the inertial range in decaying homogeneous turbulence[END_REF]. The nonstationarity function f is shown in figure 2 for Re λ = 50, 380. The comparison shows that there is reasonable agreement for low to moderate Reynolds numbers, in particular for small values of r. The EDQNM prediction appears to increase slightly faster towards the asymptotic value f = 4/3 for large r values. However, the difference between EDQNM results and grid turbulence data is small, not larger than ≈ 10% in magnitude on this nonstationarity/non-equilibrium function. This validation supports the view that one might reasonably expect EDQNM data for much higher Reynolds numbers (up to Re λ = 10 6 ) to provide an estimation of the behavior of the physical quantities investigated, such as f and The EDQNM predictions for f and F are shown in figure 3. In the left column, f is plotted against r/L, r/λ and r/η from top to bottom. Two asymptotic behaviors are observed for small and large values of r. (i) One can see that converges towards a constant value of = 4/3 for r > 2L in agreement with basic theory (see [START_REF] Obligado | The non-equilibrium part of the inertial range in decaying homogeneous turbulence[END_REF]). The way the EDQNM profiles reach the asymptote is the same in terms of r/L as the curves collapse significantly well even before reaching the asymptotic value in 3(a). (ii) On the other hand, f ∝ r 2 at the small scales which also agrees with basic theory [START_REF] Taylor | Statistical Theory of Turbulence[END_REF]. This range appears to be fully established for r ≤ 5η approximately (see figure 3(e)). However, the curves do not collapse for any value of r when plotted against r/η. Self-preservation at the small scales is instead obtained when using a r/λ horizontal axis, as shown in figure 3(c). This result confirms that the EDQNM prediction comply very well with power-law turbulence decay and the theoretical

result f ≈ (r/λ) 2 (-K ε 2
dε dt ) in the small scale region [START_REF] Taylor | Statistical Theory of Turbulence[END_REF][START_REF] Obligado | The non-equilibrium part of the inertial range in decaying homogeneous turbulence[END_REF]. The power-law decay implies that (-K ε 2 dε dt ) is independent of time and that the entire time-dependence and/or Reynolds number dependence of f at scales r ≤ 5η is therefore captured by f ≈ (r/λ) 2 , hence the collapse in figure 3(c) which EDQNM is able to reproduce.

If the Reynolds number is sufficiently high, a scaling range where f ∝ r 2/3 emerges between the large scale region and the small scale region. This range can be qualitatively observed for Re λ ≥ 10 3 . In addition, the scaling predicted by EDQNM in this range complies very well with the formula f = -1.5( dε dt /ε)(r/ √ ε) 2/3 which follows from S 2 ≈ 4 (εr) 2/3 . The analysis of the non-stationarity function F , which is reported on the right column of figure 3, allows to draw almost identical conclusions. The only observable difference is the value of the asymptotic limit in the large scale region, which is F = 4/5 in agreement with basic theory [START_REF] Obligado | The non-equilibrium part of the inertial range in decaying homogeneous turbulence[END_REF].

Results are further analyzed to obtain more accurate information about the features of the ranges approximately identified. In figure 4 the exponent n f (r) = d log(f )/d log(r) is reported: it is calculated via polynomial fitting from the EDQNM results. This parameter indicates that a scaling range for which f ∝ r 2/3 is observed for at least one decade only for Re λ ≥ 10 4 . In addition, if one excludes the curve for the lowest Reynolds number Re λ = 10 2 , the inflection point for n f is at r/λ ≈ 5 for all Reynolds numbers.

As already observed, f profiles do not collapse when plotted against r/η and they exhibit self-preservation in the small scale region when plotted against r/λ. However, given that

λ/η ∝ Re 1/2 λ , it follows from f ≈ (r/λ) 2 (-K ε 2 dε dt ) that f Re λ ≈ (r/η) 2 (-K ε 2 dε dt ) where (-K ε 2 dε dt
) is independent of time and Reynolds number for a power-law turbulence decay. The EDQNM prediction complies with this formula, which is shown in figure 5. The function f Re λ exhibits partial self-preservation in the small scale region but also in the intermediate scaling range when plotted against r/η.

We now investigate detailed EDQNM predictions for the interscale energy flux Π and the third order structure function S 3 which are central quantities for characterising the equilibrium or non-equilibrium of the turbulence cascade. The function 1 -Π/ε is plotted in logarithmic scale in figure 6(a). According to Kolmogorov equilibrium, the value of this function should be zero in some inertial range if the Reynolds number is large enough. However, the equilibrium Π/ε = 1 is never observed over an inertial (i.e. viscosity-independent) range for any Reynolds number. The value of r where Π/ε is closest to 1 turns out to be r = 0.42λ for every Re λ in the range investigated. Indeed, the minimum for 1 -Π/ε is observed for r = 0.42λ irrespective of Reynolds number in figure 7(a). As Re λ increases towards infinity, 1 -Π/ε tends to 0 only for values of r which remain a fixed factor or multiple of λ as Re λ varies. In other words the tendency towards equilibrium is only present in the vicinity of the Taylor scale λ in perfect agreement with the prediction of Lundgren [START_REF] Thomas S Lundgren | Kolmogorov two-thirds law by matched asymptotic expansion[END_REF] and Obligado & Vassilicos [START_REF] Obligado | The non-equilibrium part of the inertial range in decaying homogeneous turbulence[END_REF]. As the Taylor scale depends on viscosity we cannot say that the tendency towards equilibrium happens in some inertial range. As figure 6(a) suggests, in a range of scales which are independent of viscosity the degree of non-equilibrium persists with increasing Reynolds number and the closer the length-scale r is to L the further the turbulence is from equilibrium. Results in figure 6(a) indicate that this departure from equilibrium scales as (r/λ) 2/3 for r > λ and as r/λ for r < λ. These formulae can be used to determine span of an empirical range respecting condition Π/ε ≈ 1, down to a prescribed tollerance. In fact, if one fixes a level of tollerance L T tracing a horizontal line in figure 6(a) (let's say for example L T = 10 -2 , which corresponds to an error of 1%), an approximated triangle in the spectral space can be drawn intersecting this line with the curve for a selected Re λ . Using the scalings introduced two sentences above this one, and considering that the vertical position of the lower vertex of the triangle can be approximated by the value min(1 -Π/ε), one can determine the extension of this empirical range ∆ er :

∆ er = 10 σer , σ er = (log 10 (P r ) -log 10 (min(1

-Π/ε))) 1 + 3 2 (12)
The coefficient (1+3/2) on the right is obtained from the observation of slopes r/λ for r < λ and (r/λ) 2/3 for r > λ in figure 6(a). The span of this range increases with Re λ but it also significantly changes with the degree of tollerance imposed for the relation Π/ε ≈ 1.

For high Re λ and relatively low tollerance, one can obtain a significantly large range which will not be exactly centered around λ, as the function 1 -Π/ε exhibits different evolutions moving from λ towards the small scales and the large scales respectively. Such a coarse analysis may erroneously lead one to think that the condition Π/ε ≈ 1 holds in some bit of an inertial range, but the fact remains that Π/ε approaches 1 as Re λ → ∞ only for scales r around λ, i.e. in a range of scales which is not inertial.

In conclusion, the EDQNM model agrees with the matched asymptotic expansion analysis of the Kármán-Howarth equation of Lundgren [START_REF] Thomas S Lundgren | Kolmogorov two-thirds law by matched asymptotic expansion[END_REF] (confirmed by the wind tunnel data analysis of Obligado & Vassilicos [START_REF] Obligado | The non-equilibrium part of the inertial range in decaying homogeneous turbulence[END_REF]) which predicts that there is no inertial range with an approximate equilibrium between Π and ε and that, instead, there is a systematic departure from equilibrium as scales r are considered further and further away from λ. For any fixed r/L value, however small, Π/ε does not tend to 1 as Re λ → ∞, even if r/L is so very small that Π/ε is close to 1.

Similar conclusions can be drawn for S 3 , which is shown in figure 6 shows that this value of r is r = 0.42λ for every Re λ in the range investigated.

IV. SENSITIVITY ANALYSIS

The sensitivity of the results shown in figure 7 is now analyzed by comparing results from different simulations of our database in order to quantify the effects of the variations in the parametric description.

More specifically, the sensitivity on the parameters characterising the large scale features is investigated for min(1 -Π/ε) and for the value of r/λ where minima are observed. This sensitivity analysis is based on results from the simulations 1 to 7 in our database (see table I for more information). The large-scale features encompass the actual parameters governing the large scales, such as σ, but also the two EDQNM models and the functional forms prescribed for E(k, t = 0). The results, which are shown in figure 8, show that variations in σ, in the prescribed form for E(k, t = 0) and in the model for the eddy-damping term in EDQNM do not significantly affect the results. However, some differences may be observed for simulation 2 at Re λ < 400. For this simulation, confinement effects at the large scales are more important than in the other calculations in the database. As Re λ decreases, the reduced scale separation allows the confinement effects to affect min(1 -Π/ε) and the value of r/λ where this minimum is observed.

V. CONCLUSIONS

It is perhaps remarkable how well EDQNM of freely decaying HIT agrees with the results of Lundgren's matched asymptotic expansion approach [START_REF] Thomas S Lundgren | Kolmogorov two-thirds law by matched asymptotic expansion[END_REF] to the Kármán-Howarth equation. Both predict that a Kolmogorov inertial range equilibrium is not obtained. Equilibrium is only achieved in the vicinity of the Taylor length in the limit where the Reynolds number tends to infinity and the Taylor length is not in the inertial range given its explicit dependence on viscosity. As the length-scale r moves away from λ and towards the integral scale L, the turbulence moves progressively further and further away from Kolmogorov equilibrium. The inertial range is in fact characterised by a gradually increasing departure from equilibrium as r increases from λ to L, not by uniform near-equilibrium over much or even a significant part of this range, however large the Reynolds number may be. At a length-scale r taken to be a fixed small fraction (however small) of the integral scale, the ratio of the interscale energy flux to the turbulence dissipation does not tend to 1 a the Reynolds number tends to infinity.

Similar considerations can be extended to the case where r decreases from λ to η except that there is no inertial range at these very small length-scales. However, this last range is significantly shorter and it can be clearly observed only for very high Re λ .

This raises the question of justifying the presence of the Taylor-Kolmogorov scaling of the turbulence dissipation even though there is no equilibrium in the inertial range, particularly at the larger-scale end of it. Goto & Vassilicos [START_REF] Goto | Unsteady turbulence cascades[END_REF] introduced the concept of balanced non-equilibrium for such a justification. They considered the Lin equation integrated over wavenumbers from k to infinity and defined balanced non-equilibrium to mean that the 3 terms in this integrated Lin equation vary together in time for wavenumbers k ∼ 1/L.

Balanced non-equilibrium therefore means that

∂ ∂t +∞ 1/L E(k, t)dk ∼ 2ν +∞ 1/L k 2 E(k, t)dk ∼ Π(1/L, t) (13)
which leads to the Taylor-Kolmogorov dissipation scaling given that 2ν +∞ 1/L k 2 E(k, t)dk ≈ ε and that Π(1/L, t) can be considered to be independent of inlet/initial conditions and viscosity. The validity of balanced non-equilibrium was demonstrated in the context of EDQNM by Meldi & Sagaut [START_REF] Meldi | Further insights into self-similarity and self-preservation in freely decaying isotropic turbulence[END_REF] and therefore explains the fact that EDQNM returns the Taylor-Kolmogorov scaling for the turbulence dissipation as reported in previous publications [START_REF] Bos | Spectral imbalance and the normalized dissipation rate of turbulence[END_REF][START_REF] Meldi | Further insights into self-similarity and self-preservation in freely decaying isotropic turbulence[END_REF][START_REF] Sagaut | Homogenous Turbulence Dynamics[END_REF]. Future potential investigations could employ more sophisticated tools able to take into account anisotropy of the flow, such as specific versions of the EDQNM model reported in the literature [START_REF] Mons | Numerical investigation on the partial return to isotropy of freely decaying homogeneous axisymmetric turbulence[END_REF][START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically averaged descriptors[END_REF][START_REF] Briard | Spectral modelling for passive scalar dynamics in homogeneous anisotropic turbulence[END_REF].

12 T

 12 (k, t) 3(sin(kr) -(kr)cos(kr)) -(kr) 2 sin(kr) (kr) 5 dk (9)

∂S 2 FIG. 1 .

 21 FIG. 1. Energy spectrum and non-linear energy transfer calculated via EDQNM. (a) Normalized energy spectrum E/max(E). (b) Compensated energy spectrum E/(k -5/3 ε 2/3 ). (c) Normalized derivative in Fourier space of the non-linear energy transfer budget term, kT /max(kT ). (d) Normalized inter scale energy flux Π/ε.

5 FIG. 2 .

 52 FIG.2. Non-stationarity function f . EDQNM results are compared with experimental grid turbulence results.

FIG. 3 . 2 FIG. 4 .

 324 FIG. 3. Non-stationarity functions (left column) f and (right column) F represented against (first row) r/L, (second row) r/λ and (third row) r/η.

4 FIG. 5 .

 45 FIG. 5. Non-stationarity function f times the Reynolds number Re λ .

  (b). Here the Kolmogorov equilibrium value S 3 /(εr) = -4/5 for HIT is never reached, except asymptotically as Re λ → +∞ at a value of r/λ which remains fixed as Re λ grows. In fact, Lundgren's (2002) matched asymptotic expansion analysis of the Kármán-Howarth equation predicts that the minimum value of 4/5 + S 3 /(εr) is reached at a value of r which is a fixed multiple of λ and that 4/5 -max(-S 3 /(εr)) ∝ Re -2/3 λ (see Obligado & Vassilicos 2019). Our EDQNM results are in full agreement with these predictions and return a value r = 1.12λ (see figure 7(b)) for the scale r where 4/5 + S 3 /(εr) is minimal and a decay towards 0 of this minimal value which is indeed proportional to Re -2/3 λ as Re λ increases (see figure 7(a)). More precisely, EDQNM results suggest that the minimum can be well approximated by the relation min(1 -Π/ε) ≈ 10Re -2/3 λ . EDQNM results for the function 1 -max(Π/ε) are also shown in figure 7(a) in the range Re λ ∈ [10 2 , 10 5 ]. Reflecting the behaviour of S 3 , 1 -max(Π/ε) ∝ Re -2/3 λ as Re λ → +∞ and the value of r where 1 -(Π/ε) takes its minimum value also scales with λ: figure 7(b)

FIG. 6 .

 6 FIG. 6. (a) 1 -Π/ε and (b) 4/5 + S 3 /(εr) represented in log-log scale.

FIG. 7 .

 7 FIG. 7. (a) Calculated values (for each E and T sampled in Re λ ∈ [10 2 , 10 5 ]) of 1 -max(Π/ε) and 4/5 -max(-S 3 /(εr)). (b) Value of r/λ where minima are observed.

FIG. 8 .

 8 FIG. 8. Sensitivity of (a) calculated values of 1 -max(Π/ε) and (b) value of r/λ where minima are observed, to the large scale features imposed in the EDQNM calculations.

TABLE I .

 I Database of EDQNM calculations used in the present analysis.

	Simulation Nr.	Functional form (t = 0)	σ	EDQNM model	k L (t)/k min (t)
	1	Pope (eq. 5 -6)	2	LMFA	500
	2	Pope	2	LMFA	10
	3	Pope	4	LMFA	500
	4	Pope	4	classical	500
	5	exponential (eq. 7)	2	classical	500
	6	exponential	4	classical	500
	7	exponential	4	LMFA	500


Linqiao Zhang
email: lzhang15@cs.mcgill.ca

Hazel Everett
email: hazel.everett@loria.fr

Sylvain Lazard
email: sylvain.lazard@loria.fr

Sue Whitesides

Towards an Implementation of the 3D Visibility Skeleton

Keywords: F, 2, 2 [Nonnumerical Algorithms and Problems]: Computations on discrete structures Visibility skeleton

In this note we describe the contents of a video illustrating an algorithm for computing the 3D visibility skeleton of a set of disjoint convex polytopes. The video can be found at http://www.cs. mcgill.ca/~lzhang15/video/ with file name socg07visidemo. mov.

INTRODUCTION

The 3D visibility skeleton is a graph whose vertices correspond to the maximal free line segments that are transversal to four edges of at least three distinct polytopes and tangent to those polytopes; its arcs correspond to sets of maximal free line segments that are tangent to three polytopes [START_REF] Durand | The visibility skeleton: a powerful and efficient multi-purpose global visibility tool[END_REF]. The visibility skeleton has been used for visibility computations such as computing shadow boundaries [START_REF] Durand | Fast and accurate hierarchical radiosity using global visibility[END_REF][START_REF] Duguet | Robust epsilon visibility[END_REF].

This video demonstrates a sweep plane algorithm for capturing the vertices of the 3D visibility skeleton of a set of polytopes in 3D [START_REF] Goaoc | Structures de visibilité globales : tailles, calculs et dégénérescences[END_REF][START_REF] Brönnimann | Lines and free line segments tangent to arbitrary three-dimensional convex polyhedra[END_REF].

THE ALGORITHM

The input of the algorithm is a set of k disjoint convex polytopes in general position with n edges in total. The output of the algorithm is the set of O(n 2 k 2 ) vertices of the 3D visibility skeleton of the input polyhedra. The algorithm, which runs in O(n 2 k 2 log n) time, can also be used to compute the arcs of the skeleton.

The algorithm performs a rotational plane sweep around each edge e of each polytope, sweeping from one incident face of that edge to the other incident face. The sweep plane intersects the polytopes in at most k disjoint convex polygons, which change their shape as the sweep plane rotates. During the sweep, the algorithm maintains the 2D visibility skeleton of the intersected polytopes [14]. The 2D visibility skeleton for the convex polygons in the initial sweep plane is computed and then used to determine the initial queue of critical events that will occur during the sweep. The sweep planes at which these critical events occur are called event planes. At a critical event, the 2D visibility skeleton may change its topology and the algorithm updates it, as well as the queue of critical events.

There are three types of critical events. A V-event occurs when the sweep plane encounters a polytope vertex that supports one or more bitangents in the sweep plane. A T-event occurs when two or three bitangents become colinear. An F-event occurs when a bitangent becomes colinear with a face of a polytope. There are O(nk 2 ) events per sweep, and they can be computed and processed in O(nk 2 log n) time in the usual sweep algorithm paradigm.

The vertices of the 3D visibility skeleton are captured during the sweep, as they correspond to the V, T, F-events whose associated bitangents intersect the edge e that the sweep plane is rotating about. After the n sweeps, a description of the arcs of the 3D skeleton can be computed, although the details are not illustrated in the video.

IMPLEMENTATION ISSUES AND TECHNICAL DETAILS

Although the algorithm as described in [START_REF] Brönnimann | Lines and free line segments tangent to arbitrary three-dimensional convex polyhedra[END_REF] works for any set of possibly intersecting convex polytopes in any configuration, the current implementation requires that the polytopes satisfy certain general position assumptions 1 .

The key predicate of the sweep algorithm compares two event planes, to order their occurrences in the sweep. A detailed study of this predicate and its degree is given in [START_REF] Everett | Predicates for line transversals in 3d[END_REF].

The algorithm was implemented in C++ using the CGAL [START_REF]CGAL: Computational Geometry Algorithms Library[END_REF] library. We used the CORE library [START_REF]The CORE library[END_REF] to perform exact comparisons of algebraic numbers. To compute the 2D visibility skeleton we used the CGAL-based package due to Angelier and Pocchiola [START_REF] Angelier | CGAL-based implementation of visibility complexes[END_REF], based on the Greedy Flip Algorithm [START_REF] Angelier | A sum of squares theorem for visibility complexes and applications[END_REF]14].

The graphical output was produced using the Geomview software [START_REF] Geomview | [END_REF] through the interface provided by the CGAL library. We took snapshots of the Geomview window display, while rotating the viewpoint to provide a 3D view of the objects in the display window. Finally, we used iMovie [2] to assemble all the snapshots together into the final video. We used the Audacity [START_REF]Free software distributed under the GNU General Public License (GPL). iMovie[END_REF] software for the audio.

Figure 1

 1 shows one position of the sweep plane, drawn as a faint grid, as it rotates around edge e of polytope C. Polytopes A and B are intersected by the sweep plane. Polytope C lies above the plane, with edge e in the plane. Polytope D lies below the plane. Figure 2 shows the view inside the sweep plane. The polytopes A and B intersect the plane in convex polygons A and B, which support 4 bitangents.Figure 3 shows the 2D

Figure 1 :

 1 Figure 1: One position of the sweep plane.

Figure 2 :

 2 Figure 2: The view inside the sweep plane.

Figure 3 :

 3 Figure 3: The 2D visibility skeleton for Figure 2.

The precise definition of our general position assumptions is straightforward but lengthy. It guarantees, for example, that each critical event corresponds to a unique position of the sweep plane.

ACKNOWLEDGMENTS

We would like to thank Ethan Kim, Christophe Weibel, and Nathan Yu for their help.


P Kervella
email: pierre.kervella@obspm.fr

F Thévenin

M G Petr-Gotzens

The nearby eclipsing stellar system δ Velorum I. Origin of the infrared excess from VISIR and NACO imaging ⋆

Keywords: individual: (HD 74956, δ Vel), Stars: binaries: eclipsing, Methods: observational, Techniques: high angular resolution

 that the IR excess of δ Vel has an interstellar origin, although a cold circumstellar disk could still be present. In addition, we derive the spectral types of the three stars Aa, Ab, and B (respectively A2IV, A4V and F8V), and we estimate the age of the system around 400-500 Myr.

Introduction

The southern star δ Vel (HD 74956, HIP 41913, GJ 321.3, GJ 9278) is the 43 rd brightest star in the visible sky, with m V = 1.96, and is located at a distance of 24.4 pc (π = 40.90 ± 0.38 mas; ESA 1997). Amazingly, δ Vel was discovered only very recently by Otero et al. (2000) to host the brightest eclipsing binary system in the southern sky, one of the very few observable with the unaided eye. Its orbital period of 45 days is exceptionally long in terms of probability of occurrence, especially for such a nearby star. A historical note on the discovery and recent developments in the study of δ Vel can be found in Argyle, Alzner & Horch (2002), together with the orbital parameters of the A-B system (orbital period of ≈ 142 years). In the following, we refer to the eclipsing pair as δ Vel A, and to the fainter, visual component as δ Vel B. The two eclipsing components of A are referred to as Aa and Ab, by order of decreasing brightness. Kellerer et al. (2007) excluded the physical association of the angularly nearby pair of stars sometimes labeled δ Vel C and D with the δ Vel A-B system.

An infrared (IR) excess associated with the δ Vel A-B system had been detected by IRAS (1986). This indicated that δ Vel could belong to the "Vega-like" class of objects, i.e. a main sequence star surrounded by a optically thin debris disk. However, [START_REF] Gáspár | [END_REF] recently discovered a spectacular infrared (IR) bow shock around δ Vel at 24 and 70 µm using Spitzer/MIPS images. This very large structure (≈ 1 ′ ) is explained by these authors as the result of the heating and compression of the interstellar medium by the photons from δ Vel, as the trio moves through the interstellar medium (ISM). From a detailed modeling of the star-ISM interaction, Gáspár et al. conclude that the bow shock contribution is sufficient to explain the observed IR excess of δ Vel without resorting to a circumstellar debris disk. They also note that accretion of interstellar material (λ Bootis phenomenon) could nevertheless take place on δ Vel, although at a slow rate.

In this article, we present the results of our search for circumstellar IR excess in the inner δ Vel A-B system. Sect. 2 is dedicated to the description of our high angular resolution VISIR and NACO observations, whose analysis is presented in Sect. 3. In Sect. 4, we compute the magnitudes of each of the three stars of the system, that are used to derive their physical parameters. We present in Sect. 5 our estimates of the thermal IR excesses of δ Vel Aab and B and a short analysis of the evolutionary state of the system.

Observations and data reduction

VISIR

For the thermal-IR part of our program, we used the VISIR instrument (Lagage et al. 2004), installed at the Cassegrain focus of the Melipal telescope (UT3) of the ESO/Very Large Telescope (Paranal, Chile). VISIR is a mid-IR imager, that also provides a slit spectrometer. VISIR can in principle reach a very high spatial resolving power, thanks to the 8 m diameter of the telescope. However, under standard conditions at Paranal (median seeing of 0.8 ′′ at 0.5 µm), the 8 m telescope is not diffraction limited in the MIR (seeing ≈ 0.4 ′′ vs. 0.3 ′′ diffraction). A few moving speckles and tip-tilt usually degrade the quality of the image (see e.g. [START_REF] Tokovinin | [END_REF]. To overcome this limitation, a specific mode of the instrument, the BURST mode (Doucet et al. 2007a(Doucet et al. , 2007b)), give the possibility to record separately a large number (several tens of thousand) very short exposures (∆t 50 ms), in order to freeze the turbulence. The data processing procedure we developed to reduce the resulting data cubes is described in Kervella & Domiciano de Souza (2007). During the processing, the frames were precisely co-aligned on the position of δ Vel A.

We observed δ Vel and its three calibrators, HD 82668, HD 67582 and HD 80007 in visitor mode at Paranal during the night of 23-24 April 2008. At that time δ Vel was out of ecplipse. The sequence of the observations is presented in Table 1. We adopted standard chopping and nodding offsets of 8 ′′ , with respective periods of 4 and 90 s. The calibrators were selected from the [START_REF] Cohen | [END_REF] catalog of spectrophotometric standards for IR wavelengths, except HD 80007 (for the NeII filter). For each observation of δ Vel and its calibrators, we selected during the data processing ≈ 40% of the total number of frames, rejecting those with the lowest Strehl ratio (estimated from the peak intensity in the frame). For instance, the first two observations #A and B of δ Vel (20 000 selected frames) correspond to a total of 50 000 frames before selection. The number of frames listed in Table 1 corresponds to the result of this selection. The cubes were then averaged to obtain diffraction-limited images in our three filters: PAH1, PAH2, and NeII. that have respective central wavelengths of λ = 8.59, 11.25 and 12.81 µm. The resulting images of δ Vel and the calibrators are presented in Fig. 1. The observation #C of HD 82668 was affected by saturation of the array, and was not used for our anaysis.

NACO

We observed δ Vel at ten epochs in April-May 2008 using the Nasmyth Adaptive Optics System (NAOS, [START_REF] Rousset | Proc. SPIE[END_REF] of the Very Large Telescope (VLT), coupled to the CONICA infrared camera (Lenzen et al. 1998), abbreviated as NACO. These observations were obtained to provide high-precision differential astrometry of the eclipsing system δ Vel A relative to B, and we will analyze them along this line in a forthcoming article.

We selected the smallest available pixel scale of 13.26 ± 0.03 mas/pix (Masciadri et al. 2003), giving a field of view of 13.6 ′′ ×13.6 ′′ . Due to the brightness of δ Vel, we employed a narrow-band filter at a wavelength of 2.166 ± 0.023 µm (hereafter abbreviated as 2.17) together with a neutral density filter (labeled "ND2 short"), with a transmission of about 1.5%.

Table 2 gives the list of the observations. Each of our ten epochs consisted in approximately 50 short exposures grouped over less than 10 minutes, with a detector integration time (DIT) of 0.8 s each to avoid saturation. Four such short exposures were A few images (9 over 510 in total) were affected by the opening of the adaptive optics loop, and were removed from the processing. We obtained nine epochs outside of the eclipses, and one precisely at the phase of the primary eclipse (on 18 May 2008), only 23 minutes before the center of the eclipse. As a remark, the primary eclipse (the deeper in photometry) is when the smaller, cooler star (Ab) passes in front of the bigger, hotter star (Aa), causing an annular eclipse. The secondary eclipse (shallower in photometry) is when the cooler star is totally eclipsed by the hotter one. The epochs of the primary and secondary eclipses of δ Vel A are given by the ephemeris of Otero With E the number of elapsed orbits. To retrieve the phases listed in Table 2, we converted the modified Julian dates of our observations into heliocentric Julian dates using the tools by Dan Bruton2 . The raw images were processed using the Yorick3 and IRAF4 software packages in a standard way. A sample image of δ Vel is presented in Fig. 2.

Image analysis

VISIR

We obtained photometry of δ Vel A in the PAH1 and PAH2 bands using a circular aperture of 0.3 ′′ in radius, centered on We checked that the VISIR images of δ Vel and its calibrators present no detectable diffuse background. This verification was achieved by measuring the average flux at distances larger than 1 ′′ from the central stars, to avoid the inclusion of the highorder Airy rings. Based on this result, we did not apply a sky subtraction to our observations, and as we used the same apertures for the photometry of δ Vel and the flux calibrators, no aperture correction was required either. This non-detection of the diffuse background was expected for two reasons: firstly due to the chopping-nodding observing procedure that removes most of the large scale galactic IR background, and secondly because the sensitivity of our images is insufficient to detect the spatially extended flux from the bow shock in the VISIR bands. The total flux from the model proposed by [START_REF] Gáspár | [END_REF] is at most 3-4 mJy in this wavelength range (Sect. 5), spread over a surface of several arcmin 2 , far beyond our sensitivity limit.

The integrated fluxes are affected by different atmospheric absorption due to the different airmasses compared to the calibrators. To take this into account, we used the empirical formula by [START_REF] Schütz | Proc. ESO Workshop "High Resolution Infrared Spectroscopy in Astronomy[END_REF], that gives the multiplicative correction C(λ, AM) to apply to the raw flux to remove the atmospheric absorption:

C(λ, AM) = 1 + 0.220 - 0.104 3 (λ -8.6 µm) (AM -1) (3)
We corrected separately the different observations of δ Vel and its calibrators.

In order to absolutely calibrate the measurements, the PAH1 and PAH2 fluxes from the calibrator HD 67582 were retrieved from the spectral template of [START_REF] Cohen | [END_REF] at the central wavelength of our filters. For HD 80007 (HIP 45238), we adjusted a model spectrum from [START_REF] Castelli | IAU Symposium 210, Modelling of Stellar Atmospheres[END_REF] to all the photometry available in the literature to retrieve its irradiance in the NeII filter. We obtained 5.62 Jy (1.02 10 -13 W/m 2 /µm) that is within 7% of the value published on the VISIR instrument web page5 (5.26 Jy). The statistical and calibration uncertainties of the measurements were estimated from the dispersion of the different available exposures (e.g. #ABEF for the PAH1 filter).

The uncertainty on the NeII measurement, for which only one image is available, was taken as the dispersion of the calibrated flux for apertures of 150, 300, 600 and 1200 mas. The resulting irradiances are presented in Table 3.

The ratios ρ = f (B)/ f (A) of the fluxes of δ Vel B and δ Vel A in each band can be estimated more accurately than the absolute flux of each star due to the removal of the calibration uncertainty. We obtain the following values in the three filters:

ρ 8.6 µm = 10.3 ± 0.2 %, (4)
ρ 11.25 µm = 11.1 ± 0.4 %, (5)

ρ 12.81 µm = 8.4 ± 2.3 %. (6)

NACO

We first obtained classical aperture photometry of δ Vel A using circular apertures of 1 to 40 pixels in radius (single pixel to 0.53 ′′ radius), for each of our ten epochs. Although the two components of δ Vel are well separated in the NACO images, the diffuse halo from the residuals of the adaptive optics correction of δ Vel A are not negligible at the position of B. For this reason, the photometry of B cannot be measured directly on the images. To remove the diffuse wings of the point spread function (PSF) of A, we subtracted from each point of the image its ring median at the corresponding radius from δ Vel A. With this procedure, we cleanly subtracted the diffuse background from A, mostly without introducing additional noise. Very little residuals from star A are present on the subtracted image, thanks to the good circularity of the PSF (Fig. 4). We then measured aperture photometry using the same 1-40 pixels apertures as for A, centered on B, and computed the flux ratio ρ = f (B)/ f (A) of the pair. We do not need to correct for the airmass, as both A and B are affected by the same atmospheric absorption. We checked in the images that the sky background is negligible due to the short exposures and narrow-band filter.

The choice of the aperture radius is an important factor for the accuracy of the flux ratio measurement. In order to obtain the best match between the two stars, it is essential to have the same integration radius for both stars. Over the 0.6 ′′ distance between the two stars, the variation of the PSF shape and Strehl ratio is negligible, especially for the good seeing conditions of our observations. Fig. 5 shows the variation of the derived flux ratio as a function of the aperture radius. There is a clear inflexion point around 15 pixels in radius (0.2 ′′ ). This corresponds to the radius at which the wings of the PSF of B reach the noise level. The increase of ρ observed for larger apertures corresponds to the inclusion of the residuals from A in the photometry. Fig. 5. Ratio ρ(2.17) = f (B)/ f (A) at 2.17 µm as a function of the aperture radius (in pixels) used for the integration of the flux of A and B. An aperture of 15 pixels was selected for the analysis (square).

To estimate the error bar of the flux ratio, we added quadratically the uncertainties due to the choice of aperture radius over a 10-15 pixel range (considered as systematic, thus not averaging out), and the dispersion of the measurements over the first nine epochs (statistical). The flux ratio is assumed to be constant outside of the eclipses, in order to average the nine measurements (but see also Sect. 4.4). For the epoch of the eclipse, we considered as uncertainty the standard deviation of the ratios obtained at the other epochs.

From this procedure, we obtain outside of the eclipse: ρ 2.17 µm = 9.66 ± 0.05%

corresponding to a magnitude difference of ∆m = 2.537 ± 0.005.

During the primary eclipse of δ Vel A, the flux ratio becomes:

ρ 2.17 µm (Eclipse) = 14.50 ± 0.12% (8)
giving ∆m(Eclipse) = 2.097 ± 0.009. The eclipse depth at 2.17 µm is therefore d(∆m) = 0.440 ± 0.011 mag. As a remark, we did not deconvolve our NACO images, as we did not obtain a PSF calibrator for this purpose, and we leave the astrometric analysis of these images for a future article.

Photometric analysis

Visible magnitudes

If not specified otherwise, we assume here that the photometric measurements of the δ Vel system available in the literature always include the A and B components, as their angular separation has been smaller than 2 ′′ for at least five decades. Bessell (1990) gives: Fig. 6. Measured value of f (B)/ f (A) at 2.17 µm for our 10 measurement epochs. Note the primary eclipse occuring at the last epoch.

m V (Aab + B) = 1.
with an associated uncertainty of ±0.02 mag.

The photometric observations from Hipparcos (ESA 1997) give the differential magnitude of Aab (combined) and B in the H p band, that are sufficiently close to V to consider the values identical within their error bars:

∆m H p (Aab/B) ≈ ∆m V (Aab/B) = 3.58 ± 0.07 (10)
Considering Bessel's data, this translates immediately into the following magnitude for Aab and B:

m V (Aab) = 2.00 ± 0.02 (11) m V (B) = 5.54 ± 0.08 (12)
In addition, the combined visible photometry available in the literature (see the compilation by Otero et al. 2000) gives us the combined Aab+B magnitudes, the depth of the primary (I) and secondary (II) eclipses: As shown by the curve presented by Otero et al. (2000), the secondary eclipse of δ Vel Aab is total, and we can thus derive immediately the visible magnitudes of Aa and Ab:

m V (Aa) = 2.33 ± 0.03 (17) m V (Ab) = 3.44 ± 0.03 (18)

Infrared magnitudes

In addition, our NACO observations (Sect. 3.2) give us the relative fluxes of Aab and B in the narrow band filter band at 2.17 µm with a FWHM of 0.023 µm: ρ 2.17 µm = 9.66 ± 0.05%. As this wavelength corresponds to an absorption line of hydrogen that is present in the spectra of hot stars, the conversion of the measured narrow-band ratios to standard K band ratios requires taking the position of the quasi-monochromatic wavelengths within the bands and the shape of the observed spectra into account. We thus proceeded in two steps to obtain the conversion between the narrow-band flux ratio and the broadband flux ratio. Starting from the narrow-band fluxes measured on δ Vel's components, we computed the radius and effective temperature of each star assuming (as a first approximation) that the K flux ratios are identical to the 2.17 µm ratios, using the method descussed in Sect. 4.5. This gave us a first estimation of the spectral types of the stars. We then used the Pickles (1998)6 reference spectra corresponding to their spectral types (A2IV, A4V and F8V) to recover the corresponding broadband flux ratio. The K band standard filter profile was taken from Bessell & Brett (1988). We observed that the true ratios in the K band are very close to the ratios in the 2.17 µm filter, within 0.01 mag for all three stars, which is small compared to the other uncertainties. We thus neglect this difference in the following analysis.

The COBE/DIRBE instrument measured F K = 121.7 ± 11.4 Jy at 2.2 µm (Smith et al. 2004), corresponding to m K = 1.77 ± 0.10, and the 2MASS catalog [START_REF] Skrutskie | [END_REF] gives an apparent K band magnitude of m K = 1.72 ± 0.26 for Aab+B, equivalent to F K = 127 ± 34 Jy. Averaging these two measurements, we obtain m K = 1.76 ± 0.10 and consequently:

m K (Aab) = 1.86 ± 0.09 (19) m K (B) = 4.40 ± 0.09 (20)

Photometry during the primary eclipse

The primary eclipse of Aab in the K band observed with NACO can give us the difference in surface brightness of Ab and Aa. As this eclipse is not total, but annular, we cannot derive unambiguously the brightness ratio of the two stars as it depends on their angular diameters. In this paragraph, we consider that the eclipse is perfectly centered, thus neglecting a possible differential effect of the limb darkening (LD) of the two stars. As we observed very close to the eclipse minimum and in the K band, where the LD is small (the LD/UD correction is less than 3%), this appears as a reasonable assumption. We also neglect the possible distortion of the photospheres of Aa and Ab due to a possible fast rotation.

During the primary eclipse, the disk of Ab "replaces" part of the disk of Aa with a lower effective temperature (hence lower surface brightness) "patch". The eclipse depth in the K band is 0.44 ± 0.01 mag (Sect. 3.2). This can be compared with the depth of the primary eclipse in the visible: ∆m V = 0.53 ± 0.01. The "color" of the primary eclipse is therefore ∆m V -∆m K = 0.09±0.02 mag. This corresponds to the "differential difference" in color between the two stars between the visible and IR wavelengths. This leads to a magnitude difference of:

m K (Ab) -m K (Aa) = 1.02 ± 0.04 (21)
corresponding to a flux ratio f (Ab)/ f (Aa) = 39.2 ± 1.4 % and:

m K (Aa) = 2.22 ± 0.09 (22) m K (Ab) = 3.24 ± 0.09 (23)

Emergence of Ab from the secondary eclipse

It is interesting to notice that our 2008-04-23 NACO observation has been obtained only about 30 minutes after the end of the secondary (total) eclipse during which the δ Vel Aa component covered Ab. The apparent disks of the two stars are therefore "touching" each other. Our measurement shows a slight variation of the photometric ratio ρ = f (B)/ f (A), as we measure 9.90 ± 0.12%, compared to an average value over the eight other out-of-eclipse epochs of 9.63±0.05%. This means that the masking of part of the close environment of δ Vel Ab by Aa results in a slight loss of total flux of the Aab system. Although of marginal statistical significance (2.1 σ), this measurement is the maximum value of the out-of-eclipse ratios we obtained. It could be due to the presence of a circumstellar flux contributor of unknown nature in the immediate vicinity of δ Vel Ab, accounting for 0.27 ± 0.13% of the total flux of Aab. If Ab is a fast rotator, and its rotation axis is perpendicular to the orbital plane, one hypothesis could be that it is the signature of a circumstellar disk made of gas (Ae episode).

Angular diameters, radii and effective temperatures

From our spatially resolved photometry of the three stars of δ Vel (Aa, Ab and B), we can deduce their limb-darkened angular diameters θ LD , photospheric radii R and effective temperatures T eff using the surface brightness-color (SBC) relations calibrated by Kervella et al. (2004) in (V, V -K) and the Hipparcos parallax. The results are summarized in Table 4. This approach allows us to rely only on empirical SBC relations (calibrated by interferometric angular diameter measurements) to compute the parameters of the three stars. Another approach to determine these parameters would be to adjust directly synthetic spectra to the photometry. While the comparison with such model spectra is necessary to determine the IR excess (Sect. 5), it has at this stage the disadvantage of making the derived parameters dependent on the choice of one particular numerical model library. For this reason, we prefer to use the empirical SBC relations. In any case, the good agreement between the model spectra and the measured photometry in the visible and near-IR domains shown in the bottom parts of Fig. 7 and 8 indicates that both methods give consistent results.

The parameters of δ Vel are in good agreement with the values deduced for the combined system A+B from spectroscopy by Allende Prieto & Lambert (1999): T eff = 8700 K, R = 3.9 R ⊙ . The larger radius is naturally explained by the combination of the flux from Aa and Ab (and to a lesser extent of B), while the effective temperature is a weighted average of the values of the three stars. They also find an effective gravity of log g = 3.66 and a mass of M = 2.51 M ⊙ . The reliability of these values is uncertain due to the effect of the orbital motion of Aab on the spectral line broadening. However, we do not confirm the large radii derived by Kellerer et al. (2007) for Aa and Ab (6.0 and 3.3 R ⊙ ) from interferometric observations of the trio.

Discussion

The IR excess detected around δ Vel by IRAS was recently attributed by [START_REF] Gáspár | [END_REF] to interstellar dust heated by δ Vel. In this Section, we examine the possibility that one of the components of δ Vel A or B presents an IR excess of circumstellar origin. For this purpose, we compare the fluxes extracted from our spatially resolved VISIR observations to the expected photospheric flux from the Aab and B components.

To retrieve the synthetic stellar spectra of the three stars from the library assembled by [START_REF] Castelli | IAU Symposium 210, Modelling of Stellar Atmospheres[END_REF], hereafter CK03, we used the parameters listed in Table 4. We read the tables for log g values of 4.0 for δ Vel Aa and Ab, and 4.5 for B, and solar metallicities. The sensitivity of the models to small variations of these two parameters is minimal. We would like to emphasize that we did not fit these spectra to the available photometry, but we simply multiplied the CK03 spectra by the squared limb-darkened angular diameter of each star (Table 4).

From these spectra, we derived the expected thermal IR flux from the stellar photospheres of Aa+Ab and B in the PAH1, PAH2 and NeII filters corresponding to our VISIR observations. We can also add the CK03 spectra of the three stars in order to compare the predicted flux of the trio to the existing non spatially resolved measurements of A+B from the literature. The resulting spectrum is shown in Fig. 8. The photometric measurements are taken from [START_REF] Thompson | Catalogue of stellar ultraviolet fluxes (TD1)[END_REF] UV fluxes from TD1), Morel et al. (1978;U BVRI), [START_REF] Skrutskie | [END_REF]2MASS JHK), Smith et al. (2004; COBE/DIRBE, four wavelengths from 1.25 to 4.5 µm), VISIR (8.6, 11.25 and 12.81 µm, this work), Jayawardhana et al. (2001;N and 18.2 µm Table 4. Properties of δ Vel Aa, Ab and B derived from their V and K band magnitudes using the surface brightness-color relations from Kervella et al. (2004). The bolometric corrections for the V band were taken from Bessell et al. (1998). The approximate masses of the three stars are discussed in Sect. 5. As shown in Fig. 7, Fig. 8 and Table 5, the comparison with our VISIR spatially resolved photometry of δ Vel A and B does not show a significant circumstellar IR excess in the VISIR bands compared to the expected photospheric flux of each of the two components, at a few percent level. The excess that is visible longwards of 10 µm in Fig. 8 comes from the flux contribution of the ISM bow shock at large angular distances that has been observed by [START_REF] Gáspár | [END_REF]. We note that the IRAS and Jayawardhana et al. ( 2001) photometry shows an excess compared to their model between 12 and 25 µm. This may be a consequence of the fact that the flux contribution from the ISM is (at least partly) included in the large aperture measurements (IRAS in particular), but not in our narrow aperture photometry. We notice a relatively large ultraviolet excess on δ Vel in Fig. 8. A fast rotating star usually presents a UV excess due to its overheated polar caps (von Zeipel effect, see e.g. Aufdenberg et al. 2006). A discussion on the rotational velocity of the components of δ Vel A can be found in Argyle, Alzner & Horch (2002). As a remark, Royer, Zorec & Gomez (2007) found a rotational velocity of v sin i = 150 km/s, but this was computed for the combined spectrum of the three stars of δ Vel so a bias may be present.

The physical parameters derived for the three stars (Table 4) indicate that δ Vel Aa already evolved away from the main sequence. This is confirmed by the evolutionary tracks from Girardi et al. (2000), shown in Fig. 9, that give approximate masses of 2.5, 2.0 and 1.3 M ⊙ respectively for Aa, Ab and B, and an age of approximately 400 -500 Myrs for the system. This age compares well with the age of 330-390 Myr proposed by Rieke et al. (2005), while the masses are in relatively good agreement with Argyle, Alzner & Horch (2002), who determined a total dynamical mass of the δ Vel AB system of 5.7 +1.3 -1.1 M ⊙ . We postpone a more detailed discussion of the ages and evolutionary status of the δ Vel stars to a forthcoming paper.

Conclusion

Our VISIR photometry of δ Vel A and B does not show the presence of an excess of circumstellar origin in the thermal IR domain (8-13 µm) at a level of a few percents. This result indicates that these stars do not host a warm debris disk with a typical temperature around 200-300 K. This supports the conclusions of [START_REF] Gáspár | [END_REF] who attribute the IR excess detected with Spitzer to the surrounding interstellar medium material. The ab-sence of a warm circumstellar disk in the inner δ Vel system may be a consequence of the gravitational interactions between the three stars. They could have caused the dispersion of the residual circumstellar material from which they stars formed. However, the possibility still exists that a cold debris disk is present, as in the case of Vega (Su et al. 2005) and Fomalhaut (Stapelfeldt et al. 2004), but at large distances from the stars. In order to test this scenario, observations of δ Vel in the far IR or submillimetric domain similar for instance to those of Fomalhaut by Holland et al. (2003) would be necessary.

From the flux ratios in the V and K bands, we could derive the physical properties of the three stars Aa, Ab and B. δ Vel appears as a moderately evolved system, with the Aa component currently leaving the main sequence. As suggested by the v sin i value of Royer et al. (2007), the ultraviolet excess of δ Vel and our putative observation of near-IR circumstellar emission close to Ab may indicate that δ Vel Ab (and/or Aa) could be a fast rotator, turning this eclipsing system into a promising object to map the polar brightening associated to the von Zeipel effect (von Zeipel 1924).

Fig. 2 .

 2 Fig. 2. Median NACO image of the 51 frames of δ Vel obtained on 2008-04-01.

Fig. 3 .

 3 Fig. 3. Apertures (circles) used for the VISIR photometry of δ Vel A (left) and B (right) shown respectively on the PAH1 image #B and the subtracted image #B -#D (see text for details).

Fig. 4 .

 4 Fig. 4. Aperture (15 pixels = 0.2 ′′ in radius) used for the integration of the flux of δ Vel A (left) and B (right), the latter on the ring-median subtracted NACO image (see text for details).

  ∆m V (I, Aab + B) = 0.51 ± 0.01 (13) ∆m V (II, Aab + B) = 0.32 ± 0.01 (14) By removing the flux contribution of B computed precedently, we obtain the following eclipse depths for Aab only: ∆m V (I, Aab) = 0.53 ± 0.01 (15) ∆m V (II, Aab) = 0.33 ± 0.01 (16)

Fig. 7 .

 7 Fig. 7. Spatially resolved photometry of δ Vel Aab and B compared to Castelli & Kurucz (2003) models of Aa+Ab and B, considering the stellar properties listed in Table 4 (top) and the residuals compared the observations (bottom).

Fig. 8 .

 8 Fig.8. Total model flux of δ Vel Aab+B (solid curve), together with the available photometry from the literature (squares) and the VISIR photometry (filled circles). The dashed curve is the flux contribution of the IR excess model proposed by[START_REF] Gáspár | [END_REF]. The bottom plot shows the relative photometry with respect to the combined photospheric model of the three stars.

Fig. 9 .

 9 Fig. 9. Positions of δ Vel Aa, Ab and B in the HR diagram, with the isochrones and isomass curves from Girardi et al. (2000) overplotted.

Table 1 .

 1 Log of the observations of δ Vel and its calibrators, HD 82668, HD 67582 and HD 80007 with VISIR. MJD is the modified Julian date of the middle of the exposures on the target, minus 54 580. The Detector Integration Time (DIT) is given in milliseconds for one BURST image. θ is the seeing in the visible (λ = 0.5 µm) as measured by the observatory DIMM sensor, in arcseconds. N exp. is the number of selected and averaged exposures. The airmass (AM) is the average airmass of the observation.

	#	MJD *	Star	Filter DIT N exp. θ ( ′′ ) AM
	A 0.0134	δ Vel	PAH1	20	10 000 0.86 1.18
	B 0.0231	δ Vel	PAH1	20	10 000 1.00 1.20
	C 0.0374 HD 82668 PAH1	20	5 000	1.10 1.20
	D 0.0603 HD 67582 PAH1	16	10 000 0.84 1.30
	E 0.0774	δ Vel	PAH1	16	8 000	0.70 1.33
	F	0.0874	δ Vel	PAH1	16	8 000	0.67 1.37
	G 0.0977	δ Vel	PAH2	8	10 000 0.60 1.42
	H 0.1093	δ Vel	PAH2	8	12 000 0.53 1.48
	I	0.1206	δ Vel	PAH2	8	10 000 0.80 1.54
	J	0.1363 HD 67582 PAH2	8	10 000 1.15 1.92
	K 0.1480 HD 67582 PAH2	8	10 000 1.50 2.09
	L 0.1632	δ Vel	NeII	16	10 000 1.29 1.91
	M 0.1860 HD 80007	NeII	16	10 000 1.12 1.46
	N 0.1979 HD 80007	NeII	16	10 000 1.00 1.52

Fig. 1. Top row: VISIR average images of δ Vel in the PAH1

	Bottom row: Calibrator images in the PAH1 (HD 67582), PAH2
	(HD 67582) and NeII (HD 80007) bands.
	co-averaged providing 3.2 s integrated exposure time per image.

(#A, B, E and F in Table

1

), PAH2 (#G, H, I) and NeII (#L) bands.

Table 2 .

 2 Log of the observations of δ Vel with NACO in the NB2.17 filter. MJD * is the average modified Julian date minus 54 000, HJD

	AM

* the heliocentric Julian date minus 2 454 000 and φ is the phase of the eclipsing binary (see text). N is the number of frames. As in Table

1

, θ is the seeing in the visible and AM is the airmass. The listed figures are median values over the observations. # MJD * HJD * φ(Min I) N exp. θ ( ′′ )

Table 3 .

 3 Measured thermal IR irradiances of δ Vel A and B in S.I. units and in Jy. The relative uncertainty σ on the irradiance values is listed in the last column. Star Filter λ (µm) 10 -14 W/m 2 /µm

	Jy	σ

Table 5 .

 5 Observed fluxes and excesses of δ Vel compared to the Castelli & Kurucz (2003) models. Nσ is the measured excess expressed in number of times the measurement uncertainty.

	Band	λ (µm) F (10 -13 W/m 2 /µm) F/F mod (%) Nσ
	Aab+B
	PAH1	8.59	4.07 ± 0.27	99 ± 7	-0.2
	PAH2	11.25	1.38 ± 0.12	98 ± 9	-0.3
	NeII	12.81	7.26 ± 0.08	86 ± 11	-1.3
	Aab
	PAH1	8.59	3.69 ± 0.27	98 ± 7	-0.2
	PAH2	11.25	1.24 ± 0.12	97 ± 10	-0.3
	NeII	12.81	0.67 ± 0.07	88 ± 11	-1.1
	B
	PAH1	8.59	0.382 ± 0.028	104 ± 7	0.5
	PAH2	11.25	0.138 ± 0.014	109 ± 10	0.9
	NeII	12.81	0.056 ± 0.019	74 ± 34	-0.7

http://ar.geocities.com/varsao/delta Velorum.htm

http://www.physics.sfasu.edu/astro/javascript/hjd.html

http://yorick.sourceforge.net/

IRAF is distributed by the NOAO, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation.

http://www.eso.org/sci/facilities/paranal/instruments/visir/

http://www.ifa.hawaii.edu/users/pickles/AJP/hilib.html

Acknowledgements. We are grateful to the ESO observing team at Paranal for the perfect execution of our NACO observations of δ Vel, in particular for the time-critical eclipse on May 18th, 2008. We also warmly thank Sebastian Otero (Asociación Cielo Sur, Argentina), Brian Fraser (Sunninghill Observatory, South Africa) and Jaime Garc ía (Instituto Copérnico, Argentina) for providing the visible photometry of δ Vel. This work received the support of PHASE, the high angular resolution partnership between ONERA, Observatoire de Paris, CNRS and University Denis Diderot Paris 7. This research took advantage of the SIMBAD and VIZIER databases at the CDS, Strasbourg (France), and NASA's Astrophysics Data System Bibliographic Services.