Antonio Celani

Gautam Reddy

Massimo Vergassola

Sticky issues in turbulent transport

We investigate the role of partial stickiness among particles or with a surface for turbulent transport. For the former case, we re-derive known results for the case of the compressible Kraichnan model by using a method based on bi-orthogonality for the expansion of the propagator in terms of left and right eigenvectors. In particular, we show that enforcing the constraints of orthogonality and normalization yields results that were previously obtained by a rigorous, yet possibly less intuitive method. For the latter case, we introduce a general model of transport within the atmospheric boundary layer. As suggested by experimental observations on the transport of atmospheric tracers, both drift and diffusivity scale with the height to the ground. The strength of the drift is parameterized by a velocity V . We use the bi-orthogonality method to show that for V in the range -1 < V < 0 and 0 < V < 1 there is a one-parameter family of boundary conditions that are a priori admissible. Outside of that range, there is a single boundary condition that is admissible. In physical terms, the one-parameter family is parametrized by the degree to which particles stick to the ground.

This paper is dedicated to the memory of Krzysztof Gawȩdzki, whom we had the honor and the chance to know and collaborate with. The goal of our paper is to revisit the effect of boundary conditions on particle dispersion in turbulent flows. Two situations will be considered: i) the relative dispersion of particle pairs in a compressible flow, and ii) the single particle dispersion in an incompressible flow with a rigid boundary. These two problems have in common a singular dynamics at the boundary which leads to interesting behaviours. For pair dispersion, the spatial roughness of the velocity field and the degree of compressibility conspire to generate nontrivial time-dependencies of the rate at which particles separate. In some regimes, the statistics of pair dispersion on macroscopic scales depends on the microscopic details of the interaction between particles when they are in contact with each other. In the seemingly very different case of single-particle dispersion in an incompressible, wall-bounded flow, the fact that the velocity vanishes upon approaching the boundary leads to a similar effect: depending on the degree of stickiness between particle and boundaries, different dispersion behaviours arise. As discussed in detail below, these two problems share a common mathematical structure. Before delving into the technical part, let us offer a synopsis of the main results:

• For pair dispersion in a compressible flow, we re-derive and expand on results that were obtained in Gawȩdzki's work [START_REF] Gawȩdzki | Phase transition in passive scalar advection[END_REF]. There, the trajectories of fluid particles in a compressible generalization of the Kraichnan ensemble of turbulent velocities were considered. The Kraichnan ensemble features Gaussian stochastic flows with power-law spatial correlations and short temporal correlations. The ensemble has played a key role in the first systematic understanding of anomalous scaling in turbulent flow (see [START_REF] Shraiman | Scalar turbulence[END_REF][START_REF] Falkovich | Particles and fields in fluid turbulence[END_REF] for an introduction and review). Ref. [START_REF] Gawȩdzki | Phase transition in passive scalar advection[END_REF] showed that, depending on the degree of compressibility, trajectories could either explosively separate or implosively collapse. These two cases result in different statistical properties for scalar quantities passively advected by the flow. In the first case, a direct cascade of the energy of a scalar tracer with a short-distance intermittency and dissipative anomaly is observed. In the second case, implosive collapse leads to an inverse cascade of the tracer energy with suppressed intermittency. The relevant result here is the number of admissible boundary conditions at vanishing distance between the particles. A single boundary condition is admissible for two ranges of parameters (degree of compressibility, scaling exponent of the field and space dimension), whilst a one-parameter family of boundary conditions is a priori possible for an intermediate range of parameters. This family of boundary conditions takes into account the dynamics of particles that tend to aggregate upon contact and is parameterized by the degree of stickiness. The above result was obtained in Ref. [START_REF] Gawȩdzki | Phase transition in passive scalar advection[END_REF] by noting that the principal symbol of the operator controlling the separation between particles vanishes at the origin so that the operator loses ellipticity. The operator was then mapped onto a Schrödinger case known in the literature and the corresponding mathematical literature expounded in the book [START_REF] Reed | Methods of Modern Mathematical Analysis, II. Fourier Analysis, Self-Adjointness[END_REF] was used. While rigorous, the explanation was arguably not intuitive for a physics audience, including the co-author in Ref. [START_REF] Gawȩdzki | Phase transition in passive scalar advection[END_REF] and the present paper. In Section I we re-derive via bi-orthogonality the results that were used in [START_REF] Gawȩdzki | Phase transition in passive scalar advection[END_REF].

