CNRS David Nowak

Vlad Rusu

While Loops in Coq

While loops are present in virtually all imperative programming languages. They are important both for practical reasons (performing a number of iterations not known in advance) and theoretical reasons (achieving Turing completeness). In this paper we propose an approach for incorporating while loops in an imperative language shallowly embedded in the Coq proof assistant. The main difficulty is that proving the termination of while loops is nontrivial, or impossible in the case of nontermination, whereas Coq only accepts programs endowed with termination proofs. Our solution is based on a new, general method for defining possibly non-terminating recursive functions in Coq. We illustrate the approach by proving termination and partial correctness of a program on linked lists.

Introduction

The definition of recursive functions in the Coq proof assistant [START_REF] Bertot | Interactive Theorem Proving and Program Development -Coq'Art: The Calculus of Inductive Constructions[END_REF] is subject to certain restrictions to ensure their termination, which is essential for the consistency of Coq's underlying logic. Specifically, recursive calls must be made on strict subterms, effectively ensuring that the computation eventually reaches a base case. Alternatively, users have the option to prove that a specific quantity strictly decreases according to a well-founded order. In such cases, Coq can automatically transform the recursive calls into strict subterm calls, using a so-called accessibility proof to guarantee termination. Adhering to these constraints eliminates the risk of infinitely many calls, thereby ensuring that functions terminate.

An alternative, somewhat ad-hoc strategy is to introduce an additional natural-number argument called the fuel. The fuel's value is decremented with each recursive call, thereby guaranteeing finitely many recursive calls, hence, termination. However, a crucial concern arises as one must supply enough fuel so that termination does not disrupt the intended computation of the program by occurring too early.

In this paper we present a novel approach to defining possibly partial recursive functions in Coq while achieving separation of concerns: write the program first, and prove its properties (including termination) later. In broad terms our technique consists in providing an infinite amount of fuel for recursive functions. By doing so the function can proceed with its computations without risk of exhausting its fuel.

As a result, this approach empowers developers to focus on the core logic of the recursive function, separate from the termination concern, streamlining the development process and enhancing the modularity and readability of the code.

A key property guaranteed by our technique is that, given the functional of the recursive function under definition (i.e., an abstract description of the function's body), the resulting function is the least fixpoint of its functional. We prove this general result under mild constraints on the functional -it must be monotonic and, in some sense described precisely in the paper, must preserve continuity.

The method is applied to while loops in an imperative shallowly embedded in Coq. By proving that the functional of while loops is monotonic and continuity-preserving we obtain while loops as least fixpoints of their functionals. This enables programmers to construct imperative programs featuring arbitrary while loops and provides them with tools for reasoning about loops. Specifically, the leastfixpoint property is used for proving that a while loop terminates if and only if there exists some finite amount of fuel for which it returns the desired result; hence, a termination proof can proceed by induction on the fuel value, once an adequate instantiation for this known-to-exists value is chosen.

Subsequently, we proceed to establish a Hoare logic system, which serves as a formal tool for proving the partial correctness of programs. In essence the Hoare logic provides a systematic and rigorous approach to program verification, where one defines preconditions and postconditions that govern the state of the program before and after its execution. Through these assertions one can verify that the program's execution leads to the desired outcomes, establishing its partial correctness. Here, again, the property of being a fixpoint is used for proving the soundness of the while loop's Hoare triple.

Finally, the least-fixpoint property of the partial functions being defined ensures that the functions, defined abstractly using order theory, are, as mathematical functions, the same as the ones that Coq would have generated, had it not been constrained by its logic into rejecting the functions as partial.

Outline In Section 2 we introduce the reader, termination, and state monads, which serve the purpose of writing imperative programs in the Coq proof assistant. Moving forward to Section 3 we present our method for defining possibly nonterminating recursive functions and demonstrate its application to the definition of while loops. In Section 4 we define a monadic Hoare logic and illustrate its effectiveness by applying it to a program that computes the length of a linked list. In Section 5 we address the issue of proving termination. We compare with related work in Section 6 and conclude in Section 7.

The Coq development corresponding to this paper is available at https://tinyurl.com/2p93uwdj.

Monads for Possibly Nonterminating Stateful Computation

We consider a subset of Gallina (the programming language of Coq) expressive enough for shallowly embedding possibly nonterminating imperative programs. In purely functional languages such as Gallina the usual approach is to encode imperative features with monads [START_REF] Moggi | Computational Lambda-Calculus and Monads[END_REF]. We use a combination of the termination, state, and reader monads. The first one is used for possibly nonterminating computations, the second one for stateful computations where the state may change, and the third one for state-aware computations that do not change the state, such as checking the condition of while loops.

A monad consists of a type constructor equipped with an operation usually called ret for trivial computation, and another one, usually called bind, for sequencing computations. Each particular monad also comes with specific operations. Assume a context where a type T is declared. The type constructor for the reader monad is defined as follows:

Definition reader (A: Type): Type := T -> A.
where A is a type for values returned by a computation. Intuitively, reader is a side-effect free function from states to computed values. Unlike the state monad, shown below, it only reads from states, without modifying them; hence the monad's name reader.

Trivial computation ret consists in ignoring the state:

Definition ret {A: Type}(a: A) : reader A := fun _ => a.

Placing the parameter A between curly brackets marks it as implicit -its value can therafter be automatically inferred, relieving users from the burden of having to instantiate it. Sequencing of computations consists in passing the result of the first computation to the second one:

Definition rbind {A B: Type}(m: reader A)(f: A -> reader B): reader B := fun s => f (m s) s.

Finally, the reader monad is equipped with a primitive for reading the state:

Definition get: reader T := fun s => s.

The other monad used in the shallow embedding of an imperative language in Coq is a combination of termination and state monads, for which we use the option type constructor:

Inductive option (X: Type): Type := None: option X | Some: X -> option X.

where None encodes nontermination. Termination state monads are the programs of our imperative language, which is reflected in the name of its type constructor:

Definition program (A: Type): Type := T -> option (A * T).

The reader monad is a special case of the termination state monad. Thus we introduce a coercion that enables Coq to automatically convert a reader monad into a termination state monad when needed:

Coercion reader_to_program {A: Type}(m: reader A): program A := fun s => Some (m s, s).

The ret and the get primitive of the reader monad are thus automatically converted. But the sequencing of computations needs to be redefined in order the take into account the fact that the first of the sequenced computations might change the state, or might not terminate: The termination state monad comes with an additional primitive to change the state: Definition put (s: T) : program unit := fun _ => Some (tt, s).

Definition
The unit type (inhabited by one term -the constant tt) is used in our functional-laaguage setting for modelling imperative programs that do not return anything, encoded by returning the dummy constant tt.

As a running example we consider an imperative program that computes the length of a linked list. First we need to specify the State on which the program works. We assume a machine with two positive integer registers and an unbounded memory whose adresses and stored values are also positive integers:

Record State: Type := {reg1: nat; reg2: nat; memory: nat -> nat}.

Then, using the primitives of our monad we write operations to read/write in registers and memory:

• read_reg1, write_reg1 for reading, resp. writing in the first register, and similar functions for reading/writing in the second register. The second register is increased by one using incr_reg2;

• do next <-read_addr curr assigns to next the content of the address curr in memory. More generally, reading operations return a value that can be bound to an identifier using the standard "do" notation of monads, i.e., do x <-m; f x is a shortcut for rbind m (fun x => f x);

• operations are sequenced by double semicolons: m;;f is a notation for bind m (fun _ => f). It is assumed that the linked list1 of interest starts in memory at address addr. This address is written into the first register, then the second register (which is to contain the length of the list) is initialized to zero. Next, while the current address of the first register is not null (also encoded by zero), the second register is incremented and the first register is updated to point to the next element of the linked list. Finally, at the end of the while loop (if the end is ever reached), the value in the second register is the length of the list of interest, hence, it is the value returned by our function whenever it terminates.

What is still missing is, of course, the definition of the while loop. A first attempt uses recursion: That is, a recursive function (introduced by the keyword Fixpoint) attempts to define a while loop (with condition cond and body body of appropriate types) by first checking the condition, and if the condition holds, executing the body then recursively the while loop; otherwise, doing nothing (which is encoded by reader_to_program (ret tt)). However, this function does not always terminate. For example, if the while loop is used to navigate a linked list, like in the case of the length function above, and the list is badly linked, i.e., it contains a loop, then the while loop does not terminate. Coq rejects this definition attempt, as it rejects any recursive function whose termination it cannot infer.

In the rest of the paper we show how possibly infinite while loops (and, in general, partial recursive function) can be accepted by Coq by encoding nontermination as evaluation to a special value.

Partial Recursive Functions in Coq

The while loop is a particular case of a partial recursive function. We first sketch how partial recursive functions can be encoded in Coq before providing details.

Outline of the Approach

Assume we want to define a partial recursive function f from type A to type B. A natural way to proceed is to give f the type A -> option B, where for any a:A, f a = None encodes the fact that f is undefined for the input a. In order to define a function we further need its functional, an abstract representation of the body of the function being defined. Let F : (A -> option B)-> A -> option B be the functional for f. We say that f := F f is the fixpoint definition of f. The interesting case we here solve is when fixpoint definitions are not accepted by Coq -just like in the case of the while function above.

We proceed as follows. We define an auxiliary function f_fuel : nat -> A -> option B with an additional natural-number parameter called the fuel, as the following recursive function, which is accepted by Coq because Coq "sees" that the fuel parameter strictly decreases at each recursive call: If the functional F is monotonic then, based on results in order theory explained later in this section, the function f_fuel can be lifted to a continuous function in conat -> option B where conat is the type

Inductive conat: Type:= finite: nat -> conat | infinity
That is, the inhabitants of conat are natural numbers wrapped with the finite constructor, together with the constant infinity. Putting back the parameter a:A in the type we obtain a function f_inf of the type conat -> A -> option B. The results later in the section also ensure that, under an additional condition on F (preservation of continuity), the function (f_inf infinity) is the least fixpoint of F.

Recapitulating, we started with the intention of defining a function f : A -> option B, using its functional F : (A -> option B)-> A -> option B, via the fixpoint definition f := F f. We have assumed this is rejected by Coq. Per the results below we define f:= f_inf infinity and prove that is the least solution of the fixpoint equation f = F f -precisely the solution that Coq would have constructed had it accepted the definition f := F f -with the advantage that our definition is accepted.

Elements of Order Theory

The results in this subsection have been adapted from the textbook [START_REF] Roberto | Domains and lambda-calculi[END_REF]. We have formalized them in Coq, hence, hereafter proofs are only sketched or omitted altogether. The examples are not only used for illustration purposes: they also serve as building blocks in our approach to partial recursive functions.

Definition 1 A pointed partial order (PPO) (S , ⪯, ⊥) is a partially ordered set (S , ⪯) together with a distinguished element ⊥ ∈ S such that for all s ∈ S , ⊥ ⪯ s.

Example 1 The triple (N, ≤, 0) consisting of natural numbers N, their usual order ≤, and the least natural number 0 form a PPO.

Example 2 For any set A, the triple (A ∪ {⊥}, ⪯, ⊥) with ⊥ A, and ⪯ being defined as the smallest relation on A ∪ {⊥} such that ⊥ ⪯ a and a ⪯ a for all a ∈ A, is a PPO called the flat PPO of A.

In Coq the flat PPO of a type A is encoded using the type option A where None plays the role of ⊥.

Definition 2 Given a PPO (S , ⪯, ⊥), a set S ′ ⊆ S is directed if S ∅ and for all x, y ∈ S there exists z ∈ S such that x, y ≤ z.

Example 3 Any nonempty set of natural numbers in the PPO (N, ≤, 0) is directed.

The above example is a consequence of the more general fact that any nonempty sequence, i.e., totally ordered subset of elements in a PPO, is directed. Indeed, directed sets are generalizations of sequences. Least upper bounds of directed sets are generalizations of limits of sequences.

Example 5 Consider the PPO (N ∪ {∞}, ≤, 0) with the order ≤ on natural numbers extended such that ∞ ≤ ∞ and n ≤ ∞ for all n ∈ N. Then, (N ∪ {∞}, ≤, 0) is a CPO. Indeed, in this totally ordered set all subsets T ⊆ N ∪ {∞} are directed, and lub T is either:

• the maximum of T , if it exists • ∞, if the maximum of T does not exist.
In Coq the set N ∪ {∞} shall be encoded as the type conat seen earlier in this section.

Example 6 In a PPO (A ∪ {⊥}, ⪯, ⊥), the least upper bound of a singleton {x} is x and the least upper bound of a pair {⊥, a} is a. Those are the only directed sets in this PPO; hence, (A ∪ {⊥}, ⪯, ⊥) is a CPO.

Informally, the notion of compactness below captures what it means for an element to be finite.

Definition 4

In a CPO (S , ⪯, ⊥), an element s • ∈ S is compact whenever for all directed sets T ⊆ S , if s • ⪯ lub T then there exists t ∈ T such that s • ⪯ t.

Example 7 In the CPO (N ∪ {∞}, ≤, 0) the compact elements are exactly the (finite) natural numbers.

Example 8 In the CPO (A ∪ {⊥}, ⪯, ⊥) all the elements are compact. Some CPOs are, in the following sense, completely determined by their compact elements:

Definition 5 A CPO (S , ⪯, ⊥) having the set S • ⊆ S of compacts is algebraic if for all s ∈ S , the set {s • ∈ S • | s • ⪯ s} is directed, and s = lub {s • ∈ S • | s • ⪯ s}. Example 9 The CPO (N ∪ {∞}, ≤, 0) is algebraic. Indeed, for all n ∈ N ∪ {∞}, the set of compacts (natural numbers) {m ∈ N | m ≤ n} is directed (as is any nonempty subset of N ∪ {∞}). Moreover, • if n = ∞, then the set {m ∈ N | m ≤ n} coincides with N, and lub N = ∞; • if n ∈ N, then lub {m ∈ N | m ≤ n} = n.
Example 10 The flat CPO (A ∪ {⊥}, ⪯, ⊥) is algebraic. Indeed:

• for all a ∈ A, the set {⊥, a} of compacts in the ⪯ relation with a is directed, and a = lub {⊥, a}

• the set {⊥} of compacts in the ⪯ relation with ⊥ is directed and lub {⊥} = ⊥.

We shall use the following notion, which relates a PPO to the compact elements of an algebraic CPO: Definition 6 Consider a PPO (S • , ⪯ • , ⊥ • ) and an algebraic CPO (T, ⪯, ⊥) whose set of compacts is T • . We say that (T, ⪯, ⊥) is an embedding of (S • , ⪯ • , ⊥ • ) if there exists an injection ι : S • → T such that:

• ι ⊥ • = ⊥; • ι is monotonic; • ι maps S • to T • , written ι S • = T • .
The embedding, including the injection involved in it, is denoted by ι : (S • , ⪯ • , ⊥ • ) → (T, ⪯, ⊥).

Example 11

The embedding κ : (N, ≤, 0) → (N ∪ {∞}, ≤, 0) is induced by the canonical inclusion κ of N into N ∪ {∞}. The embedding of (A ∪ {⊥}, ⪯, ⊥) into itself is also induced by the canonical inclusion.

Remark: not all embeddings are induced by canonical inclusions. In Coq, in general, they are not. This is because Coq is based on type theory, hence, one cannot just add a new element to a type; to do this one must create a new type and wrap the old type in a constructor (which, in Coq, is always an injection). An example of this is the representation of N ∪ {∞} as the type conat with a constructor finite : nat -> conat. With the appropriate order and bottom element, conat is an embedding of nat induced by the constructor finite. The only situations when canonical inclusion induces an embedding in Coq is when it coincides with the identity function, like in the embedding of (A ∪ {⊥}, ⪯, ⊥) into itself.

Definition 7 Assume an embedding ι : (S • , ⪯ • , ⊥ • ) → (T, ⪯, ⊥). We denote by ι -1 : T → S • the (unique) function such that ι -1 (ι s • ) = s • for all s • ∈ S • , and ι -1 t = ⊥ • for t ∈ T \ (ι S • ).
Hence ι -1 is the inverse of ι on the compacts ι S • of T , and elsewhere it is given the (arbitrary) value ⊥.

The next theorem is our main ingredient for defining and reasoning about partial recursive functions. It uses the following notion of continuity: Definition 8 Given two CPOs (T, ⪯, ⊥) and (T ′ , ⪯ ′ , ⊥ ′ ), a function f : T → T ′ is continuous if for any directed set S ⊆ T , its image ( f S ) ⊆ T ′ is directed, and f (lub S ) = lub( f S ).

Theorem 1 Assume two embeddings ι 1 : (S • 1 , ⪯ • 1 , ⊥ • 1 ) → (T 1 , ⪯ 1 , ⊥ 1 ) and ι 2 : (S • 2 , ⪯ • 2 , ⊥ • 2 ) → (T 2 , ⪯ 2 , ⊥ 2 ) and a monotonic function f • : S • 1 → S • 2 .
Then there exists a unique continuous function f :

T 1 → T 2 such that f = ι 2 • f • • ι -1
1 -where • is the standard notation for function composition. If the embeddings in Theorem 1 are canonical inclusions we have a simpler version of the above result:

Corollary 1 Assume two embeddings ι 1 : (S • 1 , ⪯ • 1 , ⊥ • 1 ) → (T 1 , ⪯ 1 , ⊥ 1 ) and ι 2 : (S • 2 , ⪯ • 2 , ⊥ • 2 ) → (T 2 , ⪯ 2 , ⊥ 2 )
where ι 1 , ι 2 are canonical inclusions. Then, for any any monotonic function f

• : S • 1 → S • 2 there exists a unique continuous function f : T 1 → T 2 such that for all s • ∈ S • 1 , f s • = f • s • .

Application to Partial Recursive Functions

We use the existence part of Theorem 1 in order to define partial recursive functions and the uniqueness part in order to prove that the defined functions are least fixpoints of their respective fixpoint equations. The method has been formalized in Coq; we sketch it below in mathematical notation. Assume that we want to define a partial function f : A → (B ∪ {⊥}). We have have at our disposal the functional F : (A → (B ∪ {⊥})) → A → (B ∪ {⊥}). The following assumptions on F are required:

• monotonicity: for all f, f ′ : A → (B∪ {⊥}), if for all a ∈ A, f a ⪯ f ′ a then for all a ∈ A, F f a ⪯ F f ′ a,
where ⪯ is the flat order on (B ∪ {⊥}).

• preservation of continuity: assume an arbitary function g : (N ∪ {∞}) → A → (B ∪ {⊥}). If, for each a ∈ A, the function λ n → g n a is continuous (as a function between (N ∪ {∞}) and (B ∪ {⊥}) organized as CPOs) then, for each a ′ ∈ A the function λ n → F (g n) a ′ is continuous as well.

The method proceeds as a series of steps, grounded in the results from the previous subsection:

1. A function f • : N → A → (B ∪ ⊥
) is recursively defined by the equations: for all a ∈ A and m ∈ N,

f • 0 a = ⊥ and f • (m+1) a = F ( f • m) a;
intuitively, for all m ∈ N, ( f • m) constitute approximations of the function that we want to define, constrained by the finite amount of fuel m ∈ N;

2. the monotonicity requirement on F ensures that, for all a ∈ A, the function f • a = λ n → ( f • n a) is monotonic as a function between N and (B ∪ {⊥}) organized as PPOs;

3. the existence result of Theorem 1 ensures that, for all a ∈ A that there exists a continuous function f a : (N ∪ {∞}) → (B ∪ {⊥}), satisfying f a m = f • m a for all m ∈ N; here, we have used the fact that, in the embeddings involved in our application of Theorem 1, the injections are the canonical inclusions (cf. Example 11), hence, one can apply the simpler version of the theorem -Corollary 1;

4. using the uniqueness result of the corollary, for all a ∈ A, any continuous function

f ′ a : (N ∪ {∞}) → (B ∪ {⊥}) satisfying f ′ a m = f • (m + 1) a for all m ∈ N, also satisfies f ′ a = f a • (λ n → n + 1)
; here we have used the fact that the function λ n → n + 1 : (N ∪ {∞}) → (N ∪ {∞}) is continuous, and that the composition of continuous functions is a continuous function as well;

for all a

∈ A, let f ′ a : (N ∪ {∞}) → (B ∪ {⊥}) be defined by f ′ a = λ n → F (λ (x : A) → f x n) a.
Using the continuity preservation requirement on F, the continuity of f ′ a reduces to the continuity of

λ n → (λ (x : A) → f x n) a = λ n → f a n = f a ; since f a is continuous, f ′
a is continuous as well; 6. moreover, using the definitions of f • (first item in this list), of f a (item 2), and of f ′ a (item 5): for all m ∈ N,

f ′ a m = F (λx → f x m) a = F(λ x → f • m x) a = F ( f • m) a = f • (m + 1) a; hence, (cf. item 4), f ′ a = f a • (λ n → n + 1)
. This implies that for all n ∈ (N ∪ {∞}), f a (n + 1) = F (λx → f x n) ; 7. let now f : A → (B ∪ {⊥}) defined, for all a ∈ A, by f a = f a ∞. Then, using item 6 and

∞ = ∞ + 1: for all a ∈ A, f a = f a ∞ = f a (∞ + 1) = F (λ x → f x ∞) a = F (λ x → f x) a = F f a;
that is, we have obtained the fixpoint equation f = F f . What remains to be proved is that f is its least solution;

8. for this, we inductively define a sequence of functions in A → (B∪{⊥}) by F 0 = λ x → ⊥ and, for all m ∈ N, F m+1 = F(F m ). Using the definition of f we prove the equality f = lub {F n | n ∈ N}, where the least upper bound is taken in the CPO of functions A → (B ∪ {⊥}) ordered pointwise. Finally, we use a result that says that if F is monotonic on a CPO (it is, in our case, by the monotonicity assumption) and lub

{F n | n ∈ N} is a fixpoint of F (it is, in our case, since lub {F n | n ∈ N} = f and f = F f ) then lub {F n | n ∈ N} is the least fixpoint of F; i.e., f is the least fixpoint of F.
A comparison between the proposed approach for defining partial recursive functions and the standard one based on Kleene's fixpoint theorem is discussed in Section 6 dedicated to related works.

Instantiation to While Loops

The results from the previous subsection are now instantiated to while loops. Recall from subsection 3.1 the failed attempt at defining while loops in Coq, and notice their type: In order to instantiate the method described in the previous subsection to while loops we first need to change their type to A->option B for appropriate A, B. Remembering that program T unit is defined as T->option(unit*T), once the implicit parameter T is chosen, the type of while becomes

(reader T bool) -> (program T unit) -> T -> option(unit*T)

In order to obtain a type of the form A->option B we uncurry the above type to which concludes our construction of while loops in Coq. The next step is to provide users with means to reason about programs that contain such loops. This is the object of the next two sections. They shall be using the two following facts, which are consequences of while being the least fixpoint of its functional:

• an unfolding lemma, which is just another form of the fixpoint equation: • a lemma stating that the while loop evaluates to Some x in a state if and only if there exists a fuel-constrained version of the loop that also evaluates to the same Some x in the same state:

Lemma while_iff_while_fuel{T:Type}: forall (c:reader T bool)(b:program T unit)(s:T)(x:unit*T), while c b s = Some x <-> exists (fuel:nat), while_fuel fuel c b s = Some x.

Since evaluation to Some x models termination, the lemma can also be read as"if a loop terminates, then it terminates in finitely many steps" -where the number of steps is upper-bounded by fuel.

Partial Correctness

In this section we define a monadic Hoare logic for partial correctness [START_REF] Hoare | An Axiomatic Basis for Computer Programming[END_REF][START_REF] Swierstra | A Hoare Logic for the State Monad[END_REF]. Roughly speaking, partial correctness expresses the fact that a program returns the right answer whenever it terminates. In this paper, a program is a monadic computation. Remembering that Prop is the Coq type for logical statements, one writes the Hoare triple {{P}} m {{Q}} for the proposition hoare_triple P m Q defined in Coq by Definition hoare_triple {T A : Type} (P : T -> Prop) (m : program T A) (Q : A -> T -> Prop) : Prop := forall s s' a, P s -> m s = Some (a, s') -> Q a s'.

That is, if the program m: (program T A) is in a state s: T such that the pre-condition P: T -> Prop holds for s, if the program terminates (encoded by the fact that the program returns Some (a, s')) then the pair (a, s') satisfies the postcondition Q: A -> T -> Prop.

There are Hoare triples for all monadic instructions, but the triple of interest is this paper is the one for the while loops. It states that: if the body of the loop preserves an invariant I as long as the condition cond of the loop is true, then the loop preserves the invariant whenever it terminates. The lemma while_fuel_triple is used in the proof of while_triple, together with the lemmas while_iff_while_fuel and while_unfold shown at the end of the previous subsection.

Then, while_triple is used to prove the weakest precondition triple for our running example:

Lemma length_wp (addr: nat)(P: nat -> State -> Prop): {{ fun s => forall len, Length s addr len -> P len {| reg1 := 0; reg2 := len; memory := s.(memory) |} }} length addr {{ P }}.

where Length is an inductively defined relation that says that in a state s, the list starting at address addr has length len. Having this predicate gives us an abstract manner of defining a linked list's length: In accordance to the laws of weakest precondition, length_wp says that the postcondition P must hold in the precondition when applied to the expected return value upon termination (the length len) and to the expected state upon termination {|reg1 := 0; reg2 := len; memory := s.(memory)|}.

As usual with Hoare logic, the crux of the proof is to find the right loop invariant to be fed to the lemma while_triple. Here, it is a generalization of the precondition in length_wp: With this choice of invariant the proof of length_wp is just a matter of unfolding definitions. The interest of having proved a weakest precondition lies in its generality. As immediate corollaries of length_wp we obtain a first lemma that states that whenever length addr terminates, the register reg2 contains the length of the linked list starting at address addr: Lemma length_correct1 (s0: State)(addr: nat) : {{ fun s => s = s0 }} length addr {{ fun _ s' => Length s0 addr s'.(reg2) }}.

And another lemma stating that: if the linked list starting at address addr has length len, then this is the value that will be returned by length addr whenever it terminates.

Lemma length_correct2 (len: nat)(addr: nat) : {{ fun s => Length s addr len }} length addr {{ fun n _ => n = len }}.

Termination

In our approach, termination is modeled by evaluation to Some value. In order to succesfully prove termination we need to specify quite precisely the value to which a program evaluates. The key for the termination of the length program is the termination of its while loop, which is expressed as follows:

Lemma while_terminates: forall len s addr, Length s addr len -> forall n, while do curr <-read_reg1; ret (curr != 0) (incr_reg2;; do curr <-read_reg1; do next <-read_addr curr; write_reg1 next) {| reg1 := addr; reg2 := n; memory := s.(memory) |} = Some (tt, {| reg1 := 0; reg2 := n + len; memory := s.(memory) |}).

That is, when called in a state where the first register reg1 points to the beginning of the list and the second register reg2 is initialized with some value n, the while loop ends in a state where the first register is null and the second register contains n plus the length of the list. The memory field of the state remains unchanged because the loop does not write in it. The return value tt is the unique inhabitant of unit and encodes the fact that while loops do not actually return anything relevant.

Lemma while_terminates is proved by induction on n and uses the lemma while_unfold to unfold the loop in the inductive step. It is then used to prove the termination of the length program:

Lemma length_terminates: forall s len addr, Length s addr len -> length addr s = Some (len,{|reg1 := 0; reg2 := len; memory := s.(memory)|}).

This says that the function call length addr s terminates with value len whenever, according to the inductive relation Length, the state s has a memory field where there is a well-formed linked list of length len starting at address addr. If the list were not well-formed, i.e. its links would form a loop, then the function would evaluate to None, which, in our method denotes non-termination.

Related Work

In domain theory [START_REF] Winskel | The formal semantics of programming languages -an introduction[END_REF] partial recursive functions are typically defined as least fixpoints of their functionals using Kleene's fixpoint theorem, which states that a functional has a least fixpoint whenever it is continuous, and that the least fixpoint is obtained by infinitely many iterations of the functional starting from a function defined nowhere. This theorem is very elegant and easy to prove. Our own version of a fixpoint theorem (perhaps less elegant than Kleene's theorem, and definitely harder to prove) has the same conclusion but requires that the functional be monotonic and continuity-preserving, which roughly means that, given a certain continuous function, the functional produces a continuous function. From our own experience and that of other authors (e.g., [START_REF] Bertot | Fixed point semantics and partial recursion in Coq[END_REF]) it appears that using Kleene's theorem, which requires proving continuity, is difficult in practice, even in the simplest cases. For example, it took us hundreds of lines of Coq code just to prove the continuity of the successor function on natural numbers extended with infinity. By contrast, using our version of the theorem requires a proof of preservation of continuity, which appear to be more manageable -for while-loops in a shallow embedding of an imperative programming language in Coq the proof of continuity-preservation is remarkably simple: ten lines of Coq code. One possible reason for why continuity-preservation is easier to prove for functionals than continuity is that that the latter refers to the higher-order functional itself, whereas the former concerns the argument of the functional, which is a simpler function, one order below the order of the functional.

Other authors have explored partial recursive functions in Coq. In [START_REF] Capretta | General recursion via coinductive types[END_REF] a partial recursive function's codomain is a thunk -a parameterized coinductive type that "promises" an answer as a value of its parameter, but may postpone this answer forever, yielding nontermination. However, the functions being defined now become corecursive functions, which are restricted in the Coq proof assistant. As a result, only tail-recursive functions can be defined with this approach. This was also noted in [START_REF] Chlipala | Certified Programming with Dependent Types -A Pragmatic Introduction to the Coq Proof Assistant[END_REF]Chapter 7.3].

The same author [6, Chapter 7.2] proposes an alternative for the codomain of a partial function: a computation is a type that associates to a natural-number approximation level an approximation of the intended function. Like us the author uses the Coq option type, where None stands for nontermination and Some for termination with a return value; and a monad of computations is designed so that computations can encode imperative features. However, [6, Chapter 7.2] requires that the functional of the function being defined be continuous, for a definition of continuity equivalent to that employed in domain theory; hence this approach is subject to the general difficulty of proving continuity of functionals.

Regarding the embedding of imperative languages in proof assistants for the purpose of program verification, an alternative to the shallow embedding that we here use is deep embedding, which consist in defining the syntax, operational semantics [START_REF] Krebbers | Formal C Semantics: CompCert and the C Standard[END_REF][START_REF] Leroy | Formal verification of a realistic compiler[END_REF], and program logic for the guest language in the logic of the host. Among the many projects we here cite the Iris project [START_REF]The Iris Project[END_REF], a rich environment built around a concurrent programming language and the corresponding program logic -concurrent separation logic.

Conclusion and Future Work

Recursive functions in Coq need to terminate for the underlying logic to be sound. Coq typically ensures termination via an automatically-checkable structural decreasing of terms to which recursive calls apply. In more complex cases Coq can be helped by the user with termination proofs. The termination proof becomes a part of a function definition; a function is not defined until the termination proof is completed.

For various reasons falling under the general notion of separation of concerns it is desirable to separate function definitions from termination proofs. It is also useful to have functions that do not terminate on some inputs. This paper proposes a new approach that achieves these desirable features. Non-termination is simulated as evaluation to a special value interpreted as "undefined". Under mild conditions on the function's body encoded as a higher-order functional, a possibly non-terminating function is defined and proved to be the least fixpoint of its functional according to a certain definition order.

We instantiate the general approach to while-loops in an imperative language shallowly embedded in Coq. The shallow embedding is based on a combination of monads. The least-fixpoint property of the resulting while loops is a key property enabling termination and partial-correctness proofs on imperative programs containing them. The practicality of the approach is illustrated by proving partial correctness and termination properties on a program computing the length of linked lists. Partial correctness is expressed in Hoare logic and is proved in the standard manner, by having users provide a strong-enough invariant; and termination is proved by having users provide an upper bound for the number of iterations.

Future Work A promising line of future work is to extend our approach to defining partial corecursive function in Coq. The idea is that codomains of such functions would be encodings of coinductive types organized as algebraic CPOs, generalizing the option types that we here used for recursive functions. Initial experiments with defining some difficult corecursive functions that go beyond Coq's builtin corecursion mechanisms (a filter function on streams, a mirror function on Rose trees) are promising.

The function-definition mechanism presented in this paper critically depends on functionals preserving continuity. A deeper understanding of the relationship between continuity-preservation and continuity, and of the perimeter where the continuity-preservation property holds, is also left for future work.

A more practical future work direction is to apply the instance of our approach to imperative-program definition and verification. Our intention is to build on our experience with proving low-level programs manipulating linked lists [START_REF] Jomaa | Proof-Oriented Design of a Separation Kernel with Minimal Trusted Computing Base[END_REF]. An interesting logic to consider in this setting is separation logic [START_REF] Reynolds | An Overview of Separation Logic[END_REF], perhaps by drawing inspiration from the shallow embedding of separation logic in Coq presented in [START_REF] Sergey | Programs and Proofs: Mechanizing Mathematics with Dependent Types[END_REF].

  bind {A B: Type}(m: program A)(f: A -> program B): program B := fun s => match m s with None => None | Some (a, s') => f a s' end.

Fixpoint while{T:

  Type}(cond: reader T bool)(body: program T unit): program T unit:= do c <-cond;if c then body;;while cond body else reader_to_program (ret tt)

Fixpoint

  f_fuel (fuel: nat) (a: A) : option B:= match fuel with |S fuel' => F(f_fuel fuel') a (*S is the successor function on natural numbers*) |0 => None end

Example 4

 4 In a flat PPO (A ∪ {⊥}, ⪯, ⊥), the directed sets are exactly: the singletons {x} with x ∈ A ∪ {⊥}, and the pairs of elements of the form {a, ⊥} with a ∈ A.Definition 3 A Complete Partial Order (CPO) is a PPO (S , ⪯, ⊥) with the additional property that any directed set T ⊆ S has a least upper bound, denoted by lub T .

Fixpoint

  while{T:Type}(cond:reader T bool)(body:program T unit):program T unit:= do c <-cond; if c then body ;; while cond body else reader_to_program (ret tt)

(

  (reader T bool)*(program T unit)*T) -> option(unit*T) where * builds products between types. Next, we define a function while':A->option B where A=(reader T bool)*(program T unit)*T and B=unit*T. For this we first write the functional for the while' function as follows Definition While'{T:Type}(W:((reader T bool)*(program T unit)*T)->option unit*T) (p:(reader T bool)*(program T unit)*T): option unit*T := let (cond,body,s) := decompose p in ( do c <-cond; if c then body;;(fun (s':T)=>W (cond, body, s')) (*after ;; a function on T is expected*) else reader_to_program (ret tt) ) s (Notice how the parameter p was decomposed into three components.) After proving that While' is monotonic and preserves continuity, we obtain using the method in Subsection 3.3 the function while': (reader T bool)*(program T unit)*T -> option unit*T as the least fixpoint of While'. What remains to be done is to curry the type (reader T bool)*(program T unit)*T->option unit*T to the expected type of the while function. When this is done, we obtain while{T:Type}: (reader T bool)->(program T unit)->T->option unit*T as the least fixpoint of the functional Definition While{T:Type}(W:(reader T bool)->(program T unit)->T->option unit*T) (cond :(reader T bool))(body : (program T unit))(s :T)) : option unit*T := (do c <-cond; if c then body;;(W cond body) else reader_to_program (ret tt)) s

  Lemma while_unfold{T:Type}: forall(c: reader T bool)(b: program T unit), while c b = (do c' <-c; if c' then b;;while c b else reader_to_program (ret tt))

  Lemma while_triple {T: Type}(cond: reader T bool)(body: program T unit)(I: T -> Prop): {{ fun s => I s /\ cond s = true }} body {{ fun _ s' => I s' }} -> {{ I }} while cond body {{ fun _ s' => cond s' = false /\ I s' }}. In order to prove while_triple we first prove a triple for fuel-constrained loops by induction on fuel: Lemma while_fuel_triple {T:Type}(fuel:nat)(cond: reader T bool)(body: program T unit)(I: T -> Prop): {{ fun s => I s /\ cond s = true }}body{{ fun _ s' => I s' }} -> {{ I }} while_fuel fuel cond body {{ fun _ s' => cond s' = false /\ I s' }}.

  Inductive Length (s:State): nat -> nat -> Prop := | length_nil : forall addr, addr = 0 -> Length s addr 0 | length_cons : forall addr len, addr <> 0 -> Length s (s.(memory) addr) len -> Length s addr (S len).

  fun _ s => forall len, Length s s.(reg1) len -> P (len+s.(reg2)) {| reg1:= 0; reg2:= len + s.(reg2) ; memory:= s.(memory) |}).

  This almost allows us to write a program computing the length of a linked list. Not completely, because the program uses a while loop, which is not defined until later in the paper:

	Definition length (addr: nat): program State nat :=
	write_reg1 addr;;
	write_reg2 0;;
	while (do curr <-read_reg1; ret (curr != 0))
	(
	incr_reg2;;
	do curr <-read_reg1;
	do next <-read_addr curr;
	write_reg1 next
	);;
	do res <-read_reg2; ret res.

For simplicity we consider linked lists where each element only contains the next element's address.


Rabisankar Samanta

Romain Pierrat
email: romain.pierrat@espci.psl.eu

Rémi Carminati

Sushil Mujumdar
email: mujumdar@tifr.res.in

Photon diffusion in space and time in a second-order nonlinear disordered medium

We report experimental and theoretical investigations on photon diffusion in a second-order nonlinear disordered medium under conditions of strong nonlinearity. Experimentally, photons at the fundamental wavelength (λ = 1064 nm) are launched into the structure in the form of a cylindrical pellet, and the second-harmonic (λ = 532 nm) photons are temporally analyzed in transmission. For comparison, separate experiments are carried out with incident green light at λ = 532 nm. We observe that the second harmonic light peaks earlier compared to the incident green photons. Next, the sideways spatial scattering of the fundamental as well as second-harmonic photons is recorded. The spatial diffusion profiles of second-harmonic photons are seen to peak deeper inside the medium in comparison to both the fundamental and incident green photons. In order to give more physical insights into the experimental results, a theoretical model is derived from first principles. It is based on the coupling of transport equations. Solved numerically using a Monte Carlo algorithm and experimentally estimated transport parameters at both wavelengths, it gives excellent semi-quantitative agreement with the experiments for both fundamental and second-harmonic light.

I. INTRODUCTION

Electromagnetic wave propagation and scattering are ubiquitous in diverse fields such as optics, condensed matter physics, biology, atmospheric optics, etc [START_REF] Ishimaru | Wave Propagation and Scattering in Random Media[END_REF]. Among the various related phenomena, diffusion of light has attracted maximum attention in the last decades [START_REF] Carminati | Principles of Scattering and Transport of Light[END_REF]. Two factors have primarily motivated these studies, namely, the occurrence of diffusion of light in tissues [START_REF] Martelli | Light Propagation through Biological Tissue and Other Diffusive Media: Theory, Solutions, and Software[END_REF], and the parallels of light diffusion with electron propagation in disordered conductors [START_REF] Kittel | Introduction to Solid State Physics[END_REF]. The first sub-field, namely, diffusion in tissues, has led to significant signs of progress in imaging of inclusions in living media [START_REF] Tuchin | Tissue Optics: Light Scattering Methods and Instruments for Medical Diagnosis[END_REF], with an aim to replace hazardous high-energy radiation. The second domain overlaps with condensed matter studies and is intended to understand mesoscopic effects in light transport in disordered media. Indeed, studies on mesoscopic optics have led to better insights on electronic transport, such as the transition from the diffusion regime to localization regime [START_REF] Anderson | Absence of diffusion in certain random lattices[END_REF] upon a sufficient increase in disorder.

A primary reason for the success of these studies in the optical domain has been the inherent noninteracting nature of the photons, as well as the possibility to directly image the intensity distribution, which makes the analyses easier than for electrons. Interestingly, photon propagation allows for additional aspects of transport, such as amplification or nonlinearity. On the one hand, coherent amplification accompanied by diffusion or localization has led to the creation of novel optical sources named random lasers [START_REF] Letokhov | Generation of light by a scattering medium with negative resonance absorption[END_REF][START_REF] Lawandy | Laser action in strongly scattering media[END_REF][START_REF] Pradhan | Localization of light in coherently amplifying random media[END_REF][START_REF] Wiersma | The physics and applications of random lasers[END_REF][START_REF] Tiwari | Random lasing over gap states from a quasi-one-dimensional amplifying periodic-on-average random superlattice[END_REF][START_REF] Uppu | Exponentially tempered lévy sums in random lasers[END_REF]. On the other hand, studies coupling nonlinearity with the disorder have attracted growing attention [START_REF] Melnikov | Second-harmonic generation in strongly scattering porous gallium phosphide[END_REF][START_REF] De Boer | Transmission and reflection correlations of second harmonic waves in nonlinear random media[END_REF][START_REF] Qiao | Cavity-enhanced second-harmonic generation in strongly scattering nonlinear media[END_REF][START_REF] Fröhlich | Localization in disordered, nonlinear dynamical systems[END_REF][START_REF] Mckenna | Effects of nonlinearity on anderson localization[END_REF][START_REF] Schwartz | Transport and anderson localization in disordered twodimensional photonic lattices[END_REF][START_REF] Fishman | The nonlinear schrödinger equation with a random potential: results and puzzles[END_REF][START_REF] Lahini | Anderson localization and nonlinearity in one-dimensional disordered photonic lattices[END_REF][START_REF] Heiderich | Coherent backscattering in nonlinear media[END_REF][START_REF] Wellens | Coherent backscattering of light by nonlinear scatterers[END_REF][START_REF] Agranovich | Effects of weak localization of photons in nonlinear optics: Second harmonic generation[END_REF][START_REF] Yoo | Search for the effect of weak photon localization in secondharmonic waves generated in a disordered anisotropic nonlinear medium[END_REF][START_REF] Simon | Coherent backscattering of optical second-harmonic generation with long-range surface plasmons[END_REF][START_REF] Wang | Coherent backscattering of optical second-harmonic generation in silver films[END_REF]. For instance, a large effort has been focused on the effects of χ (3) nonlinearity on Anderson localization [START_REF] Fröhlich | Localization in disordered, nonlinear dynamical systems[END_REF][START_REF] Mckenna | Effects of nonlinearity on anderson localization[END_REF][START_REF] Schwartz | Transport and anderson localization in disordered twodimensional photonic lattices[END_REF], wherein the question of whether nonlinearity subjugates localization is still an open question [START_REF] Fishman | The nonlinear schrödinger equation with a random potential: results and puzzles[END_REF]. Experimental attempts to address this problem in one-dimensional systems showed that the localized mode continues to exist under introduced nonlinearity [START_REF] Lahini | Anderson localization and nonlinearity in one-dimensional disordered photonic lattices[END_REF]. Theoretically, the effects of χ (3) nonlinearity on coherent backscattering (or weak localization) have also been documented [START_REF] Heiderich | Coherent backscattering in nonlinear media[END_REF][START_REF] Wellens | Coherent backscattering of light by nonlinear scatterers[END_REF]. Likewise, χ (2) materials, in particular relevance to second harmonic generation (SHG) [START_REF] Franken | Generation of optical harmonics[END_REF], have been studied in the regimes of diffusion and weak localization [START_REF] Agranovich | Effects of weak localization of photons in nonlinear optics: Second harmonic generation[END_REF][START_REF] Yoo | Search for the effect of weak photon localization in secondharmonic waves generated in a disordered anisotropic nonlinear medium[END_REF][START_REF] Simon | Coherent backscattering of optical second-harmonic generation with long-range surface plasmons[END_REF][START_REF] Wang | Coherent backscattering of optical second-harmonic generation in silver films[END_REF]. Specifically, several consequences of random quasi-phase-matching in disordered χ (2) materials have been discussed in recent years [START_REF] Savo | Broadband mie driven random quasi-phasematching[END_REF][START_REF] Müller | Modeling of random quasi-phasematching in birefringent disordered media[END_REF][START_REF] Morandi | Multiple scattering and random quasi-phasematching in disordered assemblies of LiNbO3 nanocubes[END_REF]. Furthermore, speckle dynamics in these materials have been recently quantified and theoretically modeled [START_REF] Samanta | Intensity-dependent speckle contrast of second harmonic light in a nonlinear disordered medium[END_REF][START_REF] Samanta | Speckle decorrelation in fundamental and second-harmonic light scattered from nonlinear disorder[END_REF]. Interestingly, second harmonic generation in a diffusive material creates a unique paradoxical situation of incoherent transport, i.e. diffusion, coupled with an inherently coherent phenomenon of second harmonic generation, and hence is of fundamental research interest. Over the years, a few theoretical [START_REF] Kravtsov | Theory of second-harmonic generation in strongly scattering media[END_REF][START_REF] Makeev | Second harmonic generation in suspensions of spherical particles[END_REF] and experimental reports [START_REF] Faez | Experimental observation of secondharmonic generation and diffusion inside random media[END_REF] have addressed this peculiar scenario. In one of the earliest reports, Kravtsov et al. [START_REF] Kravtsov | Theory of second-harmonic generation in strongly scattering media[END_REF] theoretically studied two models of disorder, namely point scatterers in a nonlinear medium and grainy nonlinear scatterers, and observed sharp peaks in the angular distribution of backward diffuse second harmonic light. In another approach, Makeev et al. theoretically studied the diffusion of second-harmonic (SH) light in a colloidal suspension of spherical nonlinear particles, wherein they found that the average SH intensity was independent of the linear scattering properties of the medium [START_REF] Makeev | Second harmonic generation in suspensions of spherical particles[END_REF]. In the experimental effort on GaP powders, a consistent picture that described the second-harmonic intensity distribution in the sample was obtained via the diffusion equation, which invoked nonlinearity as a conversion rate [START_REF] Faez | Experimental observation of secondharmonic generation and diffusion inside random media[END_REF].

The various models described above successfully invoked the scenario of weak nonlinearity, wherein the nonlinearity did not affect the distribution of fundamental photons. However, these works focused on the spatial behavior of fundamental and second-harmonic light, while no attention was paid to the temporal behavior, under ultrashort pulse illumination. In this work, we address this question, and others, in the light of experiments and of a theoretical model. We carry out spatial and temporal investigations of light propagating through a strongly disordered KDP (potassium dihydrogen phosphate) pellet. We diagnose pulse propagation through the system and measure the transport mean free paths for the fundamental and second harmonic light using the transmitted pulses. Subsequently, we measure the longitudinal spatial distribution of both components. To give more physical insights into the processes at play, we build from first principles a theoretical model involving the coupling of two Radiative Transfer Equations (RTE). Solved numerically using a Monte Carlo scheme, the model gives reliable results compared to the experiment. In Sec. II, we introduce the samples and the experiments and we present all experimental results. In Sec. III, we present the theoretical model as well as the associated numerical simulations in the same geometry as in the experiments. The discussion and conclusion are given in Sec. IV.

II. SAMPLES AND EXPERIMENTS

We exploited commercially available KDP (@ EM-SURE ACS) as our nonlinear material. It has a large second-order optical nonlinearity (∼ 0.43 pm/V). A fine powder of KDP was realized after 10 minutes of ball milling. Figure 1(a) depicts the scanning electron micrograph image of KDP grains while the inset shows the size distribution of the KDP grains, with grain sizes in the range of 2 µm to 8 µm. The size distribution (bar graph, red) is approximately log-normal (yellow line), peaking at 3.11 µm. A pellet of length L = 5 mm and radius a = 0.5 cm was created under a hydraulic press. The pellet was baked to remove remnant moisture at a temperature of 80 °C, which was far below the transition temperature (190 °C, [START_REF] Itoh | X-ray study of hightemperature phase transitions in KH2P O4[END_REF]) of KDP. Before going to the main experiment, we evaluated the effective refractive indices and internal reflection coefficients to be used later as input parameters in our theoretical model. Using the Maxwell-Garnett formula [START_REF] Choy | Effective Medium Theory[END_REF] (see App. A regarding the use of the Marxell-Garnett formula), the real parts of the effective refractive indices of the pellet were found to be n r (2ω) = 1.45 and n r (ω) = 1.43 at λ = 532 nm and λ = 1064 nm respectively. The internal power reflection coefficients at λ = 532 nm and λ = 1064 nm are estimated to be R(2ω) = 0.53 and R(ω) = 0.52 respectively [START_REF] Zhu | Internal reflection of diffusive light in random media[END_REF] (see App. A for the detailed derivation of the internal reflection coefficient). A schematic of the experimental setups is depicted in Fig. 1 (b) and (c). Two different experiments were performed. For the temporal measurements [Fig. 1 (b)], Nd:YAG ps laser pulses (EKSPLA, PL2143B), with pulse width ∼ 30 ps at the fundamental wavelength of λ = 1064 nm (hereafter referred to as IR) and beam waist ∼ 5 mm, were made incident on the front face of the pellet. The second harmonic generated light (hereafter referred to as SHG) that was transmitted from the back face was directed into the Streak Camera (Optronics SC-10). The streak camera possesses a temporal resolution of less than 3 ps and its specified spectral sensitivity ranges from 200 to 950 nm. The spectral range of the Streak Camera did not allow for measuring the IR pulse. A small part of the second harmonic light generated by a clear nonlinear crystal was directed toward the marker pulse input of the Streak Camera to calibrate the zero of the time axis. For the spatial measurements [Fig. 1(c)], the pellet was pumped with the fundamental of the Nd: YAG ns laser (pulse width ∼ 5 ns, repetition rate 10 Hz). The pulsed IR beam was launched normally onto the front face of the pellet. The scattered IR photons were imaged from one side of the pellet by a combination of lens and InGaAs CCD, as shown in Fig. 1. Simultaneously, a Silicon CCD along with a lens was employed to measure the internally generated and scattered second harmonic photons (SHG). For comparison to be used later, in both the temporal and spatial experiments, we also externally launched frequency-doubled green light from the source laser (λ = 532 nm, hereafter referred to as incident green or IG), whose temporal transmitted profile and spatial scattering profile were measured separately by the same Streak Camera and Silicon CCD respectively. The CCD images provided the spatial variation of light intensity along the length of the pellet.

Our sample consists of a large number of KDP microcrystals with random orientations. The incident photons experience multiple scattering due to the refractive index mismatch between the microcrystals and the background medium. SHG photons are generated at random positions within the medium by the crossing of two IR photons and also undergo multiple scattering. At large optical thickness and long times, we can expect the transmitted multiply scattered IR and SHG intensities to follow a diffusion process, as confirmed by the spatial intensity profiles shown below. The experimental temporal profiles of SHG (green squares) and IG (blue dots) in the pellet, plotted in Fig. 2, exhibit a classic diffusive behavior. We observe an interesting trend in the transmission where SHG photons have shorter residence times than IG photons, given that the SHG profile peaks earlier. The temporal profiles of SHG and IG show clear exponentially decaying tails, as evidenced in the inset on a logarithmic scale. Interestingly, the decay rates for the two curves are also very similar, a fact which we will return to in our theoretical work. The tails were fit by the expression A exp(-Γt), with Γ as the decay rate. The estimated decay rates fitted to numerous temporal profiles of the same sample amount to 3.17 ± 0.13 ns -1 and 3.25 ± 0.11 ns -1 for the SHG and IG, respectively. These values are as close as a fitting routine can provide. The apparent deviation may therefore be probably due to the strong noise level. The orange dashed and red solid lines in Fig. 2 correspond to the fits for SHG and IG respectively. Nonetheless, we emphasize here that, only the IG data can be fit a priori using diffusion theory. Although we use the same analysis for the SHG light, the fact that it follows a similar diffusion equation cannot be justified based on an existing theory. The theory developed in the next section will elaborate on this part and also justify the comparable decays seen in the experiments. In the presence of absorption and side loss, the decay rate Γ in a cylindrical system can be written as [START_REF] Gaikwad | Photon transport in cylindricallyshaped disordered meso-macroporous materials[END_REF] Γ

(2ω) = D(2ω) α 2 1 a e (2ω) 2 + π 2 L e (2ω) 2 + 3 ℓ t (2ω)ℓ a (2ω)
(1) where the diffusion coefficient is given by

D(2ω) = v E (2ω)ℓ t (2ω) 3 . (2)
For non-resonant scatterers, the energy velocity can be well approximated by v E (2ω) = c/n r (2ω) where c is the light velocity in vacuum. ℓ t and ℓ a are the transport and absorption mean free paths respectively. a e and L e are the effective radius and thickness of the pellet respectively and are given by a e (2ω) = a + z e (2ω) and L e (2ω) = L+2z e (2ω) where z e is the extrapolation length given by

z e (2ω) = 2 3 ℓ t (2ω) 1 + R(2ω) 1 -R(2ω) . (3)
α 1 in Eq.1 is the first zero of the zero-order Bessel function J 0 . By fitting the experimentally measured decay rates at λ = 532 nm for different sample lengths with Eq. ( 1), we extract ℓ t (2ω) ≈ 32 µm and ℓ a (2ω) ≈ 11 mm (side loss). Due to the limitation of spectral sensitivity, the streak camera could not be used for the temporal profile of the IR light and the determination of the transport mean free path. We derived the value of ℓ t (ω) from a coherent backscattering experiment [START_REF] Maret | Multiple light scattering from disordered media: the effect of brownian motion of scatterers[END_REF]. The estimated value of ℓ t (ω) is ∼ 216 µm. Note that for both the harmonics, {L, a} ≫ ℓ t and t ≫ v E ℓ t . This signifies that our sample is clearly in the diffusive regime [START_REF] Zhang | Wave transport in random media: The ballistic to diffusive transition[END_REF].

Next, we examine the spatial diffusion profiles in the pellet. Figure 3 displays the experimentally measured spatial profiles for the IG (blue solid line), IR (orange dotted line), and SHG (green dashed line) photons. Each profile shows a characteristic diffusive peak inside the input edge. The first peak to appear is the IG one followed by the IR one. This behavior is expected since ℓ t (2ω) < ℓ t (ω). The SHG peak may be expected to appear at a location where the fundamental intensity is maximum. However, this is not the case, and the maximum intensity of the SHG is observed deeper in the sample where the IR is about 80 % of its peak intensity. This comes from the fact that the SHG beam is generated in the medium by the IR beam and then propagates with the same transport characteristics as the IG beam. This behavior will be well reproduced by the theoretical model. Regarding the wings deep inside the sample, an exponential decay is observed. This arises from the loss mechanisms at play in the sample coming from the finite transverse extent of the pellet, the linear absorption, and the nonlinearity in the case of the IR light.

III. THEORETICAL MODEL

A. Disorder model

In order to give more physical insights, we develop a theoretical model taking into account multiple scattering and second harmonic generation. To this end, we generalize the standard multiple scattering theory to derive transport equations for the fundamental and second harmonic intensities averaged over the configurations of the disorder. This derivation is similar to that of Ref. 32. However, there are three major differences: (1) we explicitly consider the time-dependent regime, (2) we drop the intensity decorrelation due to scatterer displacements and (3) we take absorption into account. This is important to note that this does not change the way the derivation is performed since these differences do not enter directly the computation of the second harmonic phase function. The interested reader can refer to App. B where the derivation is given. We present in the following only the assumptions and the main results. We first define a model of the disorder. The real material is made of a fine KDP powder containing particles of different shapes and sizes. Therefore, the most simple and natural disorder model consists of a continuous and complex fluctuating permittivity ǫ(r). The disorder microstructure is then characterized by a spatial correlation function chosen to be Gaussian, in the form

C ǫ (r -r ′ , ω) = δǫ(r, ω)δǫ * (r ′ , ω) = |∆ǫ(ω)| 2 C(r -r ′ ) (4) with C(r -r ′ ) = exp - |r -r ′ | 2 2ℓ 2 . ( 5
)
In the equations above, . . . represents an average over all disorder configurations (statistical average). δǫ(r, ω) = ǫ(r, ω) -ǫ(r, ω) is the fluctuating part of the permittivity, |∆ǫ(ω)| 2 is the amplitude of the correlation and ℓ is the correlation length. |∆ǫ(ω)| 2 depends on frequency since the permittivity ǫ is dispersive. However, ℓ involves only the geometrical structure of the disorder and thus does not depend on frequency. This also implies that the χ (2) nonlinearity is supposed to be correlated in a similar way (i.e., with the same correlation length). We thus have

C χ (r-r ′ , ω) = χ (2) (r, ω)χ (2) * (r ′ , ω) = |∆χ(ω)| 2 C(r-r ′ ). (6)

B. Transport equations

We first consider the case of the fundamental beam at frequency ω corresponding to λ = 1064 nm and also denoted by IR. The application of the standard multiple scattering theory leads to the no less standard Radiative Transfer Equation (RTE) given by [START_REF] Chandrasekhar | Radiative Transfer[END_REF] 1 v E (ω)

∂ ∂t + u • ∇ r + 1 ℓ e (ω) I(r, u, t, ω) = 1 ℓ s (ω) p(u • u ′ , ω)I(r, u ′ , t, ω)du ′ . (7)
Although the pellet is large enough, a diffusion-type equation is not adequate enough to describe the experimental results. Indeed, this kind of equation does not allow an accurate reconstruction of the fluxes at short times and/or at small depths. A transport type equation such as Eq. ( 7) is required. In this equation, I(r, u, t, ω) is the specific intensity, that can be seen as the radiative flux at position r, in direction u, at time t and at frequency ω. More precisely, the derivation from the first principles shows that it is given by the Wigner transform of the field. It reads

δ(k -k r )I(r, u, t, ω) = E r + s 2 , ω + Ω 2 × E * r - s 2 , ω - Ω 2 e -iku•s-iΩt ds dω 2π ( 8
)
where E is the electric field, E * being its complex conjugate counterpart. In this definition, the electric field is a scalar quantity. Indeed, we choose here to neglect the effect of the polarization which is a very good approximation in the multiple scattering regimes where the field can be considered to be fully depolarized [START_REF] Vynck | Polarization and spatial coherence of electromagnetic waves in uncorrelated disordered media[END_REF]. k r (ω) = n r (ω)k 0 is the real part of the effective wavevector where k 0 = ω/c is the wavenumber in vacuum. We recall that n r is the real part of the effective refractive index. It describes the phase velocity of the average field E (i.e., the field averaged over all possible configurations of the disorder). ℓ e (ω) and ℓ s (ω) are the extinction and scattering mean-free paths respectively. They are connected by the relation ℓ e (ω) -1 = ℓ s (ω) -1 + ℓ a (ω) -1 with ℓ a (ω) the absorption mean free path. p(u, u ′ , ω) is the phase function representing the fraction of power scattered from direction u ′ to direction u during a single scattering event. For the Gaussian disorder considered here, it is given by

p(u, u ′ , ω) ∝ p(|k r (ω)u -k r (ω)u ′ |)
where p(q) = exp -

q 2 ℓ 2 2 ( 9
)
and normalized such that p(u, u ′ , ω)du ′ = 1. As previously mentioned, v E (ω) is the energy velocity, well approximated by the phase velocity for non-resonant scatterers. Equation ( 7) can be easily interpreted using a random walk picture. In this picture, light undergoes a random walk whose average step is given by the scattering mean free path ℓ s (ω) and whose angular distribution at each scattering event is given by the phase function p(u, u ′ , ω). The absorption is taken into account along a path of length s through an exponential decay exp[-s/ℓ a (ω)].

Let us focus now on the second harmonic light at frequency 2ω corresponding to λ = 532 nm and also denoted by SHG. By generalizing the standard multiple scattering theory taking into account second harmonic generation, we obtain a RTE given by (see App. B for details)

1 v E (2ω) ∂ ∂t + u • ∇ r + 1 ℓ e (2ω) I(r, u, t, 2ω) = 1 ℓ s (2ω) p(u • u ′ , 2ω)I(r, u ′ , t, 2ω)du ′ + α p SHG (u, u ′ , u ′′ , ω)I(r, u ′ , t, ω) × I(r, u ′′ , t, ω)du ′ du ′′ . ( 10
)
This equation is very similar to Eq. ( 7) except for the presence of a non-linear source term. Because of the statistical average over disorder configurations, it can be shown that the dominant contribution in terms of paths for the second harmonic beam involves two SHG processes (one to generate the field at 2ω and the second to generate its complex conjugate) at the same location and time [START_REF] Samanta | Speckle decorrelation in fundamental and second-harmonic light scattered from nonlinear disorder[END_REF]. This implies that there is no phase shift to take into account between two fields, one at frequency ω and the second at frequency 2ω and there is no phase-matching condition as in standard SHG experiments without scattering. Thus the SHG process reduces to the product of two specific intensities at ω coming from two different directions u ′ and u ′′ . α is a factor containing all constants involved in the SHG process such as χ (2) . p SHG (u, u ′ , u ′′ , ω) is the SHG phase function describing the distribution of the SHG source term in direction u arising from two linear specific intensities coming from directions u ′ and u ′′ . In the case of the correlated disorder we consider here, it is given by

p SHG (u, u ′ , u ′′ , ω) ∝ p(|k r (2ω)u -k r (ω)u ′ -k r (ω)u ′′ |) (11
) where the p function is given in Eq. ( 9), and normalized such that p SHG (u, u ′ , u ′′ , ω)du ′ u ′′ = 1. Thus, the SHG phase function is directly related to the disorder correlation function in the same way as the standard phase function since the SHG and scattering processes both take place in the scatterers. We note that these equations are obtained in the perturbative approach meaning that the beam at ω is a source term for the beam at 2ω but the beam at 2ω has a negligible impact on the beam at ω.

C. Numerical simulations

Equations ( 7) and ( 10) are solved using a Monte Carlo scheme in geometries as close as possible to that in the real experiments (see Fig. 4 for details). In particular, the crystal grains have sizes ranging from 2 µm to 8 µm, which is large compared to the wavelength. The correlation length is taken such that k 0 ℓ = 3 for all simulations. Using this value and the experimental estimate of the real part n r of the effective refractive index, we get estimates for the anisotropy factors g(ω) = 0.95 and g(2ω) = 0.99. Then the values of the scattering mean-free paths are deduced from the experimental estimates of the transport mean-free paths using the relation ℓ s = ℓ t (1 -g). Experimental estimates are directly used for the absorption mean-free paths ℓ a . This shows that the disorder correlation length ℓ impacts only ℓ s . Recalling that ℓ t are small compared to the dimensions of the pellet, the diffusive regime applies. Thus ℓ t is much more important than ℓ s at large depth and time and the real value of k 0 ℓ weakly impacts the numerical results except at small depth and time where slight variations can be visible. In practice, three simulations are performed. The first is done to compute a 6D map (3 positional, 2 directional, and 1 temporal coordinates) of the specific intensity at ω. It gives also access to the linear detected flux at ω denoted by F IR (t). The second is done to evaluate the SHG detected flux at 2ω [i.e., F SHG (t)] using the previous map as a source term. And the last is performed to estimate the linear detected flux at 2ω [i.e., F IG (t)].

The results are shown in Figs. 5 and6 for the temporal and spatial profiles respectively. The specific parameters for each experiment are detailed in the captions and are chosen to mimic the experiments. Both figures show very good agreement with the experiment. Let us first consider the temporal results of Fig. 5. We clearly recover the correct exponential time decay of the flux at long times which is a feature of the diffusive regime. In particular, the long time-dependence of the flux at the fundamental frequency F IR (t) is given by exp(-Γ(ω)t).

In the diffusive regime, since all radiative quantities have the same time dependence, the specific intensity itself is also decaying according to the same law. This implies that the source term of the non-linear RTE has a long time dependence given by exp(-2Γ(ω)t). By convolving this source term with the time response of a short pulse, we simply deduce the long time behavior of the SHG flux at 2ω which is given by

F SHG (t) ∝ t 0 e -2Γ(ω)t ′ e -Γ(2ω)(t-t ′ ) dt ′ ∝    e -Γ(2ω)t -e -2Γ(ω)t 2Γ(ω) -Γ(2ω) if Γ(2ω) = 2Γ(ω), te -Γ(2ω)t otherwise. ( 12
)
For the powder considered here, 2Γ(ω) > Γ(2ω), which means that the long time behavior of the SHG flux is given by e -Γ(2ω)t . In other words, the IG and SHG beams have the same time decay at long times. This is clearly seen in Fig. 5. However, we note that a different behavior could have been observed depending on the order relation between 2Γ(ω) and Γ(2ω). We also note that Eq. ( 12) has been obtained under the assumption of exponential decay at all times for all quantities included in this calculation. Since this is valid only at long times, slight deviations could be observed numerically when comparing the IG and SHG beams. However, Fig. 5 shows that this is not the case. A quick comparison with Figure 2 reveals very good qualitative agreement with the experiments.

Regarding the spatial distribution of intensity displayed in Fig. 6, we clearly see that the green light in the linear domain shows a maximum right at the input

7LPHQV

)PD[) ,*0& ,50& 6+*0& í í FIG. 5. Numerically simulated temporal transmitted flux profiles of incident green (IG, blue solid line), fundamental infrared (IR, orange dotted line), and second-harmonic generated green (SHG, green dashed line) light for the L = 5 mm pellet. Inset: flux profiles on a semilog scale. We also observe a poorer SNR on the SHG curve. This is due to slow numerical convergence (only a fraction of the SHG light is collected on the detector, since the source term for the SHG beam is dispersed throughout the pellet and a large part of the SHG beam emerges from the other sides of the pellet). All parameters are chosen to mimic the experiment. The pellet is illuminated from the front face with an ultrashort Gaussian pulse of normalized waist k0wi = 2.95 × 10 4 (at frequency ω) and k0wi = 1.77 × 10 4 (at frequency 2ω) and under normal incidence. The detection is performed through the back face of the pellet [see Fig. 4 (a)]. The other parameters are given by k0a = 2.96 × 10 4 , k0L = 2.96 × 10 4 , k0ℓs(ω) = 1.28 × 10 3 , k0ℓs(2ω) = 1.89 × 10 2 , k0ℓa(ω) = 1 × 10 10 , k0ℓa(2ω) = 5.9 × 10 5 , nr(ω) = 1.43, nr(2ω) = 1.45 and k0ℓ = 3.

interface of the sample. The fundamental IR light, similarly, shows a maximum slightly deeper in the sample, obviously owing to the larger value of ℓ t , which is the consequence of the larger wavelength. The generated second-harmonic light peaks much deeper into the sample, reflecting the process of local generation and subsequent diffusion of the emitted light. These trends accurately reproduce the experimental observations in Fig. 3 in a quantitative manner, certifying the completeness of the numerical model. Overall, for the total spatiotemporal behavior, we can claim the agreement to be semi-quantitative.

IV. SUMMARY AND CONCLUSIONS

In conclusion, we have experimentally and numerically investigated the diffusion of light in a second-order nonlinear disordered material with a strong nonlinear coefficient. In the experiments, a cylindrical pellet of KDP microcrystals, packed in a random orientation, was employed. The spatial and temporal diffusion profiles of

zPP )PD[) ,*0& ,50& 6+*0&
FIG. 6. Numerically simulated spatial flux profiles of incident green (IG, blue solid line), fundamental infrared (IR, orange dotted line), and second-harmonic generated green (SHG, green dashed line) light for the L = 5 mm pellet. Numerical diffusion profiles show excellent agreement with the experimental results in Fig. 3. All parameters are chosen to mimic the experiment. The pellet is illuminated from the front face with a Gaussian beam of normalized waist k0wi = 2.95 × 10 4 (at frequency ω) and k0wi = 1.77 × 10 4 (at a frequency 2ω) and under normal incidence. The detection is performed along the side face of the pellet at different normalized depths k0d [see Fig. 4 (b)]. The other parameters are given by k0a = 2.96 × 10 4 , k0L = 2.96 × 10 4 , k0ℓs(ω) = 1.28 × 10 3 , k0ℓs(2ω) = 1.89 × 10 2 , k0ℓa(ω) = 1 × 10 10 , k0ℓa(2ω) = 5.9 × 10 5 , nr(ω) = 1.43, nr(2ω) = 1.45 and k0ℓ = 3. fundamental and second harmonic light have been measured, along with incident light at the second harmonic wavelength. The experimental data were used to estimate the transport parameters to the best possible accuracy. Next, the experimental results were supported by a theoretical model derived from first principles and led to the coupling of transport equations for the linear and non-linear beams. A Monte Carlo scheme has been used to solve this model numerically in the same geometry as in the experiment. Excellent agreement has been obtained in the spatial behavior seen in the experiments, along with very good qualitative agreement in the temporal behavior. Moreover, this model allowed a clear interpretation of the experimental results.

We note this work addresses spatial and temporal diffusion separately for both fundamental and second-harmonic light. In the domain of linear disorder, researchers have investigated spatio-temporal diffusion [START_REF] Pattelli | Spatio-temporal visualization of light transport in complex photonic structures[END_REF][START_REF] Sperling | Can 3D light localization be reached in 'white paint' ?[END_REF], which has very interesting implications in light transport. Naturally, spatio-temporal diffusion of second-harmonic light, simultaneously with the fundamental, will form an exciting frontier in nonlinear disorder physics. Indeed, the streak camera used for this work has the capability to simultaneously capture spatial and temporal diffusion. Besides, in our cylindrical geometry, there is a scope for exploration of transverse diffusion of second harmonic as well as fundamental light. These topics will form the object of our future interest, starting with the current numerical models. We expect the present work to lay down the platform for further modeling of light transport in nonlinear disorder.

where

R ⊥ (θ) = n 1 cos θ -n 2 1 -n1 n2 sin θ 2 n 1 cos θ + n 2 1 -n1 n2 sin θ 2 2 (A3)
and

R (θ) = n 1 1 -n1 n2 sin θ 2 -n 2 cos θ n 1 1 -n1 n2 sin θ 2 + n 2 cos θ 2 (A4)
are the Fresnel power reflection coefficients for incident light polarized perpendicular and parallel to the plane of incidence, respectively. Then the power internal reflection coefficient R averaged over the incident angle θ is given by [START_REF] Zhu | Internal reflection of diffusive light in random media[END_REF]

R = 3C 2 + 2C 1 3C 2 -2C 1 + 2 (A5)
where

C 1 = π 2 0 R(θ) sin θ cos θdθ (A6)
and

C 2 = π 2 0 R(θ) sin θ cos 2 θdθ. (A7)
Putting n 1 = n r and n 2 = 1, we get R(ω) = 0.52 (i.e. for λ = 1064 nm) and R(2ω) = 0.53 (i.e. for λ = 532 nm).

Appendix B: Derivations of the transport equations

This appendix is dedicated to the derivations of the transport equations at ω and 2ω. report here the main steps and we detail the computation of the SHG source term for the non-linear transport equation at 2ω. The derivations rely mainly on the standard multiple scattering theory, details of which can be found in various books and reviews [START_REF] Carminati | Principles of Scattering and Transport of Light[END_REF][START_REF] Rytov | Principles of Statistical Radiophysics[END_REF][START_REF] Apresyan | Radiation Transfer: Statistical and Wave Aspects[END_REF][START_REF] Van Rossum | Multiple scattering of classical waves: microscopy, mesoscopy and diffusion[END_REF]].

Scattering potential and Green function

One of the main building blocks of the multiple scattering theory is the scattering potential defined by

V (r, ω) = k 2 0 [ǫ(r, ω) -ǫ b (ω)] (B1)
where ǫ b (ω) is the permittivity of the background homogeneous medium (reference medium). In order to ensure that the perturbative method applied afterward is as accurate as possible, it is advisable to choose for ǫ b the statistical average of the permittivity (i.e., ǫ(r, ω) ) which leads to

V (r, ω) = k 2 0 δǫ(r, ω). (B2)
This potential describes the scattering process by the heterogeneities of the medium. The second ingredient is the background Green function that connects two consecutive scattering events. It is given by

G b (r -r ′ , ω) = exp(ik b |r -r ′ |) 4π|r -r ′ | (B3)
where k b = n b k 0 is the wavevector in the background medium. We note that this Green function is a scalar quantity here since we choose to neglect the polarization. Indeed, it can be shown that the polarization is completely washed out on average after a propagation distance of the order of the transport mean-free path inside a disordered medium [START_REF] Vynck | Polarization and spatial coherence of electromagnetic waves in uncorrelated disordered media[END_REF]. In particular, this expression given by Eq. (B3) corresponds to the Green function of the scalar homogeneous wave equation given by

∆G b (r -r ′ , ω) + k 2 b G b (r -r ′ , ω) = -δ(r -r ′ ). (B4)

Self-energy and intensity vertex

Two operators are used to describe the propagation of the field in disordered media: the self-energy Σ(r, r ′ , ω) entering the Dyson equation that governs the average field propagation and the intensity vertex Γ(r, r ′ , ρ, ρ ′ , ω) entering the Bethe-Salpeter equation that governs the field-field correlation evolution. Both contain an infinite series of scattering sequences that are statistically not factorizable. Considering a dilute medium (quantified further by the condition k r ℓ e ≫ 1), these scattering sequences can be seen as Taylor expansions. Thus, we can apply a perturbative approach and truncate the series to the first non-zero order. In the following, we consider only the frequency ω but the same result holds at 2ω (for the linear beams only).

For the self-energy Σ, this gives

Σ(r, r ′ ) = V (r, ω)G b (r -r ′ , ω)V (r ′ , ω) c (B5)
where • c represents a statistical average restricted to the connected part, i.e.

V G b V c = V G b V -V G b V .
Since the disorder correlation function C depends only on |rr ′ | (statistical homogeneity and isotropy), we have in the Fourier domain

Σ(k, k ′ , ω) = 8π 3 δ(k -k ′ ) Σ(k, ω). ( B6
)
This allows us to define the real part of the effective optical index by

n r (ω) = n b + Re Σ(k r , ω) 2n r k 2 0 (B7)
where we recall that k r (ω) = n r (ω)k 0 and k 0 = ω/c. This is a closed equation and the computation of n r in practice is not an easy task except in very dilute media. That being known, the extinction mean-free path is defined by

1 ℓ e (ω) = Im Σ(k r , ω) k r . ( B8
)
For the intensity vertex Γ, we have

Γ(r, r ′ , ρ, ρ ′ , ω) = V (r, ω)V * (ρ, ω) c δ(r ′ -r ′ )δ(ρ -ρ ′ ). ( B9
) By invoking one again the statistical homogeneity and isotropy, we get into the Fourier domain

Γ(k, k ′ , κ, κ ′ , ω) = 8π 3 δ(k -k ′ -κ + κ ′ )Γ(k, k ′ , κ, κ ′ , ω)
(B10) which leads to the following expressions of the phase function and of the scattering mean-free path:

1 ℓ s (ω) p(u, u ′ , ω) = 1 16π 2 Γ(k r u, k 0 u ′ , k 0 u, k 0 u ′ , ω) (B11) with the normalization p(u, u ′ , ω)du ′ = 1. (B12)
This finally leads to

1 ℓ s (ω) = k 2 0 |∆ǫ| 2 16π 2 4π C(q)du ′ (B13) p(u, u ′ , ω) = C(q) 4π C(q)du ′ -1 (B14) = k 2 r ℓ 2 exp[-q 2 ℓ 2 /2] 2π [1 -exp(-2k 2 r ℓ 2 )] (B15)
with q = k r (uu ′ ) and C(q) the Fourier transform of the disroder correlation function C. In the case of a non-absorbing medium where the imaginary part of the permittivity vanishes, we have from Eqs. (B8) and (B13)

1 ℓ e (ω) = 1 ℓ s (ω) = k 2 0 |∆ǫ| 2 16π 2 4π C(k r (u -u ′ ))du ′ (B16)
which implies that ℓ -1 a (ω) = 0. In practice for the numerical simulations performed in this study, the parameters n r , ℓ e and ℓ s are not computed using Σ and Γ but are given by experimental measurements. However, Eq. (B15) is used for the profile of the phase function.

Transport equation in the linear regime

Let us first consider the transport equation at ω in the linear regime. Considering a dilute medium such that k r ℓ e (ω) ≫ 1, we can show that the dominant term in the field-field correlation function is given by the well-known ladder diagram:

E(r, ω) E * (r, ω) E 0 E * 0 . ( B17
)
In this representation, the top line represents a path for the electric field E and the bottom line is for a path of its complex conjugate E * . Solid and dashed thick lines correspond to average Green functions (describing propagation between consecutive scattering events) and average fields respectively. Circles denote scattering events and vertical dashed lines represent statistical correlations between scattering events through Eq. (B9). This diagram shows that the field-field correlation function evolution essentially reduces to the propagation of an intensity. This implies that the field-field correlation function is governed by a transport equation, namely the Radiative Transfer Equation (RTE) that reads

1 v E (ω) ∂ ∂t + u • ∇ r + 1 ℓ e (ω)
I(r, u, t, ω)

= 1 ℓ s (ω)
p(u, u ′ , ω)I(r, u ′ , t, ω)du ′ (B18) which is Eq.( 7).

Transport equation in the second harmonic regime

Since we consider a perturbative approach for the SHG, the main point to address in that case is the expression of the source term. We can show theoretically and validate numerically [START_REF] Samanta | Speckle decorrelation in fundamental and second-harmonic light scattered from nonlinear disorder[END_REF] that the leading diagram is given by E(r, 2ω) E * (r, 2ω)

E 0 E 0 E * 0 E * 0 (B19)
where the squares denote the second harmonic processes. This diagram confirms that a transport equation for the SHG beam is still valid but with a source term given by the product of two specific intensities at ω. In terms of an equation, this source term is indeed given by S(r, ρ, 2ω) = G(rr ′ , 2ω) G * (ρρ ′ , 2ω)

× Γ SHG (r ′ , r ′′ , r ′′′ , ρ ′ , ρ ′′ , ρ ′′′ , 2ω) × E(r ′′ , ω)E * (ρ ′′ , 2ω) E(r ′′′ , ω)E * (ρ ′′′ , 2ω) × dr ′ dr ′′ dr ′′′ dρ ′ dρ ′′ dρ ′′′ (B20)

where Γ SHG is the SHG vertex given by Γ SHG (r, r ′ , r ′′ , ρ, ρ ′ , ρ ′′ , ω) = χ (2) (r, ω)χ (2) * (ρ, ω) c × δ(rr ′ )δ(rr ′′ )δ(ρρ ′ )δ(ρρ ′′ ). (B21)

In this last expression, we keep the complex conjugate notation for the second-order susceptibility χ (2) although it is a real quantity in order to remind that it corresponds to the complex conjugate field. In the Fourier domain and making use of the statistical homogeneity and isotropy again, we get Γ SHG (k, k ′ , k ′′ , κ, κ ′ , κ ′′ , ω) = 8π 3 δ(k-k ′ -k ′′ -κ+κ ′ +κ ′′ ) × Γ SHG (k, k ′ , k ′′ , κ, κ ′ , κ ′′ , ω). (B22)

From this, we can define a SHG phase function given by αp SHG (u, u ′ , u ′′ , ω) = 1 125π 5 Γ SHG (k r (2ω)u, k r (ω)u ′ , k r (ω)u ′′ , k r (2ω)u, k r (ω)u ′ , k r (ω)u ′′ , ω). (B23) α is a coefficient that takes into account all constants involved in the second harmonic generation and is such that the second harmonic phase function is normalized as p SHG (u, u ′ , u ′′ , ω)du ′ du ′′ = 1.

(B24)

Plugging the expression of the correlation function C χ leads to p SHG (u, u ′ , u ′′ , ω) ∝ exp -q 2 ℓ 2 2 (B25)

where q = k r (2ω)u -k r (ω)u ′ -k r (ω)u ′′ is the SHG scattering vector. We finally end up with the non-linear RTE given by

1 v E (2ω) ∂ ∂t + u • ∇ r + 1 ℓ e (

FIG. 1 .

 1 FIG. 1. (a) SEM image of particles. The scale bar is 10 µm. Inset: Distribution of particle size. (b) Schematic for temporal measurements. Notations: NC, Nonlinear Crystal; HBS, Harmonic Beam Splitter; M, Mirror; L, Lens; F1, Filter (for λ = 532 nm); OF, Optical Fiber; SC, Streak Camera. (c) Schematic of the experiment on spatial diffusion. Notations: L, Lens; F1; F2, Filters; SWIR, Short-wavelength Infra-Red; EMCCD, Electron Multiplying Charge Coupled Device. Inset images show the pellet with 1 cm diameter and length of 5 mm.

FIG. 2 .

 2 FIG. 2. Experimental temporal diffusion profiles of SHG(green squares) and IG (blue dots) in a pellet with L = 5 mm. The SHG signal peaks earlier in time. The orange dashed and red solid lines fit the exponential decays for SHG and IG respectively. Inset: diffusion profiles on the semilog scale show comparable decay lengths.

FIG. 3 .

 3 FIG.3. Experimental spatial intensity profiles of incident green (IG, blue solid line), fundamental infrared (IR, orange dotted line), and second-harmonic generated green (SHG, green dashed line) light for the L = 5 mm pellet. The SHG signal peaks much deeper in the sample, and IG peaks closest to the entry face inside the pellet, followed by IR.

FIG. 4 .

 4 FIG. 4. Geometry considered in the Monte Carlo simulations. The illumination is a Gaussian beam with waist wi impacting the front face of the pellet under normal incidence. The detection is performed through the back face of the pellet [(a), temporal simulation] or along the side face [(b), spatial simulation]. The cylinder contains the powder (pellet) and has a radius a and a height L. The refractive index mismatch between the cylinder (real part of the effective refractive index nr) and the external medium assumed to be air is also taken into account in the simulation.

  • u ′ , 2ω)I(r, u ′ , t, 2ω)du ′ +α p SHG (u, u ′ , u ′′ , ω)I(r, u ′ , t, ω)I(r, u ′′ , t, ω)du ′ du ′′

			2ω)	I(r, u, t, 2ω)
	=	1 ℓ s (2ω)	p(u (B26)
	which is Eq. (10).

ACKNOWLEDGMENTS R.S. and S.M. acknowledge the support from Kartik K. Iyer and Sanjay K. Upadhyay during the samplemaking. We also thank Bhagyashri Chalke for her help in the SEM measurements. R.S. appreciates the support from Ashwin K. Boddeti for the initial experiments. R.S. thanks Sandip Mondal for the fruitful discussions. We express our gratitude towards the Department of Atomic Energy, Government of India for funding for Project Identification No. RTI4002 under the DAE OM No 1303/1/2020/R&D-II/DAE/5567 dated 20.8.2020; Swarnajayanti Fellowship, Department of Science and Technology, Ministry of Science and Technology, India. This work has also received support under the program Investissements d'Avenir launched by the French Government. The authors declare no conflicts of interest.

In the pellet, the volume fraction of the KDP is 0.88. The Maxwell Garnett equation reads [START_REF] Choy | Effective Medium Theory[END_REF]

where ε r is the real part of the effective dielectric constant of the medium, ε i of the inclusions, and ε m of the matrix; δ i is the volume fraction of the inclusions. In our case, the air percentage is smaller than the KDP powder. We can therefore assume that air is the inclusion and KDP is the matrix. Putting δ i = 0.12, ε i = 1 and ε m = n 2 m = (1.5124) 2 we get n r (2ω) = 1.45 at λ = 532 nm. At λ = 1064 nm, the refractive index of KDP is 1.4938 [START_REF]KDP refractive index dataset[END_REF]. The calculated value of n r (ω) is 1.4336.

Determination of the internal reflection coefficient R

We can obtain an estimate of the reflection coefficient using Fresnel's law. We assume that the direction and polarization of diffusing light incident on the boundary from inside the sample are completely random and the sample surface is flat. For an angle of incidence θ, the power reflection coefficient R(θ) averaged over polarization is