Lutz Lesshafft
email: lesshafft@ladhyx.polytechnique.fr

Artificial eigenmodes in truncated flow domains

published or not. The documents may come

Artificial eigenmodes in truncated flow domains

Introduction

Linear instability analysis of open flows today is commonly carried out in a so-called 'global' framework, where at least two non-homogeneous spatial directions of a steady base state are numerically resolved. In contrast to 'local' analysis, where the base state is assumed to be parallel, and unbounded in the flow direction, a global discretisation of an open flow problem in a truncated numerical domain necessitates the formulation of artificial streamwise boundary conditions for flow perturbations. The question then arises in how far such boundary conditions influence the instability dynamics of the truncated flow system.

This paper investigates the effect of spurious pressure feedback, due to domain truncation, on the eigenmode spectrum of incompressible open flow problems. The investigation is motivated by observations made in recent linear instability studies of jet flows [START_REF] Garnaud | Modal and transient dynamics of jet flows[END_REF][START_REF] Coenen | Global instability of low-density jets[END_REF], where a prominent family of eigenmodes (black symbols in Fig. 1), referred to as the 'arc branch' from here on, was found to present features that suggest a resonance between the inflow and outflow boundaries. Such branches are in fact ubiquitous in many, if not all, global spectra of truncated open shear flows found in the literature: boundary layers [START_REF] Ehrenstein | On two-dimensional temporal modes in spatially evolving open flows: the flat-plate boundary layer[END_REF][START_REF] Ehrenstein | Two-dimensional global low-frequency oscillations in a separating boundary-layer flow[END_REF][START_REF] Åkervik | Global two-dimensional stability measures of the flat plate boundary-layer flow[END_REF], cylinder wakes [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF][START_REF] Marquet | Sensitivity analysis and passive control of cylinder flow[END_REF], jets [START_REF] Nichols | Global modes and transient response of a cold supersonic jet[END_REF][START_REF] Garnaud | Modal and transient dynamics of jet flows[END_REF], plumes [START_REF] Chakravarthy | Local and global dynamics of buoyant jets and plumes[END_REF], three-dimensional boundary layers with roughness elements [START_REF] Loiseau | Investigation of the roughness-induced transition: global stability analyses and direct numerical simulations[END_REF][START_REF] Kurz | Mechanisms of flow tripping by discrete roughness elements in a swept-wing boundary layer[END_REF] -all these and many others present similar characteristic branches of eigenvalues that are often described as being highly dependent on the type or position of outflow boundary conditions. Typically, no convergence with respect to the length of the numerical box can be attained for such modes.

Ehrenstein & Gallaire [START_REF] Ehrenstein | On two-dimensional temporal modes in spatially evolving open flows: the flat-plate boundary layer[END_REF] remark on the resemblance between arc-branch-type global eigenmode structures obtained for a flat-plate boundary layer and spatial modes as found in a local analysis. Åkervik et al. [START_REF] Åkervik | Global two-dimensional stability measures of the flat plate boundary-layer flow[END_REF] document the dependence of boundary layer eigenvalues on the type of outflow boundary conditions. Heaton et al. [START_REF] Heaton | Global linear stability of the non-parallel Batchelor vortex[END_REF] characterise arc branch modes in the spectrum of a Batchelor vortex as artifacts, speculating that these arise from the limited precision of their numerical scheme, not from domain truncation. Cerqueira & Sipp [START_REF] Cerqueira | Eigenvalue sensitivity, singular values and discrete frequency selection mechanism in noise amplifiers: the case of flow induced by radial wall injection[END_REF] demonstrate that such precision errors indeed lead to spurious 'quasi-eigenmodes', but that these differ from the arc branch. In their analysis, quasi-eigenmodes appear in regions of the complex frequency plane where pseudospectrum ε values [START_REF] Trefethen | Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators[END_REF] are very small, below approximately 10 -12 . These modes, also visible in the lower left corner of Fig. 1, are very sensitive to the numerical scheme, to mesh refinement and to the eigenvalue shift parameter. Arc branch modes are found to be robust with respect to those details [START_REF] Cerqueira | Eigenvalue sensitivity, singular values and discrete frequency selection mechanism in noise amplifiers: the case of flow induced by radial wall injection[END_REF], but strongly dependent on the numerical domain length. Coenen et al. [START_REF] Coenen | Global instability of low-density jets[END_REF] show that arc branch eigenfunctions of a jet are characterised by an integer number of wavelengths between the inflow and outflow boundaries, and they suggest an analogy with acoustic modes in a pipe of finite length. This analogy implies that arc branch modes are the result of unwanted resonance between the numerical boundaries, potentially leading to a spurious instability of the numerical system. The present study expands on all these observations, and it aims at a detailed characterisation of unphysical resonance due to imperfect boundary conditions in open shear flow calculations.

The possibility that global instability in truncated systems may be brought about by spurious pressure feedback from boundary conditions was probably first described by Buell & Huerre [START_REF] Buell | Inflow/outflow boundary conditions and global dynamics of spatial mixing layers[END_REF]. In direct numerical simulations of perturbations in a mixing layer, unstable perturbation growth was observed at long times, in a configuration that was only convectively unstable in a local sense. Such behaviour is inconsistent with the interpretation of local convective instability, and it was demonstrated to be caused by unphysical pressure perturbations emanating from the outflow boundary, which in turn provoked the formation of vortical perturbations at the inflow boundary. Chomaz [5, section 3.2.2] interprets pressure feedback as a non-local operator variation, arguing that strong non-normality in the spectrum of convection-dominated flows is likely to induce a high sensitivity of global eigenmodes with respect to such feedback.

While the spurious generation of acoustic pressure waves from artificial boundary conditions is an important and much-discussed problem in compressible flow simulations, especially those that aim to accurately capture the acoustic radiation from shear flows [START_REF] Colonius | Modeling artificial boundary conditions for compressible flow[END_REF], the question of how such artifacts may affect the global stability behaviour remains largely unexplored, both in compressible and in incompressible configurations.

The problem is approached here in the following manner. Section 2 presents a global instability analysis of the Ginzburg-Landau equation with explicit feedback from a downstream sensor to an upstream actuator, in order to examine the effect of such feedback on the eigenvalue spectrum in a controlled setting. This model study fully describes the suspected mechanism behind arc branch modes. In Sect. 3 the analysis is extended to a parallel jet flow in a finite-size numerical domain. Global eigenmodes are projected onto their local counterparts, the global pressure field is examined, and spurious feedback effects are analysed. The question how spurious feedback may be reduced in a practical manner is addressed in Sect. 4, for the example of a spatially developing jet.

A Ginzburg-Landau model problem

The hypothesis put forward by Coenen et al. [START_REF] Coenen | Global instability of low-density jets[END_REF], that the arc branch may be the manifestation of a non-physical upstream scattering of perturbations from the outflow boundary, is first investigated with the help of a simpli-fied model. The linear Ginzburg-Landau equation is written as

∂ t ψ = -U∂ x ψ + µ(x)ψ + γ∂ xx ψ + f (x, ψ) . (1)
The complex scalar variable ψ(x,t) is a function of time t and of one single spatial coordinate x. Constant parameters U = 6 and γ = 1 -i are chosen ad hoc, whereas the coefficient µ varies linearly in x as

µ(x) = U 2 8 1 - x 20 . (2)
With this particular variation of µ, the system is marginally absolutely unstable at x = 0, convectively unstable for 0 < x < 20, and stable for x > 20, all in a local sense. Ginzburg-Landau systems of this form, with linearly decreasing µ(x) and without feedback, f ≡ 0, have been extensively used to model the global instability behaviour of spatially developing flows in semi-infinite domains [see for instance 9]. The instability characteristics mimic, in a very simple and qualitative manner, those of a spreading jet. However, the system (1) is not intended here to predict or reproduce the dynamics of any specific flow.

A forcing term f (x, ψ) is added in (1) in order to provide an explicit closed-loop forcing between the upstream and the downstream end of the flow domain. Some aspects of closed-loop forcing in the Ginzburg-Landau equation are discussed by Chomaz [START_REF] Chomaz | Global instabilities in spatially developing flows: Non-normality and nonlinearity[END_REF]; in the present context, it is used to model a suspected spurious feedback in global shear flow computations. Taken to be of the form

f (x, ψ) = C exp - (x -x a ) 2 0.1 2 ψ(x s ), (3)
a forcing proportional to ψ(x s ) is applied in a close vicinity of x a , such that a feedback loop is established between a (downstream) sensor location x s and an (upstream) actuator location x a . The complex coefficient C governs the amplitude and phase of the feedback, and the Gaussian spreading around x a is introduced for reasons of numerical resolution. Equation ( 1) is discretised on an interval 0 x 40, with a step size ∆ x = 0.1, using an upwind-biased seven-point finite difference stencil for the spatial derivatives. A homogeneous Dirichlet boundary condition for ψ is imposed at the upstream boundary, consistent with typical jet conditions. Actuator and sensor locations are chosen close to the boundaries, at x a = 1 and x s = 39, where spatial derivatives are well resolved.

Temporal eigenmodes of (1) are sought in the form ψ(x,t) = ψ(x) e -iωt . The eigenvalues of the system without feedback, C = 0, are known analytically to be

ω n = i U 2 8 - U 2 4γ + γ 1 3 U 4 3 160 2 3 ζ n , (4)
where ζ n is the n th root of the Airy function [START_REF] Chomaz | Bifurcations to local and global modes in spatially developing flows[END_REF]. These values are represented in Fig. 2 as circles. Each frame 2(a-d) also displays the numerically computed eigenvalues of systems with feedback, shown as bullet symbols, for different non-zero values of the coefficient C. Already with very low-level feedback, C = 10 -10 , the spectrum is clearly affected: only the leading three eigenvalues of the unforced system are recovered, and the lower part of the spectrum is replaced with two new stable branches. These branches, named feedback branches in the following, move upward in the complex ω plane as the feedback coefficient is increased, masking more and more of the original unforced eigenvalues. Note that those affected original eigenvalues are not merely altered by the feedback, but they rather disappear abruptly from the spectrum as they fall below the new branches.

The spectra obtained with C = 10 -6 and 10 -2 resemble those of pure helium jets, as shown in figure 6 of [START_REF] Coenen | Global instability of low-density jets[END_REF], and those of cylinder wakes obtained by Marquet et al. [20, their figure 16]. The least stable original eigenvalue lies just above the feedback branches, apparently unaffected. In the strongest feedback case, C = 10 2 , the feedback branches have merged into one, overarching the entire spectrum of the feedback-free system, which has altogether disappeared. This picture (Fig. 2d) resembles the spectrum of the slowly developing jet of Coenen et al. [START_REF] Coenen | Global instability of low-density jets[END_REF], reproduced in our Fig. 1.

The similarity between the feedback branches in the present model and the arc branch in the jet spectrum (Fig. 1) is also manifest in the pseudospectra. If feedback is thought of as a variation of the operator [START_REF] Chomaz | Global instabilities in spatially developing flows: Non-normality and nonlinearity[END_REF], it should be possible to relate feedback-induced eigenvalues to the pseudospectrum. The pseudospectrum is defined here by the spectral norm of the resolvent operator, at any complex frequency ω [START_REF] Trefethen | Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators[END_REF]. Pseudospectrum ε-contours of the Ginzburg-Landau equation are compared between the unforced setting with C = 0, Fig. 3(a), and the forced setting with C = 10 -6 , Fig. 3(b). It is seen that the feedback branches align closely with an isocontour of the feedback-free pseudospectrum. This criterion also applies to all other cases displayed in Fig. 2, with different values of C, and it is fully consistent with observations made in three flow configurations [START_REF] Heaton | Global linear stability of the non-parallel Batchelor vortex[END_REF][START_REF] Cerqueira | Eigenvalue sensitivity, singular values and discrete frequency selection mechanism in noise amplifiers: the case of flow induced by radial wall injection[END_REF][START_REF] Coenen | Global instability of low-density jets[END_REF]. Furthermore, the pseudospectrum of the system with feedback is identical to that of the system without feedback above the feedback branches, whereas the pseudospectrum is flattened, nearly constant, below the feedback branches. The same behaviour is found in the jet pseudospectrum shown in Fig. 1. Finally, the discrete distribution of feedback modes along the branch is investigated. The strong feedback case C = 10 2 is considered for illustration. Eigenfunctions of successive feedback modes are presented in figure 4, analogous to the representation of jet results in figure 5 of [START_REF] Coenen | Global instability of low-density jets[END_REF]. The absolute value of ψ is traced in logarithmic scale as a function of x. The phase of ψ is always chosen such that the real part of ψ is zero in the sensor location x s = 39, with the exception of the first mode (Fig. 4a), which has no interior wave nodes. Only positive real frequencies are considered.

∂ t ψ = Lψ → (iωI + L) -1 = ε -1 , (5)
Clearly, the eigenfunctions are characterised by integer numbers of wavelengths between the actuator and sensor locations, starting from zero wavelengths in the case of the lowest real frequency (Fig. 4a, only amplitude variations but approximately constant phase), to one wavelength in Fig. 4(b), and so forth with a continuously increasing count. The series could be continued further along the entire arc branch. The present model seems to reproduce very well the characteristics of arc branch eigenfunctions, as displayed in figure 5 of [START_REF] Coenen | Global instability of low-density jets[END_REF].

The position and spacing of eigenvalues along the branch curve, fixed approximately by an isocontour of the pseudospectrum, appears to be determined by a fitting phase relation between ψ at the actuator location and the applied feedback. To further illustrate this mechanism, the phase of the feedback coefficient C is varied. Figure 5 shows the arc branch as it is obtained with values C = 10 2 (bullets, same as before), C = 10 2 i (crosses), C = -10 2 (plus signs) and C = -10 2 i (squares). All symbols fall onto the same curve, confirming that the position of feedback eigenmodes, which represent singularities in the pseudospectrum, is fixed by the phase of the feedback relation, such that resonance can occur between both ends of the loop.

A parallel jet

The effect that imperfect numerical boundary conditions can have on eigenmode computations in a twodimensional flow domain is now investigated for the case of axisymmetric perturbations in a parallel round jet, governed by the incompressible axisymmetric Navier-Stokes equations, Boundary conditions are prescribed as

0 = ∂ r u r + u r r + ∂ x u x , (6a)
∂ t u r = -u r ∂ r u r -u x ∂ x u r -∂ r p + 1 Re ∂ rr u r + ∂ r u r r - u r r 2 + ∂ xx u r , (6b)
∂ t u x = -u r ∂ r u x -u x ∂ x u x -∂ x p + 1 Re ∂ rr u x + ∂ r u x r + ∂ xx u x . (6c)
u r = u x = 0 at x = 0, ( 7a
)
Re -1 ∂ x u x -p = 0 at x = x max , (7b)
u r = ∂ r u x = 0 at r = 0 and r = r max . (7c)
The stress-free outflow condition is a common and convenient choice for finite-element computations. Jet radius and centreline velocity are the dimensional length scales in this formulation, and the Reynolds number is chosen as Re = 100. The advantage of using a parallel base flow profile is that the global results can be rigorously compared with local instability properties. The standard analytical model of Michalke [START_REF] Michalke | Survey on jet instability theory[END_REF] is adopted,

U(r) = 1 2 1 + tanh 1 4θ 1 r -r , (8)
and a momentum thickness θ = 0.1 is chosen for this example.

The axisymmetric perturbation equations that follow from linearization of (6) about the base flow (8) are discretized with finite elements on a domain of length x max = 20 and radial extent r max = 50, using the FreeFEM++ software. This domain is resolved with 80 equidistant elements in x, and with 360 non-equidistant elements in r. The equations are solved for eigenmodes [u r , u x , p] T (r, x,t) = [ ûr , ûx , p] T (r, x) exp(-iωt), where the eigenvalue ω = ω r + iω i contains the angular frequency ω r and the temporal growth rate ω i .

The global spectrum is shown in Fig. 6. It features a clean upper arc branch (red) with maximum growth rate at ω 1 = 1.111 -0.052i. This mode, labelled '1', is chosen for further analysis. A lower branch (blue) is also present, from which the mode labelled '2' with ω 2 = 0.912 -1.211i will be examined. Similar lower branches are visible in wake spectra [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF] and in boundary-layer calculations [START_REF] Åkervik | Global two-dimensional stability measures of the flat plate boundary-layer flow[END_REF].

Projection onto local instability modes

For comparison, the corresponding local instability problem, for a domain of infinite extent x ∈ (-∞, ∞), is solved with a standard Chebyshev collocation technique on a staggered grid [START_REF] Khorrami | A Chebyshev spectral collocation method using a staggered grid for the stability of cylindrical flows[END_REF], using a coordinate transformation adapted for jet profiles [START_REF] Lesshafft | Linear impulse response in hot round jets[END_REF]. Identical radial collocation point distributions are used in the global and local computations, such as to eliminate the need for interpolation. Consistent boundary conditions (7c) are imposed in the local problem.

Spatial local instability modes are computed for the ω values corresponding to the global modes labelled in Fig. 6. The flow is convectively unstable, with an absolute frequency ω 0 = 1.074 -0.286i. Both direct and adjoint modes are solved for. The adjoint local modes represent the dual basis associated with the set of direct modes, and they serve for projecting the spatial structure of the global mode onto the local direct modes.

This projection is carried out in the following way, similar to the procedure used by Rodríguez et al. [START_REF] Rodríguez | A study of linear wavepacket models for subsonic turbulent jets using local eigenmode decomposition of PIV data[END_REF]: at a given streamwise station x, the radial variations of the global mode perturbation quantities are extracted. Since the eigenvector of the spatial local problem contains auxiliary variables k ûr and k ûx [see 18], these must also be added to the extracted slice of the global mode. This is accomplished by computing the streamwise derivative of the global ûr and ûx fields, and by augmenting the extracted vector [ ûr , ûx , p] T with [-i∂ x ûr , -i∂ x ûx ] T . This augmented global slice is finally projected onto the local modes, via scalar multiplication with the associated adjoint modes. As usual, the adjoint modes are normalised beforehand in such a way that their scalar product with the associated direct mode is unity, whereas the direct modes are by themselves scaled to have unit norm.

0 k + 1 k + 2 k - 1 k r k i (a) 0 1 2 3 4 -4 -2 0 k + 2 k + 3 k + 4 k ± k r k i (b)
The spatial local spectrum of the parallel jet profile, shown in Fig. 7 for the frequency values of global modes 1 and 2, is composed of discrete k + and k -modes [START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF], which represent downstream-and upstreampropagating hydrodynamic perturbations inside and near the jet. The labels in Fig. 7 rank all k + and k -modes according to their spatial growth rate -k i .

Figure 8 displays absolute values of the projection coefficients, obtained for global mode 1, pertaining to the three dominant local modes. These are the first two k + modes and the first k -mode. Blue symbols represent the k + 1 mode (the only one displaying unstable spatial growth). Except very near the Dirichlet inlet, the streamwise variation of this local mode amplitude is perfectly exponential, with a spatial growth rate 0.4604, as measured by a regression fit over 1 x 19. This value matches within 0.01% the imaginary part of the local eigenvalue. The amplitude of the global mode component ûx on the jet centerline is shown as a black line for reference.

Green symbols in Fig. 8 denote the amplitude of k + 2 , and red symbols the amplitude of k - 1 . Straight lines indicate the corresponding growth rate of the local eigenvalue for comparison. It is seen that both projections follow the amplitude variations expected from local analysis in a region close to one boundary, where their amplitude is maximal. Clearly, the k + 2 mode originates at the upstream boundary, whereas the k - 1 mode is forced at the downstream end. Farther away from those boundaries, both projections approximately follow the slope of the dominant k + 1 branch. This behaviour results from imperfections in the numerical projection, which apparently only allows a clean distinction between these local modes down to amplitude ratios around 10 -3 in the present setup. Higher spatial resolution in the local and global computations does not improve this threshold. Note that these three modes are highly non-orthogonal, which makes their distinction numerically delicate.

Overall, local mode contributions other than from the k + 1 mode to the global mode 1 are tractable but rather negligible. Modulations of the centerline velocity perturbation (black line in Fig. 8) are instead attributable to global pressure modes, as will be shown later on.

Very similar results are obtained for the global mode 2, for which the local mode amplitudes are displayed in Fig. 9. However, the local spectrum in this case differs from that of mode 1, as the global frequency ω 2 = 0.912 -1.211i has an imaginary part below that of the absolute frequency ω 0 = 1.074 -0.286i. Therefore, in the analysis of global mode 2, the local spatial modes are selected from a spectrum where pinching of the k + 1 and the k - 1 has already occurred (see Fig. 7b). In this setting, it is now the k + 2 mode that displays the strongest downstream growth, followed by the k + 3 mode. A mode from the mixed branch, formed from the k + 1 and k - 1 branches after pinching, is denoted k ± . The global mode 2, as represented by a black line in Fig. 9, is clearly dominated by the k + 2 mode (blue bullet symbols), but the k + 3 and k ± modes are again discernible down to amplitudes three orders of magnitude below k + 2 . It is not clear a priori how the k ± mode is to be interpreted, in particular with regard to its up-or downstream propagation. However, the projection results plainly show that this mode is generated at the downstream end, from where it propagates upstream; thereby, it behaves as a mode of k -type.

It is stressed that the analysis of spatial branches below the absolute growth rate is indeed meaningful in the present context. Spatial analysis below the absolute growth rate, i.e. after pinching has taken place, is usually said to be in violation of temporal causality, formally expressed by the fact that no integration path can be found in the complex k-plane that separates k + and k -branches [START_REF] Huerre | Open shear flow instabilities[END_REF]. This argument however arises in the context of the asymptotic flow behaviour at long times, when indeed the system dynamics are determined by the absolute mode. In the same sense, a global system is asymptotically determined by only the most unstable eigenmode (here: global mode 1). Notwithstanding, global eigenmodes with lesser growth rate do exist, and they are observable in the transient system dynamics. The spatial local modes used in the analysis of global mode 2 are justified in the same way.

In summary, the projections onto local spatial modes demonstrate that global mode 1 is supported by a local k + 1 mode, whereas global mode 2 relies on a k + 2 mode. In all likelihood, the same holds true for any global mode of the arc branch (k + 1 ) and of the lower branch (k + 2 ). Moreover, although not shown in Fig. 6, further lower branches exist at lower growth rates ω i < -1.5, for which a match with higher local k + branches is anticipated. The link between the arc branch and the dominant k + mode has been pointed out in earlier studies [for instance [START_REF] Ehrenstein | On two-dimensional temporal modes in spatially evolving open flows: the flat-plate boundary layer[END_REF][START_REF] Åkervik | Global two-dimensional stability measures of the flat plate boundary-layer flow[END_REF][START_REF] Garnaud | Modal and transient dynamics of jet flows[END_REF]. However, the argument so far has one loose end: the presence of a local k + mode is contingent on it being forced upstream. The present results show that this forcing takes place immediately at the upstream boundary. The essential ingredient that can give rise to a global mode with a k + wave is feedback from downstream.

Global pressure feedback

The ellipticity of the global linear jet problem is contained in the pressure gradient and in the viscous terms. The latter will only be noticeable over distances much shorter than the numerical box length, and they are not considered in the following analysis. The pressure, as noted by Ehrenstein & Gallaire [START_REF] Ehrenstein | On two-dimensional temporal modes in spatially evolving open flows: the flat-plate boundary layer[END_REF], obeys a Poisson equation, which in the present case of parallel flow takes the form

∆ p = -2∂ r U∂ x ûr , (9)
with homogenous Dirichlet and Neumann conditions at the upstream and lateral boundaries, respectively, and with p = Re -1 ∂ x ûx at the outflow.

The following analysis is restricted to the arc branch mode labelled '1' in Fig. 6, but results for mode 2 are not fundamentally different. The pressure amplitude of mode 1 is shown in Fig. 10(a) as log 10 | p|. Its structure is somewhat irregular in the downstream near-field region of the jet, yet the characteristic wavelength 2π/k r = 4.1 of the k + 1 mode is apparent. The stress-free outflow boundary condition is seen to result in p ≈ 0, and the pressure at r 5 is essentially a superposition of fundamental solutions of the homogeneous equation ∆ p = 0, with Dirichlet conditions at the inflow and outflow. These are given by

p j = A j I 0 jπ L r + B j K 0 jπ L r sin jπ L x , j ∈ N, (10)
where I 0 and K 0 are the modified Bessel functions of the first and second kind, respectively, and L = 20 is the streamwise length of the numerical box. The K 0 functions are exponentially decaying in r, manifestly dominant in the present problem, whereas the I 0 functions grow exponentially in r. These only enter the global mode at very low amplitude (A j /B j 1) in order to satisfy the Neumann condition at r max = 50. Only the p1 component is clearly visible in Fig. 10a, because it experiences the slowest radial decay, but a projection confirms that at least the first five p j components enter the pressure field with comparable global norms.

It is known from the analysis in Sect. 3.1 that the global mode involves a strong k + 1 wave. In the present section, the complementary part is sought that may provide the upstream-reaching part of a feedback loop. Therefore, the k + 1 component of the global mode is subtracted from the pressure field, using the already known projection coefficients shown in Fig. 8. This 'stripped' pressure field p is presented in Fig. 10(b), as log 10 | p|. A remarkably clean structure is recovered, suggestive of a solution of the Laplace equation ∆ p = 0 which is forced at the outflow boundary near the jet axis. A small inhomogeneity is also observed at the inflow near the the jet axis.

Note that the local k + 1 mode by construction represents a particular solution of (9) in a domain of infinite streamwise extent. Within the limits of the simplifying assumption that the propagating perturbations in the region of ∂ r U = 0 are indeed given by the k + 1 mode alone, in the interior of the bounded domain, the stripped pressure field p only needs to satisfy the inhomogeneities that arise from the boundary conditions. These inhomogeneities on both ends of the domain are thus coupled through the Laplace equation in p.

As a final plausibility check, the amplitudes of the supposed downstream-propagating and upstreamreaching components of the feedback loop are compared in Fig. 11 along a path at constant r. The radial position of the critical point of the local k + 1 mode is chosen, r c = 0.92. A thin blue line represents the pressure amplitude of the k + 1 mode component, according to the projection carried out in Sect. 3.1, a thick red line represents the amplitude of the pressure feedback signal p, and black symbols mark the amplitude of the total pressure field p. The up-and downstream branches have approximately equal amplitude at the downstream boundary, nearly cancelling each other. Near the upstream boundary, the feedback signal is larger by about a factor 3 than the k + 1 wave -this appears reasonable in view of the assumption that the latter is forced by the former.

The influence of box length

Arc branch modes in the literature are consistently found to be sensitive to the length of the numerical box. In some instances [e.g. 12], a longer box yields eigenvalues with lower growth rates, in other instances [e.g. 10] the opposite effect is observed. If the present feedback model is correct, it should allow an estimation of the influence of the numerical box length on the arc branch growth rates.

Eigenmodes of the parallel jet (8) have been computed in numerical domains of streamwise lengths L = 10, 15 and 25, for comparison with the standard configuration L = 20; in all cases, the radial discretisation and the constant step size ∆ x are unchanged. Resulting modes of the arc branch are displayed in Fig. 12(a). As L is increased, the entire branch is seen to shift to higher growth rates, and the modes are more densely spaced.

The analysis so far suggests that a downstream-convecting k + mode and the elliptic pressure field are coupled in small regions near the inflow and outflow boundaries. In analogy with the Ginzburg-Landau model discussed in Sect. 2, an effective reflection coefficient C may be defined, which relates the forcing of the k + mode at the inflow to the k + amplitude at the outflow, thus lumping the narrow interaction regions into singular actuator and sensor positions in x. The effective forcing at the Dirichlet inflow boundary is modelled as

-iω q = C q e ikL → C = -iωe -ikL , ( 11
)
where ω is the global eigenmode frequency, and (k, q) denote the relevant local k + mode. The coefficient C includes the effects of coupling on both ends of the domain, as well as the upstream decay of the pressure signal, seen in Fig. 11.

Numerically obtained values of C are reported in figures 12(b,c), in terms of their absolute value and their phase, for the arc branch in all four numerical domains. In each domain, the reflection coefficient decreases with real frequency, but only weakly so for frequencies larger than one. This effect may be caused by the slower radial decay of q at low frequencies, resulting in radially more extended interaction regions on both domain ends. The phase of C varies only slightly with frequency, and it is independent of box length.

If the reflection coefficients are scaled with the cube of the box length, as shown in Fig. 12(d), they neatly collapse onto one curve. Consequently, the effective feedback imparted by boundary reflections decreases with the domain length as L 3 . This scaling factor indicates that the relevant component of the pressure signal, which couples the downstream with the upstream end of the k + wave, is of a quadrupole type (signals from monopole and dipole sources at x = L would decay as L and L 2 , respectively).

It is then straightforward to interpret the effect of an increased numerical box length on a given arc branch mode: if the exponential growth of the dominant k + mode over the added streamwise interval is larger than the algebraic decay of the upstream-reaching pressure signal, over the same added distance, then the spurious inlet forcing will increase in strength, resulting in a higher growth rate ω i . If the k + mode in the added region is locally stable, then the inlet forcing will decrease, and the global growth rate will be lower as a result.

This interpretation appears to be consistent with the two cited examples: the boundary layer investigated by Ehrenstein & Gallaire [START_REF] Ehrenstein | On two-dimensional temporal modes in spatially evolving open flows: the flat-plate boundary layer[END_REF] is locally unstable at the outflow, as confirmed by the authors, and longer domain sizes are found to result in higher global growth rates. In contrast, the rapidly spreading jet considered by Garnaud et al. [START_REF] Garnaud | Modal and transient dynamics of jet flows[END_REF] is locally stable at the outflow with respect to axisymmetric perturbations, and indeed lower global growth rates are obtained for the arc branch in longer domains. This criterion may of course be frequency-dependent: certain (complex) frequencies may be locally stable at the outflow while others are unstable, and opposite trends ought to be observed in these frequency ranges. This seems indeed to be the case in the Re = 360 setting of Coenen et al. [START_REF] Coenen | Global instability of low-density jets[END_REF].

Possible remedies, tested for a non-parallel jet

While the preceding analyses of the Ginzburg-Landau equation and of the parallel jet served to characterise the global feedback mechanism due to domain truncation, the question how such feedback may be reduced is addressed in this section for the example of a spatially developing jet. With this choice, the results from the previous section can be largely transferred to the new setting.

Each type of open flow may present particular problems in view of domain truncation. In favourable configurations, the numerical boundaries can be placed in stable flow regions, as for instance in a uniform flow upstream of a solid obstacle, like a cylinder. In other cases, the need for truncation may be avoided by prescribing consistent physical boundaries, like confining solid walls, which are then part of the flow configuration.

Jets belong to a more problematic category. They are created by some upstream source of momentum, and it is in general not feasible to include the entire upstream apparatus (fan, chamber, nozzle, etc.) in the calculations, which would furthermore defeat the purpose of any generic description of jet dynamics. The formulation of upstream flow and boundary conditions is therefore necessarily imperfect with respect to any flow realisation, as potentially important regions are not accounted for.

For the present study, the jet is modelled as issuing from an orifice in a solid wall, very similar to the configuration of Garnaud et al. [START_REF] Garnaud | Modal and transient dynamics of jet flows[END_REF]. Upstream of the orifice, the flow develops in a straight circular pipe, from where it exits with a fairly thin boundary layer. The base flow is computed as a steady, axisymmetric solution of the incompressible Navier-Stokes equations [START_REF] Chomaz | Bifurcations to local and global modes in spatially developing flows[END_REF], with an inflow condition

U(r) = tanh 5 2 1 r -r (12)
imposed at x = -10. Newton-Raphson iterations are performed in order to converge to a steady flow state. The numerical domain for these base flow calculations is truncated downstream at x max = 50 and radially at r max = 50, where stress-free conditions are applied. The Reynolds number is set to 1000 in the free jet. Table 1 Numerical configurations of the various eigenmode calculations discussed in this section.

However, in order to maintain a thin shear layer at the orifice, this value is varied exponentially inside the pipe, between Re = 10 5 at x = -10 and Re = 10 3 at x -0.2, in the base flow calculation. The resulting flow field near the orifice at x = 0 is represented in Fig. 13, together with the characteristic streamwise variations of the shear layer momentum thickness θ m and of the centreline velocity.

Eigenmode calculations are carried out on smaller domains, where portions of the base flow are cropped at the inflow, the outflow, or both, in order to probe the effect of domain truncation. All configurations are listed in table 1. The Reynolds number in these calculations is set to 1000 throughout the domain.

The influence of upstream truncation on eigenmodes is considered first. Figure 14a shows spectra for three different domains, with homogeneous Dirichlet conditions (7a) imposed at x min = (0, -5, -10), respectively. The downstream end of the numerical domain is placed at x max = 40, where stress-free conditions (7b) are applied in all three cases. Unconverged and lower-branch eigenvalues are not shown here and in the following for clarity; the criterion for convergence is that an eigenvalue could be reproduced within three-digit accuracy using two different shift values.

It is seen from Fig. 14a that the arc branch is quite insensitive to the position of the upstream Dirichlet boundary. Even a full truncation of the upstream pipe (case 1) only results in a slight stabilisation at high frequencies. Eigenvalues obtained with pipe lengths of 5 and 10 radii (cases 2 and 3) are virtually identical. Stress-free inflow conditions (cases 4, 5 and 6, spectra shown in Fig. 14b) are found to perform less favourably. When applied at x = 0, these conditions give rise to global instability; the inclusion of portions of the pipe has a stabilising effect, but even with x min = -10 the growth rates along the entire arc branch are still significantly higher than those obtained with Dirichlet conditions. It is to be expected that global pressure fluctuations induce a generation of vorticity waves at the solid corner at x = 0, and the independence of the eigenvalues in cases 2 and 3 (Fig. 14a) on the upstream boundary position suggests that this physical effect is dominant over spurious pressure-vorticity coupling at x min . It remains to be determined to what extent the downstream boundary conditions emit spurious pressure signals, and how these may be reduced.

Garnaud et al. [START_REF] Garnaud | Modal and transient dynamics of jet flows[END_REF] tested stress-free and convective outflow conditions in an almost identical jet configuration, and found no significant difference in the resulting spectra. In the present study, the potential of absorbing layers [START_REF] Colonius | Modeling artificial boundary conditions for compressible flow[END_REF] for a stabilisation of the arc branch is examined. To this end, an artificial damping term is added to the linear perturbation equations, with the purpose of reducing perturbation amplitudes before they reach the numerical boundary at x max . If the original global eigenvalue problem is written

-iωB q = L q, ( 13
)
then the eigenvalue problem with absorbing layer is defined as

-iωB q = [L -λ (x)B] q. ( 14
)
For cases 7 and 8 (see table 1), the damping parameter is prescribed as Fig. 16 Amplitude variations in x of selected eigenfunctions, as marked in Fig. 15. The amplitude is defined in Eq. ( 16). This variation of λ is continuous in its first derivative, and it has maximum values 0.5 and 0.75 at x max = 40 and x max = 50, respectively.

λ (x) =    0 for x 19, 0.00625x 2 -0.2375x + 2.
The effect of absorbing layers on the eigenvalue spectrum is documented in Fig. 15. Damping in the downstream region 19 x 40 (case 7) is seen to reduce the growth rates of arc branch modes significantly (compared to case 2). Eigenfunctions of selected modes, marked by circles in Fig. 15, are compared in Fig. 16. An integral amplitude measure is defined as

A(x) = 1 0 | ûr | 2 + | ûx | 2 r dr 1 2 , ( 16
)
and the amplitude curves in Fig. 16 are normalised with respect to their values at the outflow boundary for comparison. While the red curve (case 2, no absorbing layer) shows a monotonic growth of perturbation amplitude throughout the domain, similar to the parallel jet results in figures 8 and 9, the blue curve (case 7) reaches a maximum inside the absorbing layer and subsequently decays in x. Remarkably, the ratio A(x max )/A(0) is identical in both cases, which is consistent with the interpretation of these modes as being the result of spurious feedback from the outflow, in the same way as detailed in Sect. 3 for the parallel jet. The absorbing layer reduces the gain of the k + hydrodynamic branch of the feedback loop, and thereby the temporal modal growth rate.

Based on the discussion in Sect. 3, stronger artificial damping of the hydrodynamic branch ought to lead to ever smaller modal growth rates; longer domains should furthermore lead to a reduced amplitude of the incident pressure signal upstream, due to its cubic decay. Both conditions are combined in case 8 (see table 1), of which results are included in figures 15 and 16. Indeed, the temporal growth rates are further reduced, and the spatial decay of the eigenfunction amplitude in the absorbing layer is more pronounced than in the case 7.

However, the shape of the arc branch in case 8 displays some differences with respect to all other cases, and similar mode patterns have been reported from jet calculations on long domains by Coenen et al. [START_REF] Coenen | Global instability of low-density jets[END_REF], where no artificial damping was applied. In the vicinity of the least stable mode (circled in Fig. 15), a regular spacing in ω r is still observed, but with larger distances between consecutive modes than in cases 2 and 7. The corresponding amplitude function in Fig. 16 shows that the ratio A(x max )/A(0) is smaller than in the two other cases, whereas the cubic decay of the reflection coefficient with box length, as described in Sect. 3, should have resulted in a larger ratio. Furthermore, stronger damping through higher values of λ , tried in test calculations that are not shown here, does not decrease the maximum growth rate much further, but it quickly leads to ill-conditioned system matrices. Case 8 seems indeed to mark an efficiency limit of absorbing layers for the stabilisation of the arc branch.

The above observations, in particular the increased spacing of modes in ω r , suggest that pressure feedback in the strongly damped case 8 may originate from the interior of the domain, possibly from the location of the amplitude maximum. Such feedback may be spurious, in the sense of Heaton et al. [START_REF] Heaton | Global linear stability of the non-parallel Batchelor vortex[END_REF], or it may be physical. This hypothesis is suggested for further examination.

Conclusions

A family of linear instability eigenmodes, named the arc branch, has been analysed in view of its physical or numerical origin. Branches of this type have previously been observed in a large variety of open flows. All results presented herein lead to the conclusion that arc branch modes in incompressible flow calculations are an artifact of domain truncation, due to spurious pressure feedback between numerical inflow and outflow boundaries.

A Ginzburg-Landau model was investigated first, because in this setting the effect of explicitly prescribed feedback between a downstream sensor and an upstream actuator could be examined without ambiguity. In the presence of feedback, a branch of eigenmodes was observed to arise that exhibits all typical characteristics of the arc branch described by Coenen et al. [START_REF] Coenen | Global instability of low-density jets[END_REF] for the example of light jets. These modes align with a contour of the pseudospectrum, spaced at regular intervals of the real frequency, and their spatial eigenfunction is characterised by an integer number of wavelengths between actuator and sensor locations. As the strength of the feedback is increased, the branch moves steadily upward to higher growth rates in the complex frequency plane. Eigenmodes of the feedback-free system that lie below the arc branch are no longer detectable. An important observation is that this arc branch in the Ginzburg-Landau equation with feedback has no counterpart in the feedback-free spectrum: these eigenmodes are not merely affected by the presence of feedback, indeed without it they do not exist.

An incompressible parallel jet in a truncated domain was examined next. No explicit feedback was prescribed in this setting, but a similar arc branch of eigenmodes was nonetheless found to dominate the spectrum. A lower branch of eigenmodes with stronger temporal decay was also described, which may correspond to subdominant branches observed in wakes and boundary layers. Both the arc-branch and the lower-branch modes in the parallel jet were shown to be composed primarily of one downstream-propagating k + wave, as computed from a local spatial analysis. Each of the two branches involves a different k + mode. This observation raises the question how a k + wave is generated at the upstream boundary of the numerical domain. The analogy with the Ginzburg-Landau model from Sect. 2 suggests feedback from downstream.

The most plausible mechanism for global feedback in a truncated domain is the generation of pressure perturbations due to the artificial downstream boundary condition, as described by Buell & Huerre [START_REF] Buell | Inflow/outflow boundary conditions and global dynamics of spatial mixing layers[END_REF]. The pressure field of the least stable arc branch mode was decomposed into one portion associated with the prominent k + wave, which cannot be involved in upstream feedback, and into a residual part that is essentially governed by a Laplace equation forced at the domain boundaries. The latter appears to be strongly dominated by a source region at the outflow, near the jet axis. A numerical evaluation of effective reflection coefficients, in analogy with the Ginzburg-Landau model of Sect. 2, established as a principal result that the strength of upstream feedback decays with the numerical box length to the third power. This scaling is indicative of a spurious pressure quadrupole situated at the outflow. The algebraic nature of the pressure decay allowed a prediction of the effect of box length variations on arc-branch growth rates, depending on the local stability or instability of the flow, apparently consistent with common observations. While this study has been limited to incompressible flow settings, similar mechanisms will also be present in compressible calculations, with the difference that spurious pressure signals are then propagated by a wave equation. This may in fact result in stronger feedback, as the far-field acoustic pressure only decays with the first power of the distance from the source.

Finally, possible strategies for a reduction of spurious pressure feedback have been examined for the example of a spreading jet. On the one hand, upstream boundary conditions must be chosen that minimise the unphysical conversion of pressure feedback from downstream to vortical perturbations. It has been found, in this particular flow example, that Dirichlet conditions perform much better than stress-free conditions in this regard. On the other hand, the generation of spurious pressure signals from the outflow boundary must be reduced. It has been demonstrated that artificial damping in a downstream absorbing layer provides an efficient means to achieve this. An advantage of this technique is that it is straightforward to implement in any usual open shear flow problem. It has been noted, however, that increasingly strong damping does not reduce the growth rate of the arc branch to arbitrarily low levels.

Fig. 1

 1 Fig. 1 Spectrum and pseudospectrum of a jet at Re = 1000, showing a typical arc branch. Eigenmodes are drawn as black symbols. The definition of ε values, shown as colour contours, is given in Sect. 2. The figure is made from the data of Coenen et al. [7].

a) C = 10 - 10 b) C = 10 - 6 c) C = 10 Fig. 2

 1010106102 Fig. 2 Spectra of the Ginzburg-Landau equation (1) with various feedback coefficients C. Eigenvalues with closed-loop feedback (•), and analytical eigenvalues (4) of the case without feedback (•).

Fig. 3

 3 Fig. 3 Pseudospectra of the Ginzburg-Landau equation; (a) without feedback, C = 0, and (b) with feedback, C = 10 -6 . Contours of ε, defined in (5), are represented in logarithmic colour scale.

Fig. 4

 4 Fig. 4 Eigenfunctions of the feedback branch for C = 10 2 . Left column: position of the eigenvalue in the spectrum; right column: corresponding eigenfunction log 10 |Re(ψ)| along x. Zero wavelengths are present in the top frame, one in the second, and so forth up to four wavelengths in the bottom frame. The series continues further along the branch.

Fig. 5

 5 Fig. 5 Eigenvalues of the feedback branch for |C| = 10 2 with varying phase: (•) C = 10 2 , (×) C = 10 2 i, (+) C = -10 2 , ( ) C = -10 2 i.

Fig. 6

 6 Fig.6Spectrum of the parallel jet in a finite domain with inflow and outflow boundary conditions. Red: arc branch; blue: lower branch. Labelled modes are discussed in the following.

Fig. 7

 7 Fig. 7 Spatial eigenvalues from a local analysis of the parallel jet. Lines represent branches obtained by varying the real part of the frequency. a) global mode 1: ω = 1.111 -0.052i; b) global mode 2: ω = 0.912 -1.211i.

Fig. 8 Fig. 9

 89 Fig. 8 Projection of global mode 1 onto spatial local modes of frequency ω = 1.111-0.052i: absolute values of the projection coefficients as functions of x. Legend: (•) k + 1 ; (+) k + 2 ; (×) k - 1 ; ( ) global mode component ûx (0, x) on axis; colored lines indicate the growth rates according to local analysis.

Fig. 10 Fig. 11

 1011 Fig. 10 Pressure perturbation fields of global mode 1. a) total pressure perturbation log 10 | p|; b) pressure perturbation log 10 | p|, obtained by subtraction of the k + 1 component.

Fig. 12

 12 Fig. 12 Influence of the numerical box length on eigenvalues. In all figures, different symbols denote different box lengths: ( ) L = 25; (•) L = 20; ( ) L = 15; ( ) L = 10. a) Arc branch eigenvalues. b) Effective reflection coefficients according to (11), absolute values. c) Effective reflection coefficients, phase. d) Effective reflection coefficients, absolute value scaled with L 3 .

Fig. 13

 13 Fig. 13 Steady state of a spatially developing jet, used as a base flow in all calculations of Sect. 4. a) Contours of axial velocity near the orifice at x = 0. Contour values are distributed between 0 and 1, with constant spacing of 0.1. The base flow is computed on a larger domain, with -10 x 50 and 0 r 50. b) Streamwise development of the shear layer momentum thickness. c) Streamwise development of the centreline velocity.

  max = 0.75

Fig. 14 Fig. 15

 1415 Fig. 14 Arc branch spectra obtained in calculations with (a) Dirichlet and (b) stress-free upstream boundary conditions.

  25625 for 19 < x < 21, 0.025(x -20) for x 21.

Acknowledgements Discussions with Xavier Garnaud and Wilfried Coenen greatly helped shape the ideas presented in this paper. The study was supported by the Agence Nationale de la Recherche under the Cool Jazz project, grant number ANR-12-BS09-0024.


domain: spi phys

text:
Guillaume Walter

Jean-Baptiste Dherbecourt

Jean-Michel Melkonian

Myriam Raybaut

Cyril Drag

Antoine Godard

Spectral Modulations in a Picosecond OPO Based on a Chirped Quasi-Phase Matched Crystal

Spectral Modulations in a Picosecond OPO Based on a

Chirped Quasi-Phase Matched Crystal Aperiodic-or chirped-quasi-phase-matching (QPM) nonlinear crystals offer unique features compared to conventional quadratic nonlinear materials. Indeed, one can engineer the phase-matched wavelengths as functions of the position z along the propagation axis in the material, the QPM grating period Λ(z) being a function of z.

Used as the amplification medium in an optical parametric oscillator (OPO), a chirped QPM crystal can provide a broad parametric gain spectrum without the need to use short crystals with reduced conversion efficiency. Such OPOs based on chirped QPM were reported in various temporal regimes [START_REF] Tillman | Low-threshold, high-repetition-frequency femtosecond optical parametric oscillator based on chirped-pulse frequency conversion[END_REF][START_REF] Tillman | Monolithic optical parametric oscillator using chirped quasi-phase matching[END_REF][START_REF] Descloux | Spectrotemporal dynamics of a picosecond OPO based on chirped quasi-phase-matching[END_REF].

In picosecond synchronously pumped OPOs (SPOPOs) based on chirped QPM crystals and broadband mirrors (without intracavity filter), we observe signal spectra consisting of a main peak with lower side lobes [START_REF] Descloux | Wide and fast dispersion tuning of a picosecond OPO based on aperiodic quasi-phase matching using an axially chirped volume Bragg grating[END_REF]. These spectral modulations are not observed for similar picosecond SPOPOs based on uniform QPM crystals.

Here, we report on the experimental and numerical investigations of these spectral modulations in OPOs based on chirped QPM, which still need to be explained. We show that the spectral modulations are induced by cascaded three-wave mixing between different parts of the interacting-waves spectra, quasi-phase matched at different positions in the nonlinear crystal. As presented in Fig. 1(a), the SPOPO is based on an aperiodically poled lithium niobate (APPLN) crystal (HC Photonics) and is pumped by a Nd:YVO 4 laser with a repetition rate of 76 MHz and a pulse duration of 8 ps. The signal oscillating in the OPO cavity can be analyzed at the output of a plane mirror with a reflectivity of 98%. Here, experiments are carried out with an apodized chirped QPM grating, with a central chirp rate of 1.17 × 10 4 m -2 , leading to a parametric gain bandwidth of 40 nm and an oscillation threshold of 1 W.

By the use of an optical spectrum analyzer (OSA), the pump and signal spectra can be studied. Fig. 1(b) and Fig. 1(c) present the evolution of the spectra of the signal and the depleted pump at the output of the nonlinear crystal, as functions of the pump power. The higher is the pump power, the broader is the emitted signal spectrum with more and more side lobes. In the same time, modulations are also observed in the depleted pump spectrum. A previously validated home-made numerical code has been used to understand the development of such modulations [START_REF] Descloux | Spectrotemporal dynamics of a picosecond OPO based on chirped quasi-phase-matching[END_REF]. For that purpose, we model the evolution of the interacting parametric waves while they pass through the nonlinear crystal. Fig. 2 shows the evolution the signal and pump spectra during a propagation through the APPLN crystal when the OPO is in steady-state. From this figure, we can identify that wavelengths located either on peaks or on valleys of the output signal spectrum experience very different evolution while they propagate through the crystal. Indeed, the part of the spectrum which corresponds to modulation peaks experiences an exponential growth over a characteristic length, expected from theory [START_REF] Charbonneau-Lefort | Optical parametric amplifiers using chirped quasi-phase-matching gratings {I}: practical design formulas[END_REF]. Then, this part has an oscillatory evolution through the rest of the crystal. The part which corresponds to modulation valleys of the spectrum exhibits a different evolution. These wavelengths first experience a power decrease before being amplified. As observed in Fig. 2(b), this decrease corresponds to the generation of new spectral sidebands in the traveling pump spectrum where signal spectrum modulation valleys exhibit a power decrease.

Such a behavior is consistent with cascaded difference-and sum-frequency generation (DFG-SFG) conversion: idler wavelengths are generated by quasi-phase-matched DFG between incident signal and pump pulses, then quasiphase-matched SFG occurs at other positions in the crystal between the generated idler and other signal wavelengths which generate pump sidebands. Such cascaded processes are enabled by the variation of the QPM period over the length of the aperiodic nonlinear crystal.

The origin of the spectral modulations in OPOs based on chirped QPM crystals are thus explained for the first time owing to our study. Further details, additional experimental and numerical results will be presented.

Fig. 1 .

 1 Fig. 1. Schematic diagram of the SPOPO setup (a). Evolution of the signal spectrum (b), and of the pump spectrum (c), as a function of the pump power.

Fig. 2 .

 2 Fig. 2. Simulation at steady state of the parametric generation evolution of the signal spectrum (a), and of the evolution of the pump spectrum (b), through the APPLN crystal

Guillaume Walter 1 , Jean-Baptiste Dherbecourt 1 , Jean-Michel Melkonian 1 , Myriam Raybaut 1 , Cyril Drag 2 , Antoine Godard 1,*

  We investigate and model spectral modulations that are specific features of OPOs based on chirped quasi-phase matched crystals. Their occurrence is related to cascaded three-wave mixing processes that are quasi-phase matched at different positions..

	1 DPHY, ONERA, Université Paris Saclay, F-91123 Palaiseau, France
	2 Laboratoire de Physique des Plasma, Ecole Polytechnique-CNRS-Univ Paris-Sud-UPMC, Université
	Paris-Saclay, 91128 Palaiseau, France
	OCIS codes: (190.4970) Parametric oscillators and amplifiers, (190.4223)	Nonlinear wave mixing,
	(190.4410) Nonlinear optics, parametric processes

* antoine.godard@onera.fr Abstract:

Guillaume Walter, Jean-Baptiste Dherbecourt, Jean-Michel Melkonian, Myriam Raybaut, Cyril Drag, et al.. Spectral Modulations in a Picosecond OPO Based on a Chirped Quasi-Phase Matched Crystal. CLEO 2019, May 2019, SAN JOSE, France. OSA, 10.1364/CLEO_AT.2019.JTu2A.93. hal-

Funding

Agence Nationale de la Recherche (ANR) (ANR-14-CE26-0036); "Investissements d'Avenir" LabEx PALM (SYCLOP and SORA) (ANR-10-LABX-0039).


Denis S Grebenkov
email: denis.grebenkov@polytechnique.edu

Ralf Metzler

Gleb Oshanin

Leonardo Dagdug

Alexander M Berezhkovskii

Alexei T Skvortsov

Affiliations

https://doi.org/10.1063/1.5098390

The problem of trapping of diffusing particles by identical, perfectly absorbing rings periodically distributed along an otherwise reflecting infinite cylindrical tube is very difficult for analytical treatment due to the mixed boundary conditions on the tube surface. Analytical solutions to this problem are considered in two recent papers. 1,2 Among other results reported in these works, the authors propose approximate expressions for the mean lifetime τ of a particle whose starting point is uniformly distributed over the tube surface, and the tube is surrounded by an outer concentric tube which is a reflecting boundary for diffusing particles. These expressions can be written in a unified form as

τ = L R 2 -ρ 2 (2πDρ) f , (1)
where ρ and L are the radius and half-period of the inner tube, R is the outer tube radius, and D is the particle diffusivity. The factor f is, in general, a dimensionless function of three dimensionless geometric parameters: ρ/L, R/L, and ε = πl/L, where l is the half-width of the absorbing ring, and hence, l/L is the surface fraction of the inner tube covered by the rings. The two papers give different expressions for this factor, denoted by f 1 and f 2 . In this note, these expressions are tested against exact values of f obtained by solving the mixed boundary-value problem numerically by the finite element method.

Being motivated by search problems in biology, the authors of Ref. 1 studied the binding of a particle diffusing inside a finitelength outer cylinder (a "bacillus-shaped bacteria") to a partially absorbing ring located on the surface of a concentric inner cylinder in the geometry shown in Fig. 1. Such trapping problems arise, for example, in modeling protein binding to specific sequences on DNA molecules or on specific locations of the nucleolus that generates cascades of biochemical reactions that support the functioning of living systems. [START_REF] Alberts | Molecular Biology of the Cell[END_REF][START_REF] Berg | [END_REF][5][6][7][8] The focus of Ref. 1 is on the mean particle lifetime considered as a function of the particle initial position, geometric parameters of the system, the intrinsic reactivity of the ring surface, and the particle diffusivity. To solve the problem, the authors adopt the approximate approach proposed in Ref. 9. This allows them to obtain an analytical expression for the mean lifetime. In the special case where the ring is perfectly absorbing and the particle starting point is uniformly distributed over the surface of the inner cylinder, their general result reduces to the expression for τ in Eq. ( 1) with the factor f given by

f ≈ f 1 = 2 ∞ n=1 (Gn n)[sin nε (nε)] 2 , (2)
with

Gn = I 1 (πnR L)K 0 (πnρ L) + I 0 (πnρ L)K 1 (πnR L) I 1 (πnρ L)K 1 (πnR L) -I 1 (πnR L)K 1 (πnρ L) , (3)
where Iν(z) and Kν(z) are the modified Bessel functions of the first and second kind. By symmetry, the trapping problem analyzed in Ref. 1 is equivalent to that in the case of two infinite concentric cylindrical tubes with partially absorbing rings of width 2l periodically arranged on the inner tube, with period 2L. Trapping by infinite periodically striped cylindrical surfaces is studied in Ref. 2, assuming that the stripes are perfectly absorbing. The analysis is performed for three orientations of the stripe direction with respect to the tube axes: perpendicular, parallel, and at the angle π/4 to the axis. This is done for both internal and external problems, where the particles diffuse inside and outside the striped tube, respectively. The key approximation used in Ref. 2 is the so-called boundary homogenization which is the replacement of nonuniform boundary conditions on the surface by an effective radiation boundary condition with properly chosen reactivity that is uniform over the surface. Since a striped cylindrical surface can be formed by "rolling" a corresponding flat striped surface, it is hypothesized that its reactivity is the same as that of the flat surface, for which the effective surface reactivity is known. 10 For the mean particle lifetime discussed above, this leads to the expression for τ in Eq. ( 1) with the factor f given by f ≈ f 2 = 2 ln(1 sin(ε 2)), (4) which is much simpler than its counterpart in Eq. (2). In contrast to the factor f 1 , the factor f 2 depends only on ε = πl/L, i.e., on the strip surface fraction, and is independent of two other parameters, ρ/L and R/L. Note that f 1 and f 2 have the same asymptotic behavior in the limiting case of ε → 0, where the strip width tends to zero, and both expressions for the factor f diverge as 2ln(1/ε). As might be expected, both f 1 and f 2 vanish when ε → π since the entire inner tube surface becomes absorbing.

Below we compare the expressions in Eqs. ( 2) and ( 4) with the exact solution for the factor f obtained numerically with the aim to establish their accuracy and limitations. This is done for a wide range of inner and outer tube radii at two values of the ring surface fraction, when the rings occupy one tenth and one half of the inner tube surface, ε = π/10 and ε = π/2, respectively. The results of the comparison are presented in Fig. 2 that shows the contour plots of the absolute values of the relative error of the two expressions in percent, 100|1 -f i/f |, i = 1, 2, as functions of δ/L and ρ/L, where δ = Rρ is the distance between the inner and outer tubes. The relative error of Eq. ( 4) is shown in panels (a) and (b) for l/L = 0.1 and 0.5, respectively. Corresponding results for Eq. ( 2) are presented in panels (c) and (d). As ρ/L and δ/L tend to infinity, the relative errors approach their asymptotic values. For the expression in Eq. ( 4), this asymptotic value is zero (by construction) since the cylindrical geometry becomes effectively flat in this limiting case. The asymptotic value of the relative error of the expression in Eq. ( 2) is a function of ε: it is 8% at ε = π/10 and about 20% at ε = π/2. From Fig. 2, one can see that the simple expression in Eq. ( 4) fails for narrow-compared to the system period 2L-inner cylinders and when the outer tube is too close to the inner one-the distance δ is shorter or equal to the period 2L. Otherwise, i.e., when both ρ and δ exceed 2L, this expression accurately predicts the factor f. This is true for both values of the ring surface fraction. It is important that when Eq. ( 4) does not work, at small l/L, one can use Eq. ( 2) that predicts the factor f with high accuracy over the entire range of both parameters, ρ/L and δ/L, provided that the ring surface fraction is low.

To summarize, the expression for the factor f in Eq. ( 2) works better than its counterpart in Eq. ( 4) at small ring surface fractions. For example, at l/L = 0.1, the relative error of its predictions is less than 10% over the entire range of ρ/L and δ/L. As might be expected, this relative error increases with the ring surface fraction, but even at l/L = 0.5, it does not exceed 50%. As concerns the simple expression in Eq. ( 4), it is reliably applicable at an arbitrary ring surface fraction on condition that both ρ/L and δ/L are greater than unity. These two approximations complement each other.

This study was supported by the Intramural Research Program of the NIH, Center for Information Technology.

FIG. 1 .

 1 FIG. 1.Schematic view of the diffusion domain studied in Ref. 1. A particle diffuses between the two concentric cylinders of radii R and ρ ≤ R. The intercylinder space is constrained by two reflecting walls separated by distance L, which are perpendicular to the cylinder axis. The particle is trapped by a partially absorbing ring of length l located on the inner cylinder near the reflecting wall. The rest of the inner cylinder surface and that of the outer cylinder are reflective.

FIG. 2 .

 2 FIG. 2. Contour plots showing the absolute values of the relative error of the two approximate expressions for the factor f, given in Eqs. (2) and (4), in percent, 100|1 -f i /f |, i = 1, 2, as functions of the geometric parameters δ = R ρ and ρ, where δ is the distance between the inner and outer cylinders, rescaled by the system length L. The contour plots are presented for two values of the surface fraction occupied by the absorbing ring on the inner cylinder, l/L = 0.1 [panels (a) and (c)] and l/L = 0.5 [panels (b) and (d)]. The relative errors of Eqs. (4) and (2) are shown in panels (a) and (b) and panels (c) and (d), respectively. The numbers indicate the magnitudes of the relative error.

J. Chem. Phys. 150, 206101 (2019); doi: 10.1063/1.5098390 150, 206101-1

J. Chem. Phys. 150, 206101 (2019); doi: 10.1063/1.5098390 150, 206101-2


L Fruchter

V Brouet

F Brisset

H Moutaabbid

Y Klein

Growth facets of SrIrO 3 Thin Films and Single Crystals

The crystallographic orientation of SrIrO3 surfaces is decisive for the occurrence of topological surface states. We show from density functional theory computations that (001) and (110) free surfaces have comparable energies, and, correspondingly, we experimentally observe that single micro-crystals exhibit both facet orientations. These surfaces are found to relax over typically the length of one oxygen octahedron, defining a structural critical thickness for thin films. A reconstruction of the electronic density associated to tilts of the oxygen octahedra is observed. At the opposite, thin films have invariably been reported to grow along the (110) direction. We show that the interfacial energy associated to the oxygen octahedra distortion for epitaxy is likely at the origin of this specific feature, and propose leads to induce (001) SrIrO3 growth.

I. INTRODUCTION

Amongst the Ruddlesden-Popper series for iridates, R n+1 Ir n O 3n+1 where R= Sr, Ba and n = 1,2,∞, SrIrO 3 (n = ∞, SIO) occupies a peculiar place. While the n=1,2 compounds are antiferromagnetic insulators with a layered structure, SIO is a three-dimensional, nearly compensated paramagnetic semi-metal [1]. It was proposed that the semi-metallic state is stabilized as bands cross at a topologically-protected line of nodes below the Fermi level [2,3]. In SIO the Dirac line of nodes is enforced by symmetry in the nonsymmorphic space group Pbnm [4], and the key-role of n-glide operation was stressed as spin-orbit coupling is present [4,5]. As a consequence, the crystallographic nature of the surface of practical samples is decisive, when the existence of topological surface states is questioned [5][6][7].

Elaboration of bulk samples in the polycrystalline form, using high-pressure synthesis, have been reported up to now only [8,9], and epitaxy was used to stabilize the single crystal material. The perovskite structure for SIO makes it compatible with many other perovskite substrates commonly used in film growth. Indeed, SIO thin films have been grown using molecular beam epitaxy [10,11], pulsed laser deposition [5,[12][13][14][15][16][17][18], metalorganic chemical vapor deposition [19], and RF sputtering [20].

While rough structural studies of these films generally rely on the pseudo-cubic structure of the primitive cell of the material, there has been some detailed structural investigations of SIO and the parent ruthenate compound SrRuO 3 thin films. It has been reported that SrRuO 3 epitaxial layers are grown on (001) SrTiO 3 with out-of-plane [110] direction and [001] and [1][2][3][4][5][6][7][8][9][10] aligned in-plane [21][22][23][24]. A similar growth orientation was reported for SIO films on (001) SrTiO 3 , elaborated with PLD, as evidenced by electron diffraction observations [25], and for SIO grown on GdScO 3 , from a complete structure determination [5]. It was found that the substrate vicinality systematically orients the in-plane orientation of the films, with the [001] direction lying nearest to the vicinal steps one.

In order to gain some understanding on the mechanisms that drive the surface growth selectivity, we have investigated SIO free-standing single crystal surfaces, and computed surface energies and configurations, using ab initio density functional theory (DFT).

II. EXPERIMENTAL

A. Single crystals

We obtained single crystals from a SIO sintered pellet, grown under high pressure. The high pressure synthesis was realized using a Paris-Edinburgh VX3 hydraulic press, as described in details elsewhere [26]. In brief, powder of pristine hexagonal SIO was pressed into a 2.9-mm-diameter and 5.5-mm-high pellet and wrapped into a thin platinum foil. The latter was then introduced in a MgO capsule sleeve which was in turn inserted into a graphite cylinder acting as a heater. We applied a quasi-hydrostatic pressure of 5 GPa and a temperature of 1000 • C +/-50 • C. The capsule was kept at constant pressure and constant temperature for 30 min and then quenched to room temperature in a few seconds. We then slowly released the pressure at a rate of 0.5 GPa/h. To separate single crystals, the pellet was first finely grounded in dilute polymethyl methacrylate (PMMA) resist. A silicon wafer was then spin-coated, with a drop of this PMMA/SIO mixture. To eliminate the PMMA resist from the surface of the SIO grains, the wafer was finally etched using an oxygen plasma. As can be seen in Figs. 2 and3, several single crystals with well defined facets were glued on the wafer, with their upper flat surface parallel to it. Such a configuration, where the surface orientation is approximately known in the electronic diffraction geometry, allowed to determine the surface crystallographic orientation, using Kikuchi line maps patterns.

The patterns from the back scattered electron diffraction (EBSD) were obtained with the wafer normal oriented at 20 deg. from the electron beam, in a Zeiss Supra 55 VP scanning electron microscope fitted with a Hikari EBSD System from TSL-EDAX TM microscope. We assumed that the single crystals had two parallel surfaces aligned with the wafer (Figs. 2 and3), to infer the orientation of the incident beam with respect to the crystal surface. The accelerating voltage was 20 kV and the current a few nA. The Kikuchi maps were analyzed using OIM TM software.

B. DFT computations

In order to compute the surface energies, we have performed DFT computation of the electronic energy of periodic slabs, using the full-potential linearized augmented plane wave (FP-LAPW) method, as implemented in the Wien2k software [27].

As pointed out by Heifets et al [28], the surface energy is the sum of the cleavage and relaxation energies. One may then virtually cleave a bulk cell made of N structure unit cells, and compute the difference in energy between the cleaved slabs and the bulk material, for both the unrelaxed and relaxed structures. As for SrTiO 3 , we expect that, due to the covalent nature of the Ir-O bonding, as compared to ionic character of the Sr-O one, the SrO-IrO 2 stacking is weakly polar with a negative charge in the SrO plane, although formal valences assign the compound to Tasker's classification type I [29]. In order to impede the formation of an electrostatic macroscopic dipole, which could induce a spurious reconstruction of the surface, it is preferable to create two symmetric slabs (i.e. with identical terminations at both ends), each of them terminated with different atomic planes [30]. These terminating layers may be SrO or IrO 2 (Fig. 1).

To investigate the (001) surface, we built couples of slabs consisting of the stacking along the c-axis of SrO-IrO 2 -...-IrO 2 -SrO and IrO 2 -SrO-....-SrO-IrO 2 (001) layers. Depending on the number of unit cells contained in the original bulk cell (N ), the slabs could be made of N/2 unit cells (N even), or (N -1)/2 unit cells (N odd). While the bulk structure has the Pbnm symmetry [9,31] (No. 62 -8 symmetry operations), the cleaved slabs have a lower orthorhombic symmetry, obtained as the maximal non isomorphic subgroup of this group with the preserved symmetry elements: -P2 1/c (No. 14 -4 symmetry operations, including inversion center) for N even, IrO 2 -terminated cells, and N odd, SrO-terminated slabs, -Pmc2 1 (No. 26 -4 symmetry operations, including a c-axis mirror plane) for N even, SrO-terminated cells, and N odd, IrO 2 -terminated slabs. In the case of (110) surface, the [110] direction makes an angle γ = 89.6 deg. from the (110) plane. As a consequence, the slabs have a lower, monoclinic symmetry, with one mirror symmetry, Pm (No. 6). The unit cell is also larger by ≈ √ 2 in the basal plane, with unit cell parameter 7.98 Å along the [-110] direction, and 7.948 Å along the [001] one. The vacuum between each cell was chosen 10 Å wide. We used for the bulk structure the one obtained by relaxing both the cell parameters and the internal coordinates. As can be seen from Table I, the fully relaxed structure has oxygen octahedra slightly more tilted than found in the literature, and larger cell parameters. To be consistent with our free-surface cell relaxation procedure, we used this relaxed structure to build the cleaved slabs. For the consistency, also, between energy and structural relaxation computations, muffin-tin radii were systematically reduced by 8% for all computations. Structural relaxations were performed down to residual force 0.2 mRy/a.u. We found that a k-mesh with a number of k-points as small as 500 was sufficient to insure the energy convergence of the cells, with accuracy well within the measured energy differences. Given the total electronic energy for cells as defined above, the surface energy may be computed using:

E(SrO/IrO2) = E cleavage + E rel (SrO/IrO2) E cleavage = 1/4[E unrel (SrO) + E unrel (IrO2) -N E bulk ] E rel (SrO/IrO2) = 1/2[E rel (SrO/IrO2) -E unrel (SrO/IrO2)]
where E rel/unrel (SrO/IrO2) is the total electronic energy for the relaxed/unrelaxed cell, with SrO/IrO 2 termination, and E bulk is the electronic energy for one (relaxed) bulk unit cell. Very well contrasted Kikuchi maps were obtained from the selected crystals, indicating clean, metallic surfaces, with a well defined crystallographic orientation (Fig. 4). The Kikuchi diffraction lines result from the diffraction by the atomic planes of the diffused electron beam (a mechanism similar to diffused X-ray diffraction for the Kossel lines), that is typically from a few tens nanometers thick layer from the surface. The pole figure for a collection of grains (Fig. 5) only shows two dominant orientations for the top surfaces, with 75% of them in the (110) orientation, and the rest in the (001) one. Considering the degeneracy of the (±1±10) surfaces and (00±1) ones for the Kikuchi maps, the expectation ratio is 2 for equiprobable occurrence of such orientations. Due to the limited number of available samples (32), the excess of our (110) observations cannot invalidate the hypothesis of an equiprobable occurrence for each orientation with more than ≈ 60 % confidence. So, we consider our observations to be compatible with equiprobable occurrence within a reasonable error bar. To compute the surface energies, we first looked for possible finite size effects. Fig. 6 displays the surface energy for (001) slabs of different thicknesses. Energies depend very weakly on the slab thickness, down to slabs with only 3 IrO 2 /SrO layers: this is an indication that the electronic and structural properties relax over quite a short distance from the surface, as will be seen below. The relaxation energy contributes to typically 30 % of the total energy, as is commonly the case. The SrO surface energy is systematically larger than the IrO 2 one, by about 3-7 %. Similar, but less extensive computations -due to the lack of symmetry, which considerably increases the computation time and the convergence accuracy -for (110) slabs were performed. As may be seen from table II, surface energies are comparable within error bars to the ones for (001) slabs. This is in agreement with our observation of comparable occurrence for (001) and ( 110) facets for single crystals. The top layers of a relaxed (001) IrO 2 -terminated slab are displayed in Fig. 7. It may be seen that the oxygen octahedra of the vacuum-interface layer rotate around the c-axis, so as to increase the Ir-O bas -Ir angle strongly (+8.3 deg.), while they tilt from this axis, so as to decrease the Ir-O ap -Ir angle (-3.5 deg.). The structural relaxation vanishes in the bulk, with a typical length which may be defined from the approximate exponential decay of the averaged atomic displacements in each layer. As can be seen in Fig. 9, this length is as short as about 1.3 oxygen octahedra: structural relaxation effects should manifest themselves in the electronic properties of films thinner than this limit, as well as for electronic probes below this depth (such as photo-emission). Indeed, the inspection of the electronic iso-electronic density surfaces for the top IrO 2 layer reveals that the tilt of the oxygen octahedra tends to disrupt the electronic cloud at some Ir-O-Ir bonds, where this bound has the largest bending (154.5 deg. for the weak bounds, as compared to 162.1 deg. for the strong ones in Fig. 10b). This results in the formation of chains along the b-axis, which illustrates the lower local symmetry of the surface, whereas uniform Ir-O-Ir angles in the bulk insures equivalent bonding in the IrO 2 plane (Fig. 10a)). This adds a free surface structural limit, in addition to electronic confinement effects and epitaxial constraint ones in thin films [33,34]. Relaxed (001) SrO-terminated slabs display smaller and opposite displacements (Fig. 7): the octahedra rotation about the c-axis increases (-0.7 deg. for the Ir-O bas -Ir angle), while the tilt from this axis decreases (+2.4 deg. for the Ir-O ap -Ir angle).

Relaxed (110) surfaces show similar trends as for (001) surfaces: the IrO 2 -terminated slabs display an increase of the Ir-O bas -Ir angle and a decrease Ir-O ap -Ir angle, while SrO-terminated slabs display opposite rotations. An essential difference in this case is due to the presence of the (001) mirror plane, which imposes two basal oxygen to act as pivots with respect to c-axis rotations (Fig. 8). As for the (001) slabs, the rotation of the oxygen octahedra disrupts the electronic cloud, inducing in this case chains along the (001) direction (Fig. 11). However, due to the smaller rotations, the electronic constrictions are less marked in this case.

An alternative to quantify the octahedra rotations is to generalize Glazer's rotation angles [35] to the slab surface (assuming that the surface layer is the one of a bulk crystal with all symmetry elements, and neglecting octahedra distortions). The rotation angles are shown in Table III. These observations may be qualitatively understood considering the extra charge brought by the surface termination. As the charge on the Iridium site is increased (IrO 2 -terminated surface) or decreased (SrO-terminated surface), the Ir-O bond length is expected to increase (alt. to decrease). A direct consequence for this is a smaller (alt. larger) tilt of the octahedra along the axis perpendicular to the IrO 2 plane.

It must be stressed that the above computations are made within the constraint of the given slab cell periodicity and maximal symmetry. While the system does relax towards a local energy minimum in the phase space (which cannot be entirely explored with Wien2k package, which relies on a steepest descent algorithm), we implicitly discarded superstructures and lower symmetries. More stable configurations, with a different periodicity and/or a lower symmetry, may exist. However, the relaxation of a (001) slab with a doubled unit cell (2a,2b), and no crystal symmetry yielded the same configuration as in Fig. 7. So, the relaxed surface appears robust against periodicity and symmetry conditions. This is an indication, also, that symmetry elements needed for topological surface states can be preserved. Finally, we would like to discuss the consequences of the above computations for the growth of SIO thin films. One generally expects that the free surface energy, and the elastic one associated to the epitaxy, can be decisive for the growth orientation of thin films. As we have seen, the former does not discriminate in favor of one or the other orientation. Considering the elastic energy associated to the material strain, we notice that the (001) epitaxy distinguishes from the (110) one, as it requires only normal strains along the orthorhombic crystal axis, while (110) epitaxy requires both a normal strain along the [001] direction, and a shear strain in the (110) plane. We have evaluated the stiffness coefficients for SIO, by computing the energy of strained structures [36]. We find: C 11 ≈ C 22 ≈ C 33 ≈ 260 GPa, C 13 ≈ C 23 ≈ C 12 ≈ 120 GPa, and for the shear coefficient: C 66 ≈ 90 GPa. This evidences the nearly cubic character of the structure. In the framework of the generalized Hooke's law for the linear theory for small deformations, and for a cubic material, it may be shown that the elastic energy, W , for the (001) epitaxy is, assuming that the film relaxes freely along the [001] direction:

W/ε 2 = C + C ′ -2 C ′2 /C ( 1
)
where ε is the strain in the substrate plane, assumed identical in the

W/ε 2 = C (1 + 2 A 2 )/2 -C ′ A + 2 C 66 (1 + A) (2)
where A = C ′ /(C + C ′ ). Using the evaluations given above, the stiffness for the (001) growth, W/ε 2 , is found 270 GPa from Eq. 1, to be compared with the one for the (110) growth, 360 GPa, from Eq. 2. Thus, we find that it is not energetically favorable for stressed films to involve a shear strain, as is the case for the (110) growth, rather than only normal strains, as for the (001) growth. However, the difference in elastic energy between the two orientations is quite small, as compared to typical surface energies. Using the above estimates, a 1 % mismatch for a one-unit cell thick film yields ∆W ≈ 0.5 meV/ Å2 , so we do not expect this contribution to be pertinent, for usual mismatch values.

So, neither a difference in the free surface energy, nor one in the elastic energy associated to the epitaxial constrain, can account for the observed orientation of epitaxial thin films. We have not yet considered, however, the contribution of the epitaxial interface to the film energy, which may be decisive for the orientation of the first layer, and of the subsequent ones. Beyond matching of the film and substrate pseudo-cubic cells, which is usually crudely considered for epitaxy (and obviously not selective for the film orientation), the tilting of the oxygen octahedra required at the interface is potentialy a large, orientation-dependent, contribution to the interface energy. The evaluation of this contribution requires to investigate bilayers, which is beyond the scope of this paper. However, we may roughly compare substrates, using the displacement of the apical oxygen atom, needed to match a given substrate, once the material is strained to accommodate the pseudo-cubic cells. As may be seen in Table IV, amongst the best matched substrates for SIO growth, SrTiO 3 shows the lowest selectivity with this criterion, while there is a large difference in favor of the (110) growth on scandates substrates with the usual (110) orientation. We included LaMnO 3 , although it is not a common substrate. Indeed, (001) LaMnO 3 is close to the pseudocubic form, with Mn-Mn angle 90.02 deg. (while Ir-Ir = 90.3 deg. for (001) SIO), allowing for an alternative to the selective growth of (001) SIO, as seen in Table IV. (001) DyScO 3 and (001) GdScO 3 , which show a favorable selectivity, are quite far from the pseudocubic form (respectively, Sc-Sc = 92.8 and 92.7 deg.), but a large descrepancy between the pseudocubic cells may be acceptable, as is the case for TbMnO 3 (Mn-Mn = 95.7 deg.) grown on DyScO 3 [START_REF] Daumont | Multiferroic perovskites under epitaxial strain: the case of TbMnO3 thin films[END_REF]. In summary, we have shown that (001) and (110) SrIrO 3 surfaces have comparable surface energies, which is confirmed by the observation that single micro-crystals exhibit both facet orientations with comparable probability. The atomic displacements at the surface were found to relax in the bulk over typically one oxygen octahedron, and to determine electronic channel patterns that run on the cleaved surfaces. For films, we found that usual strains cannot induce a preferential surface orientation, and we suggest that the epitaxy is solely responsible for the observed (110) orientation.

FIG. 11: Iso-electronic density surface (d = 0.12 (a.u.) -3 ) in the IrO2 plane for a IrO2-terminated (110) slab (the red lines highlight the paths with the largest electronic density). The right panel displays the surface with the same density as in Fig. 10, d = 0.1 (a.u.) -3 .

FIG. 1 :

 1 FIG. 1: Symmetric cells with (001) IrO2/SrO free surface, obtained from the cleavage of N = 4.5 bulk unit cells. The black lines delimit the unit cells.

  FIG. 2: A selection of crystals with (110) top surface. The cross indicates the diffraction point.

FIG. 4 :

 4 FIG.4: Indexed Kikuchi diagrams for bottom right crystals in Fig.2(left) and 3 (right). The labels are the crystal directions at the intersection of lattice planes.

FIG. 5 :

 5 FIG. 5: Inverse pole figure showing the orientation of the facet parallel to the substrate for 32 crystals (indices for Pbnm notation).

FIG. 6 :

 6 FIG.6: Energies for (001) surfaces. The abscissa is the number of IrO2 (alt. SrO) layers in the IrO2-terminated (alt. SrO) slab.

FIG. 8 :

 8 FIG. 8: Relaxation of the top layers for a (110) IrO2-terminated (left) and SrO-terminated slab (right). (top view)

  [100] and [010] directions, C = C 11 = C 22 = C 33 , and C ′ = C 12 = C 13 = C 23 . For the (110) epitaxy, assuming also identical strain in the [001] and [1 10] directions, we have:

3 FIG. 9 :

 39 FIG.9: Average atom displacement in an IrO2 layer, for a N =8 IrO2-terminated slab (the surface layer is indexed 0). The line is an approximation for an exponential decay of the N = 4 data.

TABLE I

 I

				: SrIrO3 structural parameters
		Ref. [9] Ref. [31] Ref. [32] This work
	a ( Å)	5.5617 "	5.6008	5.622
	b ( Å)	5.5909 "	5.5712	5.663
	c ( Å) Ir-O bas. -Ir (deg.) a 153.5 7.8821 " 154.0 Ir-Oap.-Ir (deg.) a 156.6 154.5	7.8960 156.2 156.9	7.948 153.8 153.3
	E/unit cell (eV) b 0.22	0.73	0.31	0

a O bas./ap. = basal-plane / apical oxygen b referenced to this work

TABLE II :

 II Relaxation and total surface energies for (001) and (110) slabs (meV/ Å2 ) E

		rel (001)	E(001)	E rel (110)	E(110)
	IrO2	-17.9 ± 0.1	64.1 ± 0.7	-14 ± 0.4	66 ± 4
	SrO	-14.5 ± 0.1	67.5 ± 0.7	-18 ± 1	63 ± 4

FIG. 7: Relaxation of the top layers for a (001) IrO2-terminated (left) and SrO-terminated slab (right) (top view). The length of the arrows is proportional to the displacement amplitude.

TABLE III :

 III Glazer's rotation angles for the relaxed slabs surface layer (deg.)

		θ [100]	θ [010]	θ [001]
	bulk	9.8	9.8	8.6
	(001)-IrO2	11.5	11.5	16.8
	(001)-SrO	9.4	9.4	16.8
	(110)-IrO2	13.4	6.7	17.2
	(110)-SrO	9.7	8.0	10.7

TABLE IV :

 IV Root mean square displacement of apical oxygen for SrIrO3 epitaxially grown on different substrates ( Å).

			Favorable
	situations for (001) growth are highlighted
	substrate	(001) SrIrO3	(110) SrIrO3
	(100) SrTiO3	0.41	0.46
	(110) GdScO3	0.57	0.22
	(110) DyScO3	0.94	0.22
	(001) GdScO3	0.37	0.65
	(001) DyScO3	0.39	0.65
	(110) LaMnO3	0.38	0.08
	(001) LaMnO3	0.04	0.75
		IV. CONCLUSIONS

Acknowledgements

We thank Dorothée Colson and Anne Forget from SPEC-CEA-CNRS for providing us with the monoclinic SrIrO 3 material. We acknowledge support from the Agence Nationale de la Recherche grant "'SOCRATE"' (ANR-15-CE30-0009-01), and LabEx PALM ReflectX (ANR-10-LABX-0039-PALM).