• For single-particle dispersion in a compressible flow confined by a rigid wall, we consider a model for the transport of odors [START_REF] Celani | Odor landscapes in turbulent environments[END_REF] in the atmospheric boundary layer (ABL) [START_REF] Kaimal | Atmospheric Boundary Layer Flows: Their Structure and Measurement[END_REF][START_REF] Garratt | The Atmospheric Boundary Layer[END_REF][START_REF] Wyngaard | Turbulence in the Atmosphere[END_REF]. The olfactory system of many animals and insects is exquisitely sensitive to odors that are transported in the environment by atmospheric turbulence [START_REF] Wyatt | Pheromones and Animal Behavior[END_REF][START_REF] Riffell | Olfactory ecology and the processing of complex mixtures[END_REF]. For instance, male moths respond to minute amounts of pheromones emitted by females, and are sensitive to the finescale structure of turbulent plumes where pheromone concentration is detectible [START_REF] Murlis | Odor Plumes and How Insects Use Them[END_REF][START_REF] Reddy | Olfactory sensing and navigation in turbulent environments[END_REF]. While pheromone whiffs were known to be intermittent, a quantitative characterization of their statistical properties was lacking, which motivated the study of this challenging fluid dynamics problem. This problem is also relevant for entomology, neurobiology and the design of olfactory stimulators that reproduce physiologically relevant odor signals in well-controlled laboratory conditions [START_REF] Celani | Odor landscapes in turbulent environments[END_REF]. In this work, a Lagrangian approach was used to model the transport of odors by turbulent flows and predict the statistics of odor detections during olfactory searches. General scaling arguments were developed to compute the probability distributions of odor intensity, the durations of odor detections and the durations between odor detections. In addition, an explicit model for the dynamics of odor transport in the ABL was developed. In the original model of the ABL in Ref. [START_REF] Celani | Odor landscapes in turbulent environments[END_REF], drift and diffusivity both scale linearly with the height from the ground. This choice was motivated by experimental observations on the transport of atmospheric tracers. However, an additional assumption relating the magnitudes of drift and diffusivity was made, which led to simplified, time-reversible Lagrangian dynamics in the height z.

Here, we remove this simplifying assumption and introduce a general model for the dynamics of tracers (namely odors) transported in the ABL. Drift and diffusivity are kept linear in z but their coefficients are not constrained. The lack of this constraint leads to an additional drift term in the Fokker-Planck equation for the model, which we parameterize by a velocity V . Our first goal is to show that this term plays a role analogous to the role played by compressibility in Ref. [START_REF] Gawȩdzki | Phase transition in passive scalar advection[END_REF]. Specifically, for V in the range -1 < V < 0 and 0 < V < 1 there is a one-parameter family of boundary conditions that are a priori admissible. For V = 0 and for V outside of that range, there is a single boundary condition. Our second goal is to reformulate the derivation of such results by showing that they emerge from the request of a bi-orthogonal expansion for the propagator of the particle dynamics.

After having introduced the model for transport in the ABL in Section II, we derive the equation for the marginal probability and one possible solution in Section II A. The issue of multiple solutions is tackled in Section III, where we demonstrate the constraints imposed by bi-orthogonality for the ABL model. The selection of the solution in the range where multiple ones are a priori possible is discussed in Section IV.

I. RE-DERIVING CLASSICAL RESULTS BY THE REQUEST OF BI-ORTHOGONALITY

The procedure followed in Ref. [START_REF] Gawȩdzki | Phase transition in passive scalar advection[END_REF] is to map the problem of the transport of particles by a Kraichnan compressible flow onto the following equation for the (to be determined) function R(z) :

∂ 2 z - b 2 -1 4 z 2 R (z) = 0 , (1)
with |b| > 0 (see Section 35 in [START_REF] Landau | Quantum Mechanics[END_REF] for its physical discussion). The equation has been studied in the mathematical literature and admits rigorous results on its possible boundary conditions [START_REF] Reed | Methods of Modern Mathematical Analysis, II. Fourier Analysis, Self-Adjointness[END_REF]. Specifically, the mathematical statement is that Eq. ( 1) admits a single possible boundary condition for |b| ≥ 1 whilst for |b| < 1 it admits a one-parameter family of possible solutions. Our goal here is to re-derive this result by showing that it amounts to requesting biorthogonality of the expansion for the propagator of the dynamics. The formalism will then be used in the following Sections to analyze transport of particles in the ABL. A bi-orthogonal expansion for the propagator P(z, z 0 ; t) of a linear operator L reads [START_REF] Risken | The Fokker-Planck Equation. Methods of Solution and Applications[END_REF] :

P(z, z 0 ; t) = dEe -Et φ E (z 0 ) ψ E (z) , (2)
where φ and ψ are the left and right eigenfunction, respectively, that correspond to the eigenvalue -E of L. Using the expression of P, an initial condition f (z 0 , 0) is propagated at time t into

f (z, t) = dEe -Et f (z 0 ) |φ E (z 0 ) ψ E (z) , (3)
where | denotes the scalar product relevant to the problem under consideration. The equations for the right/left eigenvectors ψ/φ are

Lψ E = -Eψ E ; L † φ E = -Eφ E , (4)
where L † is the adjoint operator defined with respect to the scalar product mentioned above. The condition of bi-orthogonality for the expansion Eq. ( 2) amounts to the request that the scalar product

φ E1 |ψ E2 = δ (E 1 -E 2 ) . (5)
Let us now consider the case of the classical Eq. ( 1). Since the operator L defined by Eq. ( 1) is self-adjoint, its left and right eigenvectors are identical. The corresponding equation ( 4) is :

Lψ E ≡ ∂ 2 z - b 2 -1 4 z 2 ψ E (z) = -Eψ E (z) , (6)
which can be mapped onto a Bessel form by the transformation ψ E = z 1/2 χ √ Ez . The corresponding equation for χ is indeed [START_REF] Gradshteyn | Tables of Integrals, Series and Products[END_REF] :

χ + χ z + 1 - b 2 z 2 χ = 0 , (7)
which is the Bessel equation for J ±b (z). It follows that left and right eigenvectors are

ψ E = φ E = z 1/2 J ±b √ Ez .

A. Conditions for bi-orthogonality

Let us now consider conditions for bi-orthogonality of the expansion ( 5) and show that they yield the aforementioned results on possible boundary conditions.

Convergence of the integrals.

We first consider the range of b values that ensure convergence of the integrals in Eq. ( 5), with the scalar product:

φ E1 |ψ E2 = φ E1 (z) ψ E2 (z) dz.
Since J b (z) z b at the origin [START_REF] Gradshteyn | Tables of Integrals, Series and Products[END_REF], it follows that the pair z 1/2 J b (z) for φ and ψ is acceptable only for b > -1. Similarly, the pair z 1/2 J -b for both φ and ψ is acceptable only for b < 1. Mixed pairs z 1/2 J b and z 1/2 J -b behave as z at the origin and thus have no convergence constraints.

We now consider integer b's, where J b and J -b are not independent. For b = n positive integer, one solution is J n and the other independent solution is Y n . The first behaves as z n close to the origin, while the second behaves as z -n [START_REF] Gradshteyn | Tables of Integrals, Series and Products[END_REF]. It follows that the pair Y n for both φ and ψ is ruled out. The other three pairs are a priori acceptable. The case b = -n, yields the same constraints due to the identities J -n = (-1) n J n and Y -n = (-1) n Y n .

Orthogonality.

Let us now inquire which pairs of Bessel functions ensure orthogonality. By considering φ E1 |Lψ E2 -Lφ E1 |ψ E2 , we obtain that

(E 2 -E 1 ) φ E1 (z) ψ E2 (z) dz = φ E1 ψ E2 -φ E1 ψ E2 ∞ 0 . (8)
Non-trivial conditions arise from the behavior at the origin, as discussed hereafter.

Let us first consider the "homologous" pair φ

E1 (z) = z 1/2 J b √ E 1 z and ψ E2 (z) = z 1/2 J b √ E 2 z
. Each boundary term on the r.h.s. of Eq. ( 8) scales as zz b z b-1 ∝ z 2b , which is not vanishing for b ≤ 0. However, the two non-zero terms appearing on the r.h.s. are identical and therefore they cancel out. Subdominant terms then become relevant. Since the next subdominant term in the expansion of J b is ∼ z b+2 , we obtain that the r.h.s. in Eq. ( 8) behaves as ∝ z 2b+2 , which vanishes in the range of convergence of the integrals b > -1. The same argument holds for the pairs

φ E1 (z) = z 1/2 J -b √ E 1 z and ψ E2 (z) = z 1/2 J -b √ E 2 z
, which are orthogonal for b < 1. We now consider integer values of b, for which Y b is the second independent solution. For positive integers b = n, the homologous pairs of J n are the same as above and the homologous pair of Y n does not ensure convergence of the integrals, as it was shown above. Finally, negative integers b = -n are like positive integers due to the identities

J -n = (-1) n J n and Y -n = (-1) n Y n .
Let us now consider "heterologous" pairs for non-integer b's, that is pairs φ E1

(z) = z 1/2 J -b √ E 1 z and ψ E2 (z) = z 1/2 J b √ E 2 z
(or with the sign of b inverted). Explicit calculation of the boundary term of Eq. [START_REF] Wyngaard | Turbulence in the Atmosphere[END_REF] shows that the r.h.s. does not vanish. Finally, we consider again integer values of b. For positive integers b = n, the pair of eigenvectors is φ

E1 (z) = z 1/2 Y n √ E 1 z and ψ E2 (z) = z 1/2 J n √ E 2 z (or J and Y inverted)
. By using Y n (z) z -n and J n (z) z n , we check that the boundary term in Eq. ( 8) does not vanish. Finally, negative integers b = -n are like positive integers due to the identities J -n = (-1) n J n and Y -n = (-1) n Y n .

In summary, by imposing conditions for bi-orthogonality, we have recovered the known result [START_REF] Reed | Methods of Modern Mathematical Analysis, II. Fourier Analysis, Self-Adjointness[END_REF] that a single solution is admissible for |b| ≥ 1, while two possible choices are possible for |b| < 1. We shall next proceed to exploit this procedure to analyze the dynamics of particles transported in the atmospheric boundary layer.

II. A GENERAL MODEL FOR THE TRANSPORT OF ODORS IN THE ABL

A relevant physical environment for olfactory search is the near-neutral atmospheric surface layer [START_REF] Kaimal | Atmospheric Boundary Layer Flows: Their Structure and Measurement[END_REF][START_REF] Garratt | The Atmospheric Boundary Layer[END_REF][START_REF] Wyngaard | Turbulence in the Atmosphere[END_REF]. The turbulent wind flow in the neutral boundary is characterized by the following features :

(i) The mean wind depends logarithmically on the height above the ground. In practice, as the logarithmic factor varies slowly, we will consider a uniform mean wind U hereafter ;

(ii) The intensity of velocity fluctuations u is nearly constant with height and is small compared with the mean wind ;

(iii) The size of the turbulent eddies at height z is proportional to z and their correlation time is z/u .

These properties motivate two approximations (a) The action of turbulent eddies on a transported particle can be subsumed in an effective diffusion coefficient and an effective drift that depend linearly on z. This implies that Eulerian time-correlations are neglected.

(b) Fluctuations in the position along the mean wind are neglected.

According to these approximations the transport equation for the probability distribution p(y, z, t) for the position (U t, y, z) of a particle advected by the flow takes the non-dimensional form

∂P ∂t + V ∂P ∂z = ∂ ∂y z ∂ ∂y + ∂ ∂z z ∂ ∂z P (9)
This equation must be equipped with an initial condition and boundary conditions at the ground z = 0. The choice of boundary conditions is the most delicate aspect since the diffusivity is vanishingly small as z → 0 and the flux to the boundary is affected by the drift V . In Ref. [START_REF] Celani | Odor landscapes in turbulent environments[END_REF], the problem with V = 0 was considered. The simplification brought by this choice is best grasped by noting that the action corresponding to Eq.( 9) is :

S(T ) = T 0 ẏ2 (t) + ( ż(t) -V ) 2 4z dt , (10)
with time and the downwind coordinate related as t = (x -x 0 )/U . The choice V = 0 then corresponds to a reversible dynamics, a symmetry which is broken by the more general (and realistic) choice of leaving V free.

In the sequel, we will study in detail the major effect that V has on the possible boundary conditions. To avoid clutter, we focus on the solution for the marginal distribution Q(z, t) = dyP (y, z, t), that obeys :

∂ t Q + V ∂ z Q = ∂ z (z∂ z Q) . (11)
The reduction to a one-dimensional problem simplifies the analytical calculations without sacrificing the nontrivial physics at the boundary.

A. One possible solution to the equation for the marginal probability Equation ( 11) admits an analytical solution, as we proceed to show by reducing it to a heat equation. To that purpose, we define

Q (z, t) = z α q √ z, t , (12)
and the variable ξ ≡ √ z. The equation ( 11) for α = V reduces then to

∂ t q (ξ, t) = 1 4 1 ξ 2V +1 ∂ ξ ξ 2V +1 (∂ ξ q) , (13)
which is indeed the heat equation in a space having dimensionality d = 2V + 2. A radially-symmetric initial condition solution q (ξ , 0) (see below) is propagated via the (d-dimensional) heat kernel into a radially-symmetric solution (which does not necessarily satisfies the boundary conditions at ξ = 0) at time t as

q (ξ, t) = q (ξ , 0) e -(ξ-ξ ) 2 t (πt) d/2 dξ , (14)
where ξ = |ξ| and ξ = |ξ |. The proportionality factor in the expression of the spherically-symmetric initial condition q (ξ , 0) ∝ δ (ξ -ξ 0 ) is determined by the normalization condition of the original Q(z, t). Requesting that Q (z, 0) dz = 2 ξ 2V +1 q (ξ, 0) dξ = 1 fixes the proportionality factor to 1/2ξ 2V +1 0 . Inserting this expression into Eq. ( 14) and expressing the integral in spherical coordinates gives

q (ξ, t) = e -ξ 2 +ξ 2 0 t 2 1 πt d/2 π 0 e 2ξξ 0 cos(θ) t (sin θ) d-2 dθ × Ω d-1 = e -ξ 2 +ξ 2 0 t √ πt d 1 Γ d-1 2 π 0 e 2ξξ 0 cos(θ) t (sin θ) d-2 dθ , (15)
where

Ω d = 2π d/2 Γ(d/2
) is the surface of the unit sphere in d dimensions. The integral in Eq. ( 15) is proportional to the integral representation of the modified Bessel function I (Eq. 8.431.1 in [START_REF] Gradshteyn | Tables of Integrals, Series and Products[END_REF]) :

I ν (z) = (z/2) ν √ π Γ (ν + 1/2) π 0 e z cos(θ) (sin θ) 2ν dθ . (16)
By identifying ν = d/2 -1 and using d = 2V + 2, i.e., V = ν, we conclude that

Q (z, t) = e -ξ(z) 2 +ξ 0 (z) 2 t t ξ(z) ξ 0 (z) V I V 2ξ(z)ξ 0 (z) t , (17)
where ξ(z) = √ z and ξ 0 (z) = √ z 0 .

While the solution (17) was a priori derived for integer d, i.e., for integer values of V , its final expression holds for any value of V . This is a consequence of the definition of I V or it can be directly verified as follows. The ansatz

Q (z, t) = e -z+z 0 t t z z 0 V /2 Q (z, t) , (18)
is immediately inspired by the solution Eq. ( 17). Inserting it into Eq. [START_REF] Murlis | Odor Plumes and How Insects Use Them[END_REF] gives

Qt + V 2 4z + z 0 t 2 Q = 1 - 2z t Qz + z Qzz . ( 19
)
The change of variable Q(z, t) = q 2 √ z0z t ≡ q (ρ) then reduces the above equation to

ρ 2 q ρρ + ρ q ρ -V 2 + ρ 2 q = 0 , (20)
which is the equation of the modified Bessel functions I V valid for any V (and the independent solutions K V or I -V , see below).

The distribution ( 17) is normalized for any V > -1 :

Q (z, t) dz = e -z 0 2t z 0 t -V +1 2 M -V +1 2 , V 2 z 0 t = e -z 0 t Φ V + 1, V + 1 ; z 0 t = 1 . (21)
Here, the first equality is formula 6.643.2 in [START_REF] Gradshteyn | Tables of Integrals, Series and Products[END_REF], the second is the definition 9.220.2 of the Whittaker function M in terms of the confluent hypergeometric Φ, and the last equality is formula 9.215.1 in [START_REF] Gradshteyn | Tables of Integrals, Series and Products[END_REF] : Φ (α, α ; z) = e z . Note that the inequality V > -1 reflects the request of integrability stemming from the behavior close to the origin :

Q(z, t) 1 Γ (V + 1) z t V e -z/t t , (22)
obtained by keeping only the first term of the expansion for I as z → 0. This expression becomes exact and valid over the entire range of z if z 0 = 0 as all other terms of the expansion of I are then vanishing.

III. OTHER SOLUTIONS AND THEIR SELECTION BY THE REQUEST OF BI-ORTHOGONALITY

The formulation that leads to Eq. ( 20) for modified Bessel functions, as well as the existence of the constraint V > -1 for Eq. ( 17), suggest that the original problem a priori admits solutions other than Eq. ( 17). It is therefore important to characterize them, and identify the corresponding boundary conditions that select them. That is the scope of this and the following Section, where we will show that the solution in Eq. ( 17) is a zero-flux solution and that in a range of parameters a non-zero flux solution will pop out of the bi-orthogonality, thereby providing a way to fix sticky boundary conditions.

A. The equations for the eigenfunctions in the ABL model

In the case of our model Eq. ( 11), the Fokker-Planck operator L is

L = [-V ∂ z + ∂ z (z∂ z )] , (23)
which is not self-adjoint. It follows that its left and right eigenfunctions do not coincide. The equation for the right eigenfunctions ψ

Lψ E = -Eψ E , (24)
is reduced to the Bessel form

ρ 2 ψ + ρ ψ + Eρ 2 -V 2 ψ = 0 , (25)
by defining ψ E (z) = z V /2 ψ (2 √ z), making the change of variable ρ ≡ 2 √ z and indicating by primes derivatives with respect to ρ. It follows that ψ is a combination of the two independent solutions

J V 2 √ Ez and Y V 2 √ Ez (or J -V 2 √ Ez for V non-integer).
Similarly, the equation for the left eigenfunctions involves the adjoint L † , which flips the sign of

V , V → -V . It follows that φ E (z) = z -V /2 φ 2 √
Ez , where φ is a combination of J V and Y V (or J -V for V non-integer).

B. Conditions for bi-orthogonality

We consider the conditions for bi-orthogonality of the expansion (5) for the ABL operator defined by Eq. (23).

Convergence of the integrals.

A first point in the analysis of Eq. ( 5) is the range of V values that ensure convergence of the integrals. We start by non-integer V 's and remark that J V ( √ z) z V /2 at the origin. It follows that the pair J V for φ and J V for ψ is acceptable only for V > -1. Similarly, the pair J -V J -V is acceptable only for V < 1. Mixed pairs J V J -V are constants around the origin and have no convergence constraint.

Let us now consider the case of integer V 's. For V = 0, one solution J 0 behaves as a constant close to the origin and the other independent solution Y 0 behaves as a logarithm, so that all three possible pairs do not feature any convergence constraints. For V = n positive integer, one solution is J n 2 √ Ez and the other independent solution is Y n 2 √ Ez . The first behaves as z n/2 close to the origin, while the second behaves as z -n/2 . It follows that the pair Y n 2 √ Ez for φ and Y n 2 √ Ez for ψ is ruled out. The case of V = -n yields the same constraints due to the identities J -n = (-1) n J n and Y -n = (-1) n Y n .

No additional constraints stem from the behavior at infinity as all the functions J V (ρ), as well as Y V (ρ) (or J -V (ρ)), behave ∝ 1/ √ ρ with additional rapidly oscillating factors and therefore all the integrals converge.

Orthogonality.

The second point in the analysis of Eq. ( 5) is about the pairs of Bessel functions that ensure orthogonality. To that aim, it is convenient to recast Eq. ( 25) in the form

L ψE ≡ ψ E + 1 ρ ψ E - V 2 ρ 2 ψE = -E ψE . (26)
The two derivative terms do not combine to a total derivative because the factor of the first derivative term p 1 (ρ) ≡ 1/ρ is not the derivative of the factor p 0 (ρ) ≡ 1 of the second derivative term. However, the operator can be recast in a self-adjoint form by a proper choice of the weight function w (ρ) that appears in the definition of the scalar product w(ρ)f (ρ)g(ρ) dρ between two generic functions f (ρ) and g(ρ). The equation for w follows from integrating w(ρ)g(ρ)Lf (ρ) dρ by parts and imposing the self-adjoint condition (for the boundary terms, see next) or requesting (p 0 w) = p 1 w. Either way, we obtain w = w/ρ, that is w(ρ) = ρ. Note that the weight is proportional to the Jacobian for the change of variable

z → ρ = 2 √ z. It follows that integrals of the form f (2 √ z)g(2 √
z) dz that appear in Eq. ( 5)

for φ and ψ of the form

J V 2 √ Ez or Y V 2 √
Ez are naturally recast in terms of a weighted scalar product as

1 2
w(ρ)f (ρ)g(ρ) dρ.

We can use the above result for quantities of the type ρ φE1 (ρ) L ψE2 (ρ) dρ and obtain

(E 2 -E 1 ) ρ φE1 (ρ) ψE2 (ρ) dρ ∝ ρ φE1 ψ E2 -φ E1 ψE2 ∞ 0 . ( 27
)
At infinity, each term on the r.h.s. vanishes as all the functions J V (ρ) and Y V (ρ) (or J -V (ρ)) behave ∝ 1/ √ ρ.

Non-trivial conditions arise from the behavior at the origin, as discussed hereafter.

Let us first consider the "homologous" pair φE1 (ρ) = J V √ E 1 ρ and ψE2 (ρ) = J V √ E 2 ρ . Each boundary term on the r.h.s. of Eq. ( 27) scales then as ρρ V ρ V -1 ∝ ρ 2V , which is not vanishing for V ≤ 0. However, the two non-zero terms appearing on the r.h.s. are identical and therefore they cancel out. Subdominant terms then become relevant. Since the next subdominant term in the expansion of J V is quadratically smaller, we obtain that the r.h.s. in Eq. ( 27) behaves as ∝ ρ 2V +2 , which vanishes in the range of convergence of the integrals V > -1. The above argument yields obviously the same conclusion for the pairs φE1 (ρ) = J -V √ E 1 ρ and ψE2 (ρ) = J -V √ E 2 ρ , which are orthogonal for V < 1.

To conclude the analysis of homologous pairs, we need to consider integer values of V , for which J -V is not independent of J V , and Y V is the second independent solution.. For V = 0, the new case to consider is the pair φE1

(ρ) = Y 0 √ E 1 ρ and ψE2 (ρ) = Y 0 √ E 2 ρ
. The behavior of Y 0 at the origin is Y 0 (z) log (z/2). A simple calculation shows that the term in brackets on the r.h.s. of Eq. ( 27) is ∝ 1 2ρ log E1 E2 and therefore does not vanish for E 1 = E 2 . For positive integers V = n, the homologous pairs of J n are the same as above and the homologous pair of Y n does not ensure convergence of the integrals, as it was shown above. Finally, negative integers V = -n are like positive integers due to the identities J -n = (-1) n J n and Y -n = (-1) n Y n .

Let us now consider "heterologous" pairs. We start with non-integer V 's and consider pairs of the form φE1 (ρ) = J -V √ E 1 ρ and ψE2 (ρ) = J V √ E 2 ρ (or with the sign of V inverted). The explicit calculation of the boundary term of Eq. ( 27) shows that the r.h.s. is ∝ 2V (E 2 /E 1 )

V /2 and therefore does not vanish.

To conclude the analysis, we need to consider again integer values of V . For V = 0, the pair of eigenvectors is φE1 (ρ) = Y 0 √ E 1 ρ and ψE2 (ρ) = J 0 √ E 2 ρ (or J and Y inverted). By using Y 0 (z) log (z/2) and J 0 (z) z 0 , we can check that the boundary term in Eq. ( 27) does not vanish. For positive integers V = n, the pair of eigenvectors is φE1 (ρ) = Y n √ E 1 ρ and ψE2 (ρ) = J n √ E 2 ρ (or J and Y inverted). By using Y n (z) z -n and J n (z) z n , we can check that the boundary term in Eq. ( 27) does not vanish. Finally, negative integers V = -n are like positive integers due to the aforementioned identities J -n = (-1) n J n and Y -n = (-1) n Y n .

We summarize the conclusions in table I.

-1 < V < 0 z -V /2 JV 2 √ Ez z V /2 JV 2 √ Ez dEe -Et z -V /2 0 JV 2 √ Ez0 z V /2 JV 2 √ Ez z -V /2 J-V 2 √ Ez z V /2 J-V 2 √ Ez dEe -Et z -V /2 0 J-V 2 √ Ez0 z V /2 J-V 2 √ Ez 0 < V < 1 z -V /2 JV 2 √ Ez z V /2 JV 2 √ Ez dEe -Et z -V /2 0 JV 2 √ Ez0 z V /2 JV 2 √ Ez z -V /2 J-V 2 √ Ez z V /2 J-V 2 √ Ez dEe -Et z -V /2 0 J-V 2 √ Ez0 z V /2 J-V 2 √ Ez V = 0 J0 2 √ Ez J0 2 √ Ez dEe -Et J0 2 √ Ez0 J0 2 √ Ez V ≥ 1 z -V /2 JV 2 √ Ez z V /2 JV 2 √ Ez dEe -Et z -V /2 0 JV 2 √ Ez0 z V /2 JV 2 √ Ez V ≤ -1 z -V /2 J-V 2 √ Ez z V /2 J-V 2 √ Ez dEe -Et z -V /2 0 J-V 2 √ Ez0 z V /2 J-V 2 √ Ez

Normalization

The last item to be determined is the normalization conditions, i.e., the possible factors needed to normalize the integrals in Eq. ( 5) to a δ.

For integer values of V , one can directly use the formula 6.512.8 in [START_REF] Gradshteyn | Tables of Integrals, Series and Products[END_REF], which we report here

ρJ n (aρ) J n (bρ) dρ = δ (b -a) a . (28)
Applying the formula to our case with a = √ E 1 and b = √ E 2 , and keeping in mind the factor 1/2 of the Jacobian for the change of variable ρ = 2 √ z, give that the functions in table I are normalized without any additional factor. For non-integer values of V , the same conclusion holds, as demonstrated by using the formula 6.633.2 in [START_REF] Gradshteyn | Tables of Integrals, Series and Products[END_REF] ρe

-2 ρ 2 J ν (aρ) J ν (bρ) dρ = 1 2 2 e -a 2 +b 2 4 2 I ν ab 2 2 , (29)
valid for a > 0, b > 0, > 0 and ν > -1. We apply this formula with a = √ E 1 , b = √ E 2 and take the limit → 0. By using the asymptotic expansion I ν (z)

e z √ 2πz , we obtain ρ 2 J ν E 1 ρ J ν E 2 ρ dρ = lim →0 1 2 (E 1 E 2 ) 1/4 1 √ 4π 2 e -( √ E 1 - √ E 2 ) 2 4 2 = 1 2 √ E 1 δ E 1 -E 2 = δ (E 1 -E 2 ) . (30)
C. Explicit forms of the propagator.

The integrals over the energies E that appear in table I can be performed by using again the formula used in Eq. (29). Eq. ( 17) is recovered as

z z 0 V /2 dEe -Et J V 2 Ez 0 J V 2 √ Ez = 2 z z 0 V /2 ue -u 2 t J V (2u √ z 0 ) J V 2u √ z du = e -z+z 0 t t z z 0 V /2 I V 2 √ zz 0 t . (31)
"unusual" is with respect to the classical heat equation as similar limits appear in the footnote 6 in Ref. [START_REF] Gawȩdzki | Phase transition in passive scalar advection[END_REF] for the case of the compressible Kraichnan's model. We first provide some intuition by using the flux J, which vanishes for Eq. ( 31) and takes the value in Eq. (34) for Eq. (32). A term proportional to a is obtained by taking a derivative in (36) (the next term ∝ b scales as z and is therefore subdominant with respect to z V ) and compensating the remaining powers of z and t, i.e. z t -V z∂ z Ψ = aV t . In sum, by fixing the value µ at the origin z = 0 in the relation

J = µ lim z→0 t 2 zz 0 V z∂ z Q , (37)
we fix the ratio b/a = µ Γ(1-V ) Γ(1+V ) and span the various linear combinations. The no-flux solution with b = 0 is obviously µ = 0.

Let us now explicitly consider the behavior (36) of the linear combination itself close to the origin. For V > 0, the second term is dominant around the origin and lim z→0 Ψ(z) = b t V -1 z V 0 e -z 0 /t Γ(1-V ) . Vice versa, for V < 0, we have the first term dominating and lim z→0 z -V Ψ(z) = a 1 t 1+V e -z 0 /t Γ(1+V ) . Let us then consider the derivative. For V > 0, the dominant constant term drops out and we can obtain lim z→0 z 1-V ∂Ψ(z) = aV t 1+V e -z 0 /t Γ(1+V ) . For V < 0, we want to cancel out the dominant term, which requires the multiplication of Ψ by z -V . We obtain then lim z→0 z 1+V ∂ z -V Ψ(z) = b|V |t V -1 z V 0 e -z 0 /t Γ(1-V ) . In summary, the one-parameter family of solutions is obtained by prescribing the following quantities :

V > 0 b a = z V 0 t 2V V Γ (1 -V ) Γ (1 + V ) lim z→0 Ψ(z) lim z→0 z 1-V ∂Ψ(z) ; V < 0 a b = t 2V z V 0 |V |Γ (1 + V ) Γ (1 -V ) lim z→0 z -V Ψ(z) lim z→0 z 1+V ∂ (z -V Ψ(z)) , (38)
which is the analogous to the expression in footnote 6 of Ref. [START_REF] Gawȩdzki | Phase transition in passive scalar advection[END_REF] for the compressible Kraichnan model.

In conclusion, inspired by K. Gawȩdzki's work in Ref. [START_REF] Gawȩdzki | Phase transition in passive scalar advection[END_REF], we introduced an alternative, possibly more intuitive approach, for the selection of solutions and boundary conditions in models of turbulent transport. The basic idea is to impose the request of a bi-orthogonal expansion of the propagator. Integrability and orthogonality yield conditions that are analogous to those obtained via the methods in [START_REF] Reed | Methods of Modern Mathematical Analysis, II. Fourier Analysis, Self-Adjointness[END_REF]. This was first established for the well-studied example [START_REF] Reed | Methods of Modern Mathematical Analysis, II. Fourier Analysis, Self-Adjointness[END_REF] that was used for the compressible Kraichnan model originally discussed in Ref. [START_REF] Gawȩdzki | Phase transition in passive scalar advection[END_REF], and then shown for a general model of transport of odors in the atmospheric boundary layer. We wish we could have had the chance of discussing this alternative method and our results with Krzysztof Gawȩdzki, whom we remember with deep emotion.

TABLE I

 I

			. Possible solutions and forms of the propagator
	V	φE(z)	ψE(z)	P(z, z0 ; t)

Acknowledgements. This work was partially supported (MV) by the NIH Grant 1RF1NS128865-01.

It is verified that the corresponding flux, [V -z∇] •, vanishes at the origin z = 0. Conservation of normalization was already shown by Eq. ( 21).

The alternative possible expression is made explicit by the same formula Eq. (29) as

This alternative expression, which is the other independent solution to the Bessel Eq. ( 20), is a non-zero flux solution and does not conserve normalization. Indeed, the integral is

Here, the first equality is formula 6.643.2 in [START_REF] Gradshteyn | Tables of Integrals, Series and Products[END_REF] whilst the second is the definition 9.220.2 of the Whittaker function M in terms of the confluent hypergeometric Φ = 1 F 1 . The third equality uses the integral representation 9.211.2 of Φ. The incomplete Γ function is defined as Γ(a, z) ≡ ∞ z t a-1 e -t dt. The mass M (t) at any z 0 = 0 is then unity at t = 0 and eventually vanishes as t → ∞. The corresponding probability goes into a δ function at the origin, where particles are absorbed and progressively pile up. This is also seen more intuitively by noting that the flux J at the origin reduces to V times the value of the function at the origin, i.e.,

Integrating the flux in time, one gets the change of mass, which coincides with (33).

While the analysis of the dynamics in the presence of a finite diffusivity is beyond the scope of the present work, we make a short final remark on that case. The solution with a non-vanishing molecular diffusivity κ, i.e., to the equation

is obtained by noting that κ + z ≡ ζ reduces the equation in ζ to the same above form [START_REF] Murlis | Odor Plumes and How Insects Use Them[END_REF]. It is therefore sufficient to replace z by z + κ and z 0 by z 0 + κ in (31) and (32) to obtain solutions in the presence of a finite molecular diffusivity.

The role of κ is to regularize the behavior at the origin, as it can be explicitly verified by a direct expansion at the origin of those solutions.

IV. MIXED BOUNDARY CONDITIONS AND THEIR SELECTION

In the intermediate range -1 < V < 1 (and V = 0), table (I) shows that two solutions are possible and the propagator can be a linear combination of them, that is of the expressions (31) and (32) with respective coefficients a and b. The resulting dependence of the linear combination at small z (which we indicate by Ψ(z)) is

The sum of the coefficients a and b is fixed by normalization to unity at the initial time whilst their difference is free and yields a one-parameter family. In analogy with the case of Robin's conditions for the heat equation [START_REF] Gustafson | The third boundary condition -was it Robin's?[END_REF], we want to express the ratio b/a as a function of limiting values of the propagator and its derivative at the origin. Two major differences with respect to Robin's conditions appear, though : first, the function and/or its derivative can be singular at the origin ; second, contrary to the case of the heat equation, it is not the case that one component of the propagator vanishes at the origin while the other one has a vanishing derivative. It is not possible then to identify the two coefficients by simply taking the value of the propagator and the value of its derivative at the origin. The combination of these two aspects entails the "unusual" limits that will appear below. It is worth noticing that