V Dorier

J Lampart

S Guerin

H R Jauslin

S Guérin

Canonical quantization for quantum plasmonics with finite nanostructures

Keywords: numbers: 42.50.Nn, 71.45.Gm, 42.50.Ct Phys.Rev

published or not. The documents may come

Canonical quantization for quantum plasmonics with finite nanostructures

I. INTRODUCTION

Quantum effects in plasmonic systems are expected to open a whole branch of new perspectives both in fundamental physics and in technological applications [START_REF] Tame | Quantum plasmonics[END_REF][START_REF] Bozhevolnyi | The case for quantum plasmonics[END_REF][START_REF] Barnes | Surface plasmon subwavelength optics[END_REF][START_REF] Altewischer | Plasmon-assisted transmission of entangled photons[END_REF][START_REF] Martino | Observation of Quantum Interference in the Plasmonic Hong-Ou-Mandel Effect[END_REF][START_REF] Dheur | Marquier Remote preparation of single-plasmon states[END_REF][START_REF] Novotny | Principles of nano-optics[END_REF], like e.g. in quantum information, sensors, and generally in new types of circuits involving combinations of electronic and photonic effects. In the theories for quantum plasmonics there are two main approaches that have been widely used: (a) a phenomenological approach formulated in terms of quantum Langevin equations [START_REF] Gruner | [END_REF][START_REF] Vogel | Quantum Optics[END_REF][START_REF] Scheel | Macroscopic quantum electrodynamics -concepts and applications[END_REF][START_REF] Buhmann | Dispersion forces I[END_REF] and (b) microscopic oscillator models for the medium coupled to the electromagnetic field [START_REF] Huttner | Quantization of the electromagnetic field in dielectrics[END_REF][START_REF] Suttorp | Field quantization in inhomogeneous absorptive dielectrics[END_REF][START_REF] Suttorp | Fano diagonalization of a polariton model for an inhomogeneous absorptive dielectric[END_REF][START_REF] Philbin | Canonical quantization of macroscopic electromagnetism[END_REF]. The latter are variations of models of the type first proposed by Hopfield [START_REF] Hopfield | Theory of the Contribution of Excitons to the Complex Dielectric Constant of Crystals[END_REF]. The quantization of a model for a homogeneous bulk medium was first treated by Huttner and Barnett [START_REF] Huttner | Quantization of the electromagnetic field in dielectrics[END_REF]. Extensions to inhomogeneous media were treated in [START_REF] Suttorp | Fano diagonalization of a polariton model for an inhomogeneous absorptive dielectric[END_REF][START_REF] Philbin | Canonical quantization of macroscopic electromagnetism[END_REF]. The main criterion for the choice of the microscopic model is that, if one integrates the equations for the medium and one inserts the obtained currents into the microscopic Maxwell equations, one should obtain the macroscopic Maxwell equations. The elementary quanta of the model are called plasmon-polaritons.

Both approaches led to the following three basic relations:

[A] The diagonal plasmon-polariton Hamiltonian is of the form

Ĥ = d 3 r ∞ 0 dν ν Ĉ † ( r, ν) • Ĉ( r, ν), (1)
where Ĉ † ( r, ν) and Ĉ( r, ν) are bosonic creationannihilation operators.

[B] The electric field is expressed in terms of the bosonic operators as

Ê( r) = µ 0 πc 2 ∞ 0 dν d 3 r ν 2 1/2 i ( r ,ν)G( r, r ,ν) Ĉ( r ,ν)+h.c., (2)
* jauslin@u-bourgogne.fr where i ( r , ν) is the imaginary part of the dielectric coefficient, and G( r, r , ν) is a Green tensor satisfying

∇ ∧ ∇ ∧ -( r, ν) ν 2 c 2 Ḡ( r, r , ν) = 1 3 δ( r -r ). (3)
[C] The Green tensor satisfies the following identity

ν 2 c 2 d 3 r i ( r) G T ( r, r A )G * ( r, r B ) = Im G( r A , r B ). ( 4
)
These three relations may be justified in the case of a bulk medium, i.e. when the dispersive and dissipative medium extends over all space R 3 [8-11, 14, 15]. In this article we show that within the microscopic oscillator models neither of these three relations is valid for a finite medium, e.g. in particular for a system with finite metallic or dielectric nano-structures. We present a procedure to construct the complete diagonalization of the Hamiltonian, Eq. ( 22), and the electric field operator, Eq. ( 25), adapted to finite media, which replace Eqs. ( 1) and ( 2). Although we do not use the Green tensor (3) explicitly in our construction, we remark that for finite media a boundary term has to be added in Eq. (4), as described in [START_REF] Drezet | Equivalence between the Hamiltonian and Langevin noise description of plasmon-polaritons in a dispersive and lossy inhomogeneous medium[END_REF]Appendix A].

It has been stated in several instances in the literature [START_REF] Di Stefano | Mode expansion and photon operators in dispersive and absorbing dielectrics[END_REF][START_REF] Drezet | Description of spontaneous photon emission and local density of states in the presence of a lossy polaritonic inhomogeneous medium[END_REF] that the formula (2) for the electric field observable cannot be generally correct, and that it has to be completed by including a term reflecting the contribution of the free electromagnetic field, while Eq. ( 2) describes only the electric field produced by the charges of the medium. A simple argument showing that (2) cannot be the complete formula for a finite medium is that if we take the limit where the size of the medium goes to zero, according to (2) the electric field observable would disappear. There must be another contribution that converges to the free electric field in the limit of vanishing medium.

Some extensions of Eq. (2) were proposed in Refs. [START_REF] Di Stefano | Mode expansion and photon operators in dispersive and absorbing dielectrics[END_REF][START_REF] Drezet | Description of spontaneous photon emission and local density of states in the presence of a lossy polaritonic inhomogeneous medium[END_REF], consisting in adding the free electric field E free :

Êtot = Êfree + Êmedium , (5)
with Êmedium given by (2). This approach is not quite satisfactory, since Êfree was added without a justification arXiv:1810.08014v2 [quant-ph] 2 Oct 2019 from the diagonalization of the microscopic model. This might explain why this procedure was not widely adopted in the literature [START_REF] Khanbekyan | Input-output relations at dispersing and absorbing planar multilayers for the quantized electromagnetic field containing evanescent components[END_REF].

In this article, we show that the electric field operator has indeed two components: one that converges to the free electric field when the size of the medium goes to zero, and another one that vanishes in that limit. However, our result shows that the first component is not just the free electric field, but a modification induced by the presence of the medium. The second component is not equal to Eq. ( 2), but to an expression that is modified due to the presence of the first term. We obtain this result by diagonalizing the microscopic model for a finite medium, as opposed to an infinite bulk medium as in Refs. [START_REF] Huttner | Quantization of the electromagnetic field in dielectrics[END_REF][START_REF] Suttorp | Field quantization in inhomogeneous absorptive dielectrics[END_REF][START_REF] Suttorp | Fano diagonalization of a polariton model for an inhomogeneous absorptive dielectric[END_REF][START_REF] Philbin | Canonical quantization of macroscopic electromagnetism[END_REF]. Furthermore, we show that the diagonalized Hamiltonian is not of the form [START_REF] Tame | Quantum plasmonics[END_REF]. It can be shown that, in the limit of vanishing coupling between the medium and the electromagnetic field, the operator (1) tends toward the Hamiltonian of the uncoupled medium, and the Hamiltonian of the uncoupled electromagnetic field would be missing. Finally, we compare in perturbation theory the electric field observable derived from our model with the one proposed in the literature.

II. THE MICROSCOPIC PLASMON-POLARITON MODEL

We consider the microscopic model for the classical electromagnetic field in interaction with a linear medium used in Ref. [START_REF] Philbin | Canonical quantization of macroscopic electromagnetism[END_REF]. We consider a non-magnetic medium for simplicity, but our construction can easily be extended to magnetic media and to the models of similar type considered e.g. in Refs. [START_REF] Suttorp | Field quantization in inhomogeneous absorptive dielectrics[END_REF][START_REF] Suttorp | Fano diagonalization of a polariton model for an inhomogeneous absorptive dielectric[END_REF]. It can be written in Hamiltonian form by choosing the following pairs of canonically conjugate variables: Π X ( r, ν), X( r, ν) for the harmonic oscillators describing the medium, the vector potential A( r), and its canonically conjugate variable

Π A ( r) := -D( r) = -0 E( r) - ∞ 0 dν α( r, ν) X( r, ν). ( 6
)
The Hamiltonian is

H = H em + H med + H int (7)
with

H em = 1 2 d 3 r 1 0 Π 2 A ( r) - 1 µ 0 A( r) • ∆ A( r) , (8a)
H med = 1 2 Vm d 3 r ∞ 0 dν Π 2 X ( r, ν) + ν 2 X 2 ( r, ν) , (8b)
H int = 1 0 Vm d 3 r Π A ( r) • ∞ 0 dν α( r, ν) X( r, ν) + 1 2 0 Vm d 3 r ∞ 0 dν α( r, ν) X( r, ν) 2 , (8c)
where V m is the volume occupied by the medium, and α( r, ν) is the coupling. It is related to the imaginary part i of the dielectric coefficient of the medium by α 2 ( r, ν) = 2 0 ν i ( r, ν)/π. This microscopic model produces the macroscopic Maxwell equations when the medium variables are integrated [START_REF] Philbin | Canonical quantization of macroscopic electromagnetism[END_REF]. We will use the following real basis of transverse generalized eigenfunctions of the Laplacian,

ϕ k,σ,π ( r) =    1 2π 3/2 ε σ ( k) cos( k • r), π = c, 1 2π 3/2 ε σ ( k) sin( k • r), π = s, (9)
where σ = ± is the index for the basis of linear polarization, and ε ± ( k) are two real unit vectors orthogonal to k and to each other.

In a first step, we apply the simple canonical transformation ( Π A , A) → (p, q)

Π A ( r) = dκ √ 0 ω κ ϕ κ ( r)q κ , (10a)
A( r) = -dκ 1 √ 0 ω κ ϕ κ ( r)p κ . (10b)
Here, we use the abbreviation for the labels κ := ( k, σ, π), and we denote by ω κ = c| k| the eigenfrequencies. The integral over dκ is an abridged notation for a combination of integrals and sums. In these variables, the coupled Hamiltonian takes the following form

H = 1 2 P • P + 1 2 Q • Ω 2 Q, (11)
where

P := p Π X , Q := q X , (12)
Ω 2 = Ω 2 0 + V. (13)
Ω 0 and Ω are called the frequency operators of the uncoupled and coupled systems, and V is the coupling operator which can be explicitly extracted from Eq. (8c). Since V acts as a multiplication by α( r, ν) and an integration over the volume of the medium, it is uniquely defined for every dielectric coefficient(hence for every material) and for every geometry of the medium.

III. GENERAL QUANTIZATION PROCEDURE

The quantization of a Hamiltonian of the general form [START_REF] Buhmann | Dispersion forces I[END_REF] and its expression in terms of creation-annihilation operators can be formulated in a systematic way [START_REF] De Bièvre | Where's that quantum?[END_REF][START_REF] De Bièvre | Local states of free bose fields; proceedings of "Summer School on large Coulomb systems[END_REF] in terms of an orthonormal basis of generalized eigenfunctions of

Ω 2 ψ λ,d λ = λ 2 ψ λ,d λ , (14)
where d λ is a label for the degeneracy indices. This quantization procedure can be thought off as a generalization of the usual quantization of the free electromagnetic field, in which the plane waves, which are eigenfunctions of the Laplacian, are replaced by the eigenfunctions of Ω 2 . The uncoupled Ω 2 0 has continuous spectrum. Under some conditions on the coupling V this is also the case for the coupled operator Ω 2 . The diagonal quantized Hamiltonian can then be written as

Ĥ = ∞ 0 dλ d λ λ B † ψ λ,d λ Bψ λ,d λ . ( 15
)
The notation d λ for the degeneracy is an abridged notation for a combination of sums and integrals. The creation-annihilation operators act in the bosonic Fock space F B (P) [START_REF] De Bièvre | Where's that quantum?[END_REF][START_REF] De Bièvre | Local states of free bose fields; proceedings of "Summer School on large Coulomb systems[END_REF][START_REF] Reed | Methods of Modern Mathematical Physics[END_REF], with

B † ψ |∅ = |ψ , ( 16
)
where |∅ is the vacuum state and ψ are functions corresponding to the coordinates of the classical phase space P that enters in the definition of the Fock space. Their interpretation is that B † ψ creates one quantum on the phase space mode ψ. The orthonormality relation of the generalized eigenfunctions of Ω 2 implies the canonical commutation relations. The canonical observables can be expressed in terms of the creation-annihilation operators as

Q = 2 ∞ 0 dλ d λ λ -1/2 ψ λ,d λ Bψ λ,d λ + h.c. , (17a)
P = -i 2 ∞ 0 dλ d λ λ 1/2 ψ λ,d λ Bψ λ,d λ -h.c. . ( 17b
)

IV. QUANTIZATION OF THE PLASMON-POLARITON MODEL

We can now apply this general procedure to the model ( 7)-( 8). As we show below, the (squared) frequency operator Ω 2 has, in general, a continuous spectrum with two distinct families of generalized eigenfunctions, ψ e and ψ m ,

Ω 2 ψ e ω,d ω = ω 2 ψ e ω,d ω , (18a)
Ω 2 ψ m ν,d ν = ν 2 ψ m ν,d ν , (18b)
where the degeneracy indices are d ω := (ϑ, φ, σ, π) with ϑ and η the angles in spherical coordinates of k, and d ν := ( r, j) with j = 1, 2, 3 labeling the three components of the fields.

The eigenfunctions have the same structure as Eq. ( 12), i.e. they are composed of two blocks:

ψ e ω,d ω = u e ω,d ω v e ω,d ω , ψ m ν,d ν = u m ν,d ν v m ν,d ν . ( 19
)
The functions u in the upper block depend on the variables ( k, σ, π) associated with the electromagnetic field, and the functions v in the lower block depend on the (ν, r, j) of the oscillators. In the limit where the coupling V vanishes, the ψ e tend to eigenfunctions of the free electromagnetic field and the ψ m tend to eigenfunctions of the uncoupled medium. They satisfy the Lippmann-Schwinger equations [23, Sect XI.6, p.98]:

ψ e ω,d ω = φ e ω,d ω -Ω 2 0 -ω 2 ∓ i0 + -1 V ψ e ω,d ω , (20a) ψ m ν,d ν = φ m ν,d ν -Ω 2 0 -ν 2 ∓ i0 + -1 V ψ m ν,d ν , (20b)
with ∓i0 + two possible ways to avoid the singularity in the complex plane. The functions φ are eigenfunctions of the uncoupled operator Ω 2 0 :

φ e ω,d ω (ω , d ω ) = δ(ω -ω )δ d ω ,d ω 0 , (21a)
corresponding to a plane wave in Fourier representation, and

φ m ν, r,j (ν , r , j ) = 0 δ(ν -ν )δ( r -r )δ j,j , (21b)
corresponding to an oscillator at position r oriented in direction e j . This gives a one-to-one correspondence between the coupled and the uncoupled generalized eigenfunctions. We will show below that the correspondence is unitary, which implies that the spectrum of Ω 2 has the same degeneracy structure as the one of the uncoupled operator Ω 2 0 . The diagonalized Hamiltonian takes the form

Ĥ = ∞ 0 dω d ω ω D † ω,d ω Dω,d ω + ∞ 0 dν Vm d 3 r ν Ĉ † ν, r • Ĉν, r , (22)
where we have introduced the notation Ĉν, r,j := Bψ m ν,d ν ,

Dω,d ω := Bψ e ω,d ω .

This expression resembles the Hamiltonian (1) to which one would have added a free-field contribution. However, the creation-annihilation operators that appear in Eq. ( 22) are a deformation of the ones of the uncoupled model.

In the exterior of the medium, the last term in Eq. ( 6) vanishes, and Eq. (10a) gives

Ê( r) = - 1 √ 0 dκ ω κ ϕ κ ( r)q κ . (24)
Using Eq. ( 12), the electric field observable can be written in terms of creation-annihilation operators on both sets of eigenmodes:

Ê( r) = Êe ( r) + Êm ( r), (25) Êe
( r) = dκ β κ ( r) ∞ 0 dω √ ω d ω u e ω,d ω(κ) Dω,d ω + h.c. , Êm ( r) = dκ β κ ( r) ∞ 0 dν √ ν d ν u m ν,d ν(κ) Ĉν,d ν + h.c. ,
with β κ ( r) := -2 0 ω κ ϕ κ ( r), and the sums over degeneracy indices correspond to the following combinations:

d ω := ω 2 c 3 π/2 0 dϑ sin ϑ 2π 0 dη σ=± π=c,s , (26)
d ν := Vm d 3 r 3 j=1 . ( 27
)
Equation ( 25) is the main result of this construction. In particular, it exhibits the appropriate convergence in the limit of zero coupling. In summary, the electric field observable at a point r in the exterior of the medium can be written according to Eq. ( 25) as the sum of two terms. When the coupling goes to zero, one can show that the first term Êe goes to the free electric field:

Êe ( r) → -dω d ω ω 2 0 ϕ ω,d ω ( r) D0 ω,d ω + D0 † ω,d ω , ( 28
)
and the second term Êm vanishes. However, for an arbitrary coupling, the first term is not equal to the free electric field, but to a deformation due to the presence of the medium.

V. STRUCTURE OF THE EIGENFUNCTIONS

The uncoupled operator Ω 2 0 has absolutely continuous spectrum of infinite degeneracy, labeled by the indices d ω for the electromagnetic field and d ν for the oscillators of the medium. We will now explain why the spectrum of the coupled operator Ω 2 = Ω 2 0 +V has the same structure, including degeneracy, for the type of couplings that occur in plasmon-polariton models for a finite-size medium.

The methods and the intuition in the theory of continuous spectra are closely related to scattering theory. A central object in scattering theory is the Møller wave operator [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF]Sect XI.3 p.18]

W ± (Ω 2 , Ω 2 0 ) = lim t→∓∞ e iΩ 2 t e -iΩ 2 0 t . ( 29
)
If this operator is well defined (i.e. if the strong limit t → ∓∞ exists) and is unitary we have

Ω 2 = W ± Ω 2 0 W † ± , (30)
so Ω 2 is unitarily equivalent to Ω 2 0 and its spectrum must also be continuous with the same degeneracy. The generalized eigenfunctions are then given by ψ λ,d λ = W ± φ λ,d λ . The unitarity of W ± ensures that the eigenfunctions satisfy the usual relations of orthogonality and completeness.

The coupling V in the plasmon-polariton models belongs to a class that has been thoroughly studied since pioneering works by Friedrichs [START_REF] Friedrichs | Über die Spektralzerlegung eines Integraloperators[END_REF][START_REF] Friedrichs | On the perturbation of continuous spectra[END_REF][START_REF] Faddeev | On the Friedrichs Model in the Theory of Perturbations of a Continuous Spectrum[END_REF]. The wave operator W ± is an integral operator that can be constructed using the methods of [START_REF] Yafaev | Mathematical Scattering Theory[END_REF] (note [START_REF]Under appropriate continuity assumptions on α(ν) and with an ultraviolet cut-off of arbitrary size on the electromagnetic field, one can apply [27[END_REF]). The choice of sign of W ± induces two possible sets of eigenfunctions verifying the Lippmann-Schwinger equation ( 20) either with -i0 + (for W + ) or with +i0 + (for W -). The two sets of eigenfunctions are related by the unitary map S = W † + W -, which is the scattering operator [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF]. Further analysis reveals that Ω 2 has no eigenvalues and purely absolutely continuous spectrum, provided the coupling function α(ν) does not vanish for any ν > 0. This assumption is satisfied for the usual models for dissipative and dispersive media. Consequently, W ± is unitary and this shows that the generalized eigenfunctions ψ λ,d λ have the structure we claim. Figure 1 is a sketch of this preservation of the spectral structure.

VI. COMPARISON IN PERTURBATION THEORY

For a precise comparison between our construction and the results of Refs [START_REF] Gruner | [END_REF][START_REF] Vogel | Quantum Optics[END_REF][START_REF] Scheel | Macroscopic quantum electrodynamics -concepts and applications[END_REF][START_REF] Buhmann | Dispersion forces I[END_REF][START_REF] Huttner | Quantization of the electromagnetic field in dielectrics[END_REF][START_REF] Suttorp | Field quantization in inhomogeneous absorptive dielectrics[END_REF][START_REF] Suttorp | Fano diagonalization of a polariton model for an inhomogeneous absorptive dielectric[END_REF][START_REF] Philbin | Canonical quantization of macroscopic electromagnetism[END_REF][START_REF] Hopfield | Theory of the Contribution of Excitons to the Complex Dielectric Constant of Crystals[END_REF][START_REF] Di Stefano | Mode expansion and photon operators in dispersive and absorbing dielectrics[END_REF][START_REF] Drezet | Description of spontaneous photon emission and local density of states in the presence of a lossy polaritonic inhomogeneous medium[END_REF][START_REF] Khanbekyan | Input-output relations at dispersing and absorbing planar multilayers for the quantized electromagnetic field containing evanescent components[END_REF] one must be able to solve the Lippmann-Schwinger equations [START_REF] De Bièvre | Where's that quantum?[END_REF] in a given situation. A complete numerical solution which includes the high degree of degeneracy of the model will be the subject of future work. We present here a preliminary result where the equations [START_REF] De Bièvre | Where's that quantum?[END_REF] are solved in first order perturbation theory, i.e., where the volume of the medium and the dissipation are small. The details of the calculation are described in the appendix. Once the obtained functions u e and u m are inserted in the electric field operator [START_REF] Friedrichs | On the perturbation of continuous spectra[END_REF], it reads

Ê( r) = Ê0 ( r) + µ 0 πc 2 ∞ 0 dν Vm d 3 r ν 2 1/2 i (ν, r ) G ± 0 (ν; r , r) + F nf (ν; r , r) Ĉν, r + h.c., (31)
with r a position in the exterior of the medium. Ê0 is of the same form as the electric field operator in vacuum [START_REF]Under appropriate continuity assumptions on α(ν) and with an ultraviolet cut-off of arbitrary size on the electromagnetic field, one can apply [27[END_REF]. It is however not strictly the same, since the annihilation operator D acts on the plasmonic Fock space, which has a modified ground state. G ± 0 are two possible Green functions in vacuum:

G ± 0 ( r , r) = c 2 d k σ,π ϕ k,σ,π ( r ) ⊗ ϕ k,σ,π ( r) ω 2 -ν 2 ∓ i0 + , (32)
and F nf is a near-field term, which tends quickly to zero for r far from the surface of the medium:

F ij nf ( r , r) = -c 2 4πν 2 | r -r | 3 δ ij - 3(r i -r i )(r j -r j ) | r -r | 2 . ( 33
)
The two possible Green tensors G ± 0 can be obtained for a given choice of the Møller operator W ± . One can show that G - 0 satisfies the Sommerfeld radiation condition at infinity, with only outgoing waves to infinity, whereas G + 0 contains only incoming waves from infinity. It is customary to use the Sommerfeld radiation condition. Both representations can be used, however, since they are linked by a unitary transformation.

The expression ( 31) is obtained at first perturbative order in √ i . The link with Eq. ( 2) can be done by also developping the Green equation (3) to first order in √ i , which corresponds to replacing the Green tensor in (2) by its vacuum expression G 0 . We conclude that to first order in √ i , the expression (31) takes the same form as the formulas proposed in Refs. [START_REF] Di Stefano | Mode expansion and photon operators in dispersive and absorbing dielectrics[END_REF][START_REF] Drezet | Description of spontaneous photon emission and local density of states in the presence of a lossy polaritonic inhomogeneous medium[END_REF] (i.e., an electric field of the form (2) but with an extra term associated with the free field). There is also an additional near-field term in (31); this term diverges at the surface, and it is thus beyond the range of validity of perturbation theory.

It is worth mentioning that the "free-field term" Ê0 is modified at higher orders of the perturbative development, which contradicts the idea that adding the contribution of the free field only is sufficient to achieve a complete quantization of the model.

VII. CONCLUSION AND OUTLOOK.

In summary, unlike in Eq. ( 1), the diagonalized coupled Hamiltonian ( 22) is of the same form as the uncoupled one: it has absolutely continuous and positive spectrum with degeneracy labeled by the indices d ω and d ν .

Since the starting point for the diagonalization procedure in Refs. [START_REF] Suttorp | Fano diagonalization of a polariton model for an inhomogeneous absorptive dielectric[END_REF][START_REF] Philbin | Canonical quantization of macroscopic electromagnetism[END_REF] is the assumption of a diagonal Hamiltonian of the form (1), the preceding arguments imply that the results of Refs. [START_REF] Suttorp | Fano diagonalization of a polariton model for an inhomogeneous absorptive dielectric[END_REF][START_REF] Philbin | Canonical quantization of macroscopic electromagnetism[END_REF] cannot be applied to finite size media, in particular to nano-structures, where the Hamiltonian has the form [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF]. The Fanotype diagonalization technique [START_REF] Fano | Effects of Configuration Interaction on Intensities and Phase Shifts[END_REF], which was adapted in Refs. [START_REF] Suttorp | Fano diagonalization of a polariton model for an inhomogeneous absorptive dielectric[END_REF][START_REF] Philbin | Canonical quantization of macroscopic electromagnetism[END_REF] was conceived initially to treat a discrete mode coupled to a continuum. A characteristic phenomenon is that, under some hypothesis, the discrete mode "dissolves" in the continuum, as was first proven by Friedrichs [START_REF] Friedrichs | On the perturbation of continuous spectra[END_REF] and rederived by Fano [START_REF] Fano | Effects of Configuration Interaction on Intensities and Phase Shifts[END_REF], and the coupled system has only continuous spectrum. This is an intuition that may have led to postulate a diagonal Hamiltonian of the form (1). However, in the present models for finite media we have two coupled continua instead of a discrete mode and a continuum, and the Fano-Friedrichs mechanism does not apply. Instead we need an analysis of what happens when two continua interact, which is what we presented in this article.

We remark that the fact that the frequency operator Ω is diagonalized by a unitary transformation, given by the wave operator W + , guarantees that the operators Ĉν, r , Dω,d ω automatically satisfy the bosonic commutation relations. In the method of diagonalization used in Refs. [START_REF] Suttorp | Fano diagonalization of a polariton model for an inhomogeneous absorptive dielectric[END_REF][START_REF] Philbin | Canonical quantization of macroscopic electromagnetism[END_REF], the commutation relations have to be imposed a posteriori, which is a technically non-trivial operation, in particular if one wants to extend it to include the two families of boson operators as in Eq. [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF].

The relations (1)-(4) had been constructed implicitly for bulk systems, but they have been widely used to analyze phenomena in systems with finite size media, like e.g. enhanced spontaneous emission, Purcell, Casimir and Polder effects, superradiant emission, among many others [START_REF] Buhmann | Dispersion forces I[END_REF]31]. These applications have to to be reassessed taking into account the modified perspective including the missing terms, e.g. in the electric field. A natural question is whether the bulk formulas (1)-( 4) can be a good approximation to the complete expressions, in some particular regimes and beyond the perturbative approach described in Section VI.

For applications it will be necessary to develop or adapt efficient numerical methods to solve the Lippmann-Schwinger integral equations for media of different geometries (see e.g. Ref [START_REF] Atkinson | Numerical solution of Fredholm integral equations of the second kind[END_REF]), e.g. to treat the fields of nanoantennas [START_REF] Sauvan | Theory of the Spontaneous Optical Emission of Nanosize Photonic and Plasmon Resonators[END_REF][START_REF] Franke | Quantization of quasinormal modes for open cavities and plasmonic cavity quantum electrodynamics[END_REF]. An analytical and numerical study of these equations and the consequences for spontaneous emission will appear elsewhere. The present results apply to a large class of open quantum systems coupled to a bath, featuring dispersion and dissipation.

FIG. 1 .

 1 FIG.1. Spectral structure of the plasmon-polariton model before and after diagonalization. The two continuous spectra structure is preserved by the unitarity of the diagonalization transformation.

Acknowledgments

This work was supported by the "Investissements dAvenir" program, project ISITE-BFC / IQUINS (ANR-15-IDEX-03), QUACO-PRC (ANR-17-CE40-0007-01) and the EUR-EIPHI Graduate School (17-EURE-0002). We. also acknowledge support from the European Unions Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 765075 (LIMQUET)

Appendix A: First order perturbation theory

We show here that a perturbation development to first order of the Lippmann-Schwinger equation results in the expression (31) of the electric field operator. To lighten the notation, we show it only for G - 0 which satisfies the Sommerfeld radiation condition.

The Lippmann-Schwinger equations ( 20) take the general form

If the coupling V is small, Eq. (A1) can be developed as

In first order, we neglect all contributions of V n with n > 1, we thus have

How V acts on a general vector of the Hilbert space is deduced from the Hamiltonian of interaction (8c). It takes the following form:

where we recall the concise notation κ = ( k, σ, π). Note that we absorbed 0 by making the change of variable α → √ 0 α. We then apply V to the uncoupled eigenfunctions φ e/m given by Eq. ( 21) and we insert it into Eq. (A3) to obtain the solutions of the Lippmann-Schwinger equations in a small coupling regime. They take the form [START_REF] Khanbekyan | Input-output relations at dispersing and absorbing planar multilayers for the quantized electromagnetic field containing evanescent components[END_REF] with:

Now that the eigenfunctions of Ω 2 have been found, we can use them to calculate the electric field in the exterior of the medium. We insert u e and u m into Eqs. ( 25) and we obtain after some manipulations:

with

We express the integrals in spherical coordinates with Eq. ( 26) and we use the identity

for ω, ν > 0, which gives

The first term in (A10) can be rewritten as

Because of Eq. (A11) and the completeness of the eigenfunctions ϕ ω,d ω , we have dω g(ω; r , r) = δ T ( r -r ), with the transverse delta function defined as [START_REF] Stewart | Longitudinal and transverse components of a vector field[END_REF]

Since r is defined only inside the medium whereas r is taken in the exterior, the first term in the RHS vanishes.

Hence we obtain:

where

In (A13) we recognize the Green function in vacuum:

G - 0 (ν; r , r) = c 2 dω g(ω; r , r) ω 2 -ν 2 + i0 + = c 2 d 3 k σ,π ϕ k,σ,π ( r ) ⊗ ϕ k,σ,π ( r)

Introducing it back into Eq. (A7b) and using 0 µ 0 c 2 = 1, we obtain the electric field (31).


Shuoqi Li

Mingchao Liu

Dorian Hanaor

Yixiang Gan

Dynamics of Viscous Entrapped Saturated Zones in Partially Wetted Porous Media

Keywords: Porous media, multiphase fluid flow, drainage process, saturated cluster, Bond number

As a typical multiphase fluid flow process, drainage in porous media is of fundamental interest in nature and industrial applications. During drainage processes in unsaturated soils and porous media in general, saturated clusters, in which a liquid phase fully occupies the pore space between solid grains, affect the relative permeability and effective stress of the system. In this study, we experimentally studied drainage processes in unsaturated granular media as a model porous system. The distribution of saturated clusters is analysed by an optical imaging method under different drainage conditions, in which pore-scale information from Voronoi and Delaunay tessellation was used to characterise the topology of saturated cluster distributions. By employing statistical analyses, the observed spatial and temporal information of multiphase flow and fluid entrapment in porous media are described. The results indicate that the distributions of both the crystallised cell size and pore size are positively correlated to the spatial and temporal distribution of saturated cluster sizes. The saturated cluster size is found to follow a lognormal distribution, in which the generalised Bond number (𝐵𝑜 * ) correlates negatively to the scale parameter (𝜇) and positively to the shape parameter (𝜎). These findings can be used to connect pore-scale behaviour with overall hydro-mechanical characteristics in partially saturated porous media, using both the degree of saturation and generalised Bond number.

Introduction

Multiphase fluid flow in porous media is an important phenomenon with a broad range of applications, including risk assessment of landslides, biophysical systems, and groundwater hydrology [START_REF] Dai | Landslide risk assessment and management: an overview[END_REF][START_REF] Kueper | An overview of immiscible fingering in porous media[END_REF][START_REF] Saba | Effect of groundwater flow dimensionality on mass transfer from entrapped nonaqueous phase liquid contaminants[END_REF][START_REF] Seymour | Anomalous Fluid Transport in Porous Media Induced by Biofilm Growth[END_REF][START_REF] Simmons | On a test case for density-dependent groundwater flow and solute transport models: The Salt Lake Problem[END_REF]. Predicting the properties of such systems usually requires pore-scale information including the topological features and interactions between intrinsic phases. In a multiphase flow system, as the wetting phase is withdrawn and replaced by the nonwetting phase (gas), the interface of wetting and non-wetting fluid phases tend to recede to the outlet as a continues moving front. The displacement dynamics of the interface masks the complex interplay of viscous, capillary, and gravitational forces in the system and attracted extensive studies [START_REF] Aker | Viscous Stabilization of 2D Drainage Displacements with Trapping[END_REF][START_REF] Bear | On the Flow of Two Immscible Fluids in Fractured Porous Media[END_REF][START_REF] Bryant | Prediction of relative permeability in simple porous media[END_REF][START_REF] Lenormand | Invasion Percolation in an Etched Network: Measurement of a Fractal Dimension[END_REF][START_REF] Måløy | Viscous Fingering Fractals in Porous Media[END_REF][START_REF] Savani | A Monte Carlo Algorithm for Immiscible Two-Phase Flow in Porous Media[END_REF][START_REF] Schwarze | Rheology of weakly wetted granular materials: a comparison of experimental and numerical data[END_REF][START_REF] Youngs | The drainage of liquids from porous materials[END_REF].

Models used to study drainage processes in porous media include capillary tubes [START_REF] Dullien | Hydraulic continuity of residual wetting phase in porous media[END_REF][START_REF] Morrow | Irreducible wetting-phase saturations in porous media[END_REF][START_REF] Youngs | Drainage to a water table analysed by the Green-Ampt approach[END_REF] and two dimensional Hele-Shaw cells [START_REF] Haines | Studies in the physical properties of soil. V. The hysteresis effect in capillary properties, and the modes of moisture distribution associated therewith[END_REF][START_REF] Méheust | Interface scaling in a two-dimensional porous medium under combined viscous, gravity, and capillary effects[END_REF]. At high flow rates or in the absence of gravity, viscous forces tend to destabilise the continuous front into narrow fingers against the stabilising effect of gravitational force, which are the prerequisite conditions for the generation of limited saturated zones behind the continuous front [START_REF] Dullien | Hydraulic continuity of residual wetting phase in porous media[END_REF][START_REF] Måløy | Viscous Fingering Fractals in Porous Media[END_REF][START_REF] Méheust | Interface scaling in a two-dimensional porous medium under combined viscous, gravity, and capillary effects[END_REF][START_REF] Morrow | Irreducible wetting-phase saturations in porous media[END_REF][START_REF] Youngs | The drainage of liquids from porous materials[END_REF][START_REF] Youngs | Drainage to a water table analysed by the Green-Ampt approach[END_REF]. However viscous entrapped saturated zones left behind the moving front were neglected in these studies. As such, understanding interface dynamics between different liquids in drainage processes requires a better description of multiphase flow in porous media.

Fluid-fluid interface displacement in porous media along the front can be described at the pore scale as individual jumps at pore throats [START_REF] Haines | Studies in the physical properties of soil. V. The hysteresis effect in capillary properties, and the modes of moisture distribution associated therewith[END_REF][START_REF] Lenormand | Capillary fingering: Percolation and fractal dimension[END_REF]. Based on this, different approaches have been taken to explain fluid front morphologies, including a coupled capillaries model [START_REF] Hoogland | The formation of viscous limited saturation zones behind rapid drainage fronts in porous media[END_REF], which is established by a pair of neighbouring capillaries with different radii; a model consisting of small drops sliding down from different substrates [START_REF] Podgorski | Corners, Cusps, and Pearls in Running Drops[END_REF]; and the use of sliding liquid elements [START_REF] Su | Small-Scale Features of Gravity-Driven Flow in Unsaturated Fractures[END_REF]. The competition between viscous, capillary, and gravitational forces during drainage in a two-dimensional synthetic porous medium was considered by [START_REF] Auradou | Competition between correlated buoyancy and uncorrelated capillary effects during drainage[END_REF] and Meheust et al. (2002) to deduce interface movements. For cases of two or more immiscible fluid phases in porous media, variations of displacement structure are governed by the physical properties of the involved fluids and porous media [START_REF] Meakin | Scaling properties for the surfaces of fractal and nonfractal objects: An infinite hierarchy of critical exponents[END_REF][START_REF] Paterson | Diffusion-Limited Aggregation and Two-Fluid Displacements in Porous Media[END_REF]. These properties include the flow rates, wetting properties, viscosity ratio and density of the involved fluids, and pore structure [START_REF] Auradou | Competition between correlated buoyancy and uncorrelated capillary effects during drainage[END_REF][START_REF] Podgorski | Corners, Cusps, and Pearls in Running Drops[END_REF][START_REF] Su | Small-Scale Features of Gravity-Driven Flow in Unsaturated Fractures[END_REF].

The formation of drainage fronts (i.e., the interfaces between wetting and non-wetting fluids) and displacement dynamics, shaped by the aforementioned properties, determine the residual liquid phase distribution in partially wetted porous media and affect the relative permeability of the correlated unsaturated regions [START_REF] Auradou | Competition between correlated buoyancy and uncorrelated capillary effects during drainage[END_REF][START_REF] Scherer | New methods to measure liquid permeability in porous materials[END_REF]). Recent experimental and numerical studies concerning drainage processes in porous media tend to focus mainly on the evolution of viscous fingering at various flow rates and flow velocity distributions along the continuous moving front [START_REF] Alim | Local Pore Size Correlations Determine Flow Distributions in Porous Media[END_REF][START_REF] Holzner | Intermittent Lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF][START_REF] Lenormand | Capillary fingering: Percolation and fractal dimension[END_REF][START_REF] Matyka | Power-exponential velocity distributions in disordered porous media[END_REF][START_REF] Vishnudas | A comprehensive numerical study of immiscible and miscible viscous fingers during chemical enhanced oil recovery[END_REF][START_REF] Zhao | Wettability control on multiphase flow in patterned microfluidics[END_REF]. However, studies on the morphological characteristics of liquid phase entrapment behind the front remain relatively scarce. Accurate predictions of the fraction of liquid trapped behind the drainage front and its distribution are essential in order to quantify the relative permeability and the effective stress of unsaturated soils and rocks [START_REF] Blunt | Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow[END_REF][START_REF] Simmons | On a test case for density-dependent groundwater flow and solute transport models: The Salt Lake Problem[END_REF].

For low flow rate conditions, fluid invasion events along the displacement front occur pore by pore, and can be predicted by the classic percolation theory [START_REF] Gouyet | Fractal structure of diffusion and invasion fronts in three-dimensional lattices through the gradient percolation approach[END_REF][START_REF] Hunt | Percolation theory for flow in porous media[END_REF]. However, for high flow rate cases, pores are invaded simultaneously, exhibiting large pressure bursts with liquid ejection. This phenomenon is manifested primarily through the complex interplay of three forces within the system, namely gravitational force, viscous force and capillary force [START_REF] Or | Scaling of capillary, gravity and viscous forces affecting flow morphology in unsaturated porous media[END_REF]. To quantify the contribution of these three forces to viscous unstable flows in porous media, two typical dimensionless parameters, i.e., Capillary number (𝐶𝑎 = 𝜇𝑣𝑎 * 𝛾𝑘)

and Bond number (𝐵𝑜 = ∆𝜌𝑔′𝑎 * 𝛾), are defined [START_REF] Al-Fossail | Correlation between capillary number and residual water saturation[END_REF][START_REF] Alvarez | A non-gradient based algorithm for the determination of surface tension from a pendant drop: Application to low Bond number drop shapes[END_REF][START_REF] Birovljev | Gravity invasion percolation in two dimensions: Experiment and simulation[END_REF][START_REF] Blom | How to Include the Capillary Number in Gas Condensate Relative Permeability Functions[END_REF][START_REF] Gawlinski | Continuum percolation in two dimensions: Monte Carlo tests of scaling and universality for non-interacting discs[END_REF][START_REF] Méheust | Interface scaling in a two-dimensional porous medium under combined viscous, gravity, and capillary effects[END_REF][START_REF] Or | Scaling of capillary, gravity and viscous forces affecting flow morphology in unsaturated porous media[END_REF], where 𝜇 is the viscosity of the wetting fluid, 𝑣 the mean front velocity, 𝛾 the surface tension, and 𝑘 the permeability of the porous medium, ∆𝜌 the density difference between the wetting and non-wetting fluids, 𝑔′ the component of gravitational acceleration and 𝑎 is the typical pore size (more details can be found in [START_REF] Crisp | Direct measurement of pore-size distribution on artificial and natural deposits and prediction of pore space accessible to interstitial organisms[END_REF]. It is clear that the Capillary number describes the ratio of viscous to capillary forces, while the Bond number indicates the ratio of gravitational to capillary forces. Following Meheust et al. (2002), these two dimensionless parameters can be combined as the generalised Bond number, 𝐵𝑜 * = 𝐵𝑜 -𝐶𝑎 [START_REF] Méheust | Interface scaling in a two-dimensional porous medium under combined viscous, gravity, and capillary effects[END_REF][START_REF] Moebius | Pore scale dynamics underlying the motion of drainage fronts in porous media[END_REF], to quantify the transition between different flow regimes, i.e., stable and unstable flow. Specifically, positive values of 𝐵𝑜 * indicate stable moving fronts while negative values of 𝐵𝑜 * mark unstable liquid fronts.

The flow-rate dependency of the displacement front morphology raises the following question: how does the flowrate affect the liquid phase fraction trapped behind a moving drainage front? [START_REF] Wildenschild | Quantitative analysis of flow processes in a sand using synchrotron-based X-ray microtomography[END_REF] used Xray tomography to image air invasion into the saturated sandy soil and observed that large values of residual saturations occurred for high flow rates. The same trend was observed in glass bead micro-models in Løvoll's study (2011). The observed increase in residual saturation behind the displacement front for higher flow rates occurs due to the preferential drainage of large pores bypassing and entrapping liquid phase in non-invaded volumes [START_REF] Hoogland | The formation of viscous limited saturation zones behind rapid drainage fronts in porous media[END_REF]. It should be noted that liquid entrapment behind a rapidly moving front is not permanent, rather this liquid drains through a redistribution flow, albeit at a rate significantly lower than that of the primary moving front. This redistribution and the dominant liquid flow both affect the formation and evolution of these discontinuous entrapped liquid phase, which are termed locally saturated clusters. However, the extent and dynamics of locally saturated clusters are ripe for further exploration. An improved understanding of the behaviour drained and trapped liquid fractions in multiphase flow in porous media is of great value towards diverse engineering applications.

This study aims to establish a link between the complex mechanisms inside porous media and the morphological characteristics of a discontinuously entrapped liquid phase (i.e., saturated clusters). The properties of "viscous delayed" regions are studied in this work by establishing an experimental system where saturated clusters are obtained by the forced drainage of water from a porous medium. In this system, the gravitational and viscous forces can be tuned by changing tilting angles and drainage rates, and the competition of capillary, gravity and viscous forces can be examined by introducing the generalised Bond number (𝐵𝑜 * ) to mark the onset of unstable flows. In the present paper we describe an experimental setup for the acquisition of entrapped liquid profiles and present data analysis methods to determine the influence of pore-scale structure on the spatial extent of saturated clusters and their size distributions under various drainage scenarios. Results are interpreted with a view towards their significance in partially saturated porous systems.

Methodology

Experimental setup

The present experimental system was designed for temporally and spatially resolved imaging of interface motion in micro-models consisting of glass beads. The key elements of the setup are shown in Fig. 1(a). A micromodel was constructed by a single layer of glass beads with a diameter of 2 mm that were randomly arranged between two parallel glass plates. The plates and beads were then slightly compressed by clamps, and the model was sealed with silicone glue along the bottom and side edges, leaving the top side open to the atmosphere. The observable area was 180 mm×270 mm, and the resulting porosity was about 𝑛 ≈ 0.4. After the model was installed, the outlet channel was connected to a high precision pump (Masterflex 77200-62, USA) to provide a stable volumetric injection and withdrawal flow rate. A digital camera (Nikon D3300) controlled by a computer over digiCamControl software was used to record the experiments at given time intervals based on different flow rates. [START_REF] Scheidegger | The Physics of Flow Through Porous Media[END_REF]. The contact angle for the container wall is 44.56° (std=0.18°), which has similar wettability to glass plates. The applied flow rates ranged from 1 ml/min to 20 ml/min for different tilt angles. These resulted in the mean front velocity determined from volumetric flow rates at the cross-section area taking into account the porosity, as shown in Table 1.

Image processing

The grayscale pixel-intensity histogram of raw images shows two distinguishable peaks corresponding to water and air phases. In order to separate the air-filled region from the saturated region, all grey levels under a given threshold were set to black ('wet') while the others were white ('dry'). After thresholding, the whole observation area was separated into different regions with glass beads assigned to wet or dry regions. Here only the saturated regions that are larger than the projection area of a typical glass bead (𝜋 mm 2 in this study) were considered. After removal of trapped liquid volumes with projection areas smaller than a single glass bead, the resulting black and white images provide the required geometrical information. To allow the visualisation of separate clusters of Small amounts of water are trapped in the void spaces between glass beads and the front and back walls. Based on the micro-model geometry, the pore space in our sample can be illustrated in two ways: a 3-D sphere network and a 2-D cylinder network. These entrapped liquid rings are discarded in this study.

Statistical analysis

In this section, we use the probability density function (PDF), 𝑓 𝑦 𝑤 , specifying the data vector 𝑦 under the given parameter 𝑤, which defines the shape of the distribution profile. Here the subscript letter for the vector element 𝑦 denotes statistically independent observations. According to the theory of probability, the PDF for 𝑦 = (𝑦 G … 𝑦 > ) can be expressed as a multiplication of their PDFs. The probability density function of the lognormal distribution considering two parameters, which describes the relative likelihood for the continuous random variable to take on a given value, is given as [START_REF] Aitchison | The lognormal distribution[END_REF])

𝑓 𝑋 𝜇, 𝜎 * = G *KL M N 𝑒𝑥𝑝 - QB NRS M *L M , (1)
where the scale parameter 𝜇 is interpreted as the logarithmic mean of the random variables, and the shape parameter σ is interpreted as the standard deviation of the random variables in a logarithmic scale. As the parameters change in value, different probability distributions are generated. In this case, parameter estimation, a procedure of finding parameter values of a model that best fits the data, is vital to the following analysis.

In order to find the statistical parameters that best fit the acquired experimental dataset, the maximum likelihood estimation (MLE) method stands out for its attributes of sufficiency, consistency and efficiency [START_REF] Myung | Tutorial on maximum likelihood estimation[END_REF].

By taking the product of the probability densities of the individual 𝑋 A (𝑖 = 1,2,3, … . 𝑛), we derived the likelihood function of the lognormal distribution for a series of 𝑋 A . The likelihood function 𝐿 𝜇, 𝜎 * 𝑋 is defined by reversing the roles of data vector and parameter vector:

𝐿 𝜇, 𝜎 * 𝑋 = 𝑓 𝑋 A 𝜇, 𝜎 * B A[G .
By taking the natural logarithm of the likelihood function and taking the gradient of 𝐿 with respect to 𝜇 and 𝜎 * while setting it equal to zero, we can find 𝜇 and 𝜎 * which maximise 𝐿 [START_REF] Ginos | Parameter estimation for the lognormal distribution[END_REF]. In addition, to measure the linear relationship between variables, the Pearson's correlation coefficient 𝑟 is employed. Moreover, the correlation coefficient is a rescaled covariance:

𝑟 N,] = ^_`(N,]) L a L b
, where cov(𝑋, 𝑌) is the sample covariance, and 𝜎 N and 𝜎 ] are sample standard deviations. After rescaling the covariance, it ranges between -1 and +1, where 1 is total positive linear correlation, 0 is no linear correlation, and -1 is total negative linear correlation (Lee Rodgers and Nicewander, 1988).

Results analysis

Characterisation of liquid occupied regimes

The flow of two immiscible fluids through a granular medium depends not only on the complex interplay between viscous, capillary and gravitational forces, but also the topological characteristics of the porous structure, here in particular the packing characteristics of the grains. The interaction between these forces and the pore structure yields a complex flow regime, shown in Fig. 3. It can be observed that there are four different liquid content forms, namely adsorbed, pendular, funicular and fully saturated regimes. The projection areas of these four regimes can be separated on the basis of their water content. In these four regimes, the water fraction increases from the top zoomed-in picture to the bottom one. Specifically, in the absorbed regime, liquid phase exists as thin films coating on each spheres. As the liquid phase increase, pendular bridges form between consecutive spheres. When further water is added, the liquid phase becomes continuous, while the air phase becomes discontinuous, marking the funicular regime. Finally, in the fully saturated regime, all voids are occupied by the liquid phase. In the following sections, we focus on the spatial and temporal evolution of the fully saturated regime. Though the information from other regimes may be practically useful, e.g., the percolation of the liquid bridges presented in the pore space, their study may require higher resolution at the local scale.

Spatial distribution of wetted cells

To look into the influence of packing characteristics on saturated cluster formation and evolution during drainage processes, the correlation between the temporal-spatial information of saturated clusters and packing characteristics is compared here. As can be seen in Fig. 4, the histograms in the first column (S r = 0.5) have higher probability densities than those in the second column (S r = 0.6) for the same sub windows in both directions, i.e., the lower saturation level has more wetted cells than the high saturation level at the same region in the model. This is due to the removal of the large fully saturated region behind the liquid front, and for smaller saturation levels larger accumulated saturated areas can be extracted. This can be explained by considering that saturated clusters form when the liquid front meets different conductive pathways during drainage processes, and therefore a sufficient distance along the flow direction is a prerequisite for the formation of such clusters. At lower saturation levels, a greater air content is present, more resistive pathways are bypassed and subsequently the accumulated area of viscous-detained liquid also increases. In addition, as 𝐵𝑜 * increases, the accumulated wetted cell area included in saturated clusters is reduced. This follows the aforementioned theory about the influence of 𝐵𝑜 * on geometry of liquid front in Section 2: for large negative values of 𝐵𝑜 * , the capillary fingers are widely spread along the liquid front during drainage processes; for large positive values of 𝐵𝑜 * , the liquid front is observed as a nearly flat geometry. For flat liquid fronts, gravity is always large enough to stabilise the displacement geometry leading to fewer saturated clusters trapped behind the liquid front and for unstable liquid fronts, more thin fingers can be observed, which are disconnected from the fully saturated region.

To look into the influence of packing characteristics on saturated cluster formation and evolution during drainage processes, we plot the number and probability distributions of wetted cells in Fig. 5. The wet cells are labelled with 1 and dry cells are labelled with 0. By collecting the amount of 1-labelled cells in an individual group divided by the total amount of the cells in a group, the wetted probability for each group is computed. The details can be found in Appendix A. In Fig. 5, we specified the bin space with 0.4 mm 2 and collected a list of bins and histogram heights of the values. The pink curves with empty circles in each image are computed from initially fully saturated images with different drainage conditions. The other coloured curves represent wetted cell area distributions within saturated clusters for different saturation levels associated with different drainage conditions. For these cases, the largest continuous saturated region below the liquid front is removed from every picture and only saturated clusters were left behind in each frame. By employing image processing operations on individual images, we collected the data stream from five cases with the same flow conditions and merged all sub-lists as the final step. To have a better visualisation of how the data is distributed, we depicted both the wetted cell area number distribution and wetted cell area probability distribution as functions of cell area. In Figs. 5(a) and (b), it can be seen that over half of the cells have areas in the range 4-6 mm 2 . By removing the projection area of individual glass beads (𝜋 mm 2 in this study), the void space in these cells is then found to be 1.72-5.72 mm 3 , which correlates to relatively close packed regions in the experimental sample. Table 2 shows the pore-scale size parameters and wetted probabilities for Voronoi cells computed from five pictures. The experiment shown here corresponds to a saturation level value of 0.5 and a generalized Bond number of -0.003. By grouping all cells with their individual face area, the corresponding dry/wet properties of every cell are also distributed into the same group. The relatively large values of wetted probability in Table 2 are most widely observed cells in the crystalized region. However, the smaller values are from cells with extremely large or small areas. This trend shows that in our porous medium, the most widely distributed cells have stronger liquid trapping ability than the extreme regions. Similarly, in Figs. 5(c) and (d), results have shown that more than half of the wetted cells are distributed in the relatively close packed region for all drainage conditions.

Saturated cluster distribution

To characterize the effect of drainage conditions on the evolution of saturated clusters, statistical analysis of saturated cluster geometries was conducted. The largest saturated cluster identified in each image, corresponds to the saturated region in front of the continuous liquid front, and is therefore excluded from the analysis. Here, we first look at the effective size of the saturated clusters, which is acquired by first collecting the number of pixels in a localized saturated region and then divided by the spatial resolution 30 pixels per pore (1 mm 2 ). To acquire sufficient data for statistical analysis, we performed five sets of experiments for each set of conditions. By grouping pictures with the same saturation, S r , and generalized Bond number, 𝐵𝑜 * , and rearranging the experimental data into multiple sub-lists by a given bin size, the saturated cluster distributions are obtained.

The insert in Fig. 6 grouping the dataset with unit width bins equal to 2.5 mm 2 , the dataset is separated into sub lists, based on which the probability density distribution is calculated. It is observed that the probability density tends to the peak value (0.15-0.2 / mm 2 ) around 4-8 mm 2 . As our porous medium is constructed by glass beads of 2 mm in diameter, the most frequent saturated clusters in the experimental cell can be interpreted as containing 1-2 glass beads. Different types of distributions have also been tested, while a lognormal distribution function, Eq. ( 1), was selected here to best represent the experimental data. The goodness-of-fit for three common fitting models are further compared with details provided in Appendix B.

As stated in Section 2.3, the lognormal density distribution function has two parameters, i.e., scale parameter 𝜇 and shape parameter 𝜎, which govern the shape of the distribution curve. To estimate these two parameters from the given observations, we use the maximum likelihood estimation (MLE) method, described in Section 2. For a given lognormal distribution, 𝜎 and 𝜇 are the mean and standard deviations of the variable's natural logarithm, respectively. The linear relations obtained for the lognormal density function parameters as functions of 𝐵𝑜 * are therefore correlated to the logarithmic mean and variation for saturated clusters sizes. In this case, as 𝐵𝑜 * increases, the logarithmic mean saturated cluster size decreases; meanwhile, the saturated cluster sizes are spread out over a wider range of values. This correlation can be used to predict the geometrical distributions of saturated clusters trapped behind a liquid front for known pore structures and drainage conditions.

Conclusion

In this study we experimentally evaluated the formation of a temporarily entrapped liquid phase and its subsequent evolution during drainage processes of wetted granular materials. In our experiments the competition among viscous, gravitational and capillary forces was controlled to achieve different conditions of drainage from the porous media (characterised by the generalised Bond number, Bo * ). Due to the inhomogeneous capillary threshold (governed by pore throat) in the porous medium and dominance of viscous effects, saturated clusters were trapped behind the moving front. The spatial and temporal distribution of saturated clusters can thus be observed, and the following conclusions can be made:

(1) The fraction of viscous limited pore bodies depends on the hydraulic properties of the resistive liquid phase, which is held in poorly connected pore clusters and small pore bodies due to the fact that gravitational force cannot overcome the effect of the viscous force. It is demonstrated that, under various drainage conditions and saturation stages, the statistical distribution of the liquid phase is determined by the intrinsic pore network characteristics. It is found that the pore-scale structure has a strong influence on the spatial extent of viscous delayed pores during a multiphase flow. The wetted cell size probability curve closely follows the global-scale cell size probability curve.

(2) The size of saturated clusters is found to follow a lognormal distribution, and the scale parameter 𝜇 and shape parameter 𝜎 negatively and positively correlate with Bo * , respectively. Due to the fact that a relatively strong linear dependence exists between the lognormal distribution parameters, 𝜇 and 𝜎, and generalised Bond number, Bo * , we are able to describe the cluster distribution at any given drainage conditions and saturation levels.

By extending the model towards three dimensional systems based on the explicit description of disordered pore networks, conclusions extracted from this study can be applied to a variety of situations in realistic geological scenarios involving multiphase flow. The ability to predict the presence of residual fluid phase left behind a moving liquid front on the basis of topological properties is of value in optimising water injection processes towards improved oil recovery [START_REF] Blunt | Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow[END_REF][START_REF] Vishnudas | A comprehensive numerical study of immiscible and miscible viscous fingers during chemical enhanced oil recovery[END_REF] and pollutant removal [START_REF] Saba | Effect of groundwater flow dimensionality on mass transfer from entrapped nonaqueous phase liquid contaminants[END_REF]. The observations of the present study further show how particle scale affects overall flow conditions, potentially facilitating improved formulation of the effective stress term in unsaturated soil mechanics [START_REF] Gan | A particle-water based model for water retention hysteresis[END_REF][START_REF] Likos | Effective Stress in Unsaturated Soil: Accounting for Surface Tension and Interfacial Area[END_REF][START_REF] Xu | Effective stress in soils under different saturation conditions[END_REF].

The variation of the capillary threshold here was limited, due to the topological homogeneity inherent to packings of monosized glass beads. For more complex systems, a quantitative prediction of the viscous limited saturated cluster distribution remains challenging. The results presented here provide an initial step towards a better understanding of drainage in porous media, suggesting further detailed modelling and experimental studies into saturated clusters formation and evolution.

  Each image was rescaled to 1000×1500 pixels, which corresponds to a spatial resolution of 30 pixels 1 mm 2 (~ the pore area). An electroluminescent panel provided sufficient background illumination for different shutter speeds. Both the light source and model could be tilted simultaneously to various inclination angles to control the gravitational component of the fluid flow in the porous medium. The component of gravity along the model 𝑔' is given by 𝑔 sin(𝛼), where 𝑔 is the gravitational acceleration and 𝛼 is the inclination angle. Water was used as the wetting fluid in all experiments, dyed with tracer (Cole-Parmer 00298-06, 2.6ml/100ml), and air was the non-wetting fluid invading from the open top edge. Experimental values of the wetting fluid were measured at 25°C. The surface tension between two fluids 𝛾 is 64.27 mN/m (standard deviation, std=0.15 mN/m), which was performed by a goniometer (Rame-hart, model 200 Standard Contact Angle Goniometer with DROPimage Standard) with the pendent drop method. The wetting fluid has a viscosity 𝜇 of 0.95 mPa.s (std=0.2 mPa.s) and a density 𝜌 of 998 kg/m3. The permeability of the medium 𝑘 was estimated from the Kozeny-Carman relation as 4.89×10 -9 m 2

  different sizes, analyses of morphological components have been carried out and colours were assigned to the localised saturated clusters, as shown Fig. 1(b) and Fig. 1(c). The inner rings in solid beads are shown in Fig. 1(d).

Fig. 1

 1 Fig. 1 Schematic illustration of the experimental setup. (a) The camera, electroluminescent panel and cell are all connected to a common frame that can be tilted. The liquid outlet on the bottom is connected to a syringe pump that injects and withdraws the liquid. (b) A detailed zoomed-in section from the raw image extracted from the camera showing the structure of the porous network, the non-wetting phase (green) and wetting phase (dark red). (c) The same section after applying image processing. Individual saturated clusters are coloured randomly to allow the visualisation of separated clusters for different sizes. (d) Pendular rings between the glass bead and the front glass and schematic of a glass bead with pendular rings.

  Fig. 2 (a) Schematics of Voronoi lattice (blue dashed lines) and Delaunay triangulation (red solid lines). (b) Map of Voronoi mesh for an individual fully saturated image. The bar legend is scaled in different gradient, ranging from smallest (𝑒 >AB = 0.103) to largest (𝑒 >?@ = 7.169) Voronoi cell area. (c) Map of Delaunay mesh for an individual fully saturated image. The bar legend is scaled in different gradient, ranging from smallest (𝑒 >AB C = 0.103) to largest (𝑒 >?@ C = 11.162) cell void ratio.

Fig. 3

 3 Fig. 3 An example experimental image revealing different saturation regimes

Fig. 4

 4 Fig. 4 Histogram showing the distributions of wetted cells in the model for three experiments in a Voronoi mesh: (a) and (b) the wetted cell density inside sub windows in the top-bottom direction for different saturation levels, as a function of the distance 𝑦 to the outlet; (c) and (d) the wetted cell density in the perpendicular direction for different saturation levels. The insets show the corresponding sub window orientations with respect to the model (made ≈ 4 times taller than the ones used in our analysis, to aid visualisation.

Fig. 5

 5 Fig. 5 Wetted Voronoi cell and Delaunay Triangle (within saturated clusters) area distributions for six drainage conditions. The generalized Bond number in this case is -0.002, which could be referred to the experiments in Table. 1. The wetted cell number and probability distribution are calculated using Voronoi, (a) and (b), and Delaunay tesselations, (c) and (d).

  (a) shows some observed elementary units of saturated clusters containing one to three glass beads during drainage processes. (I) In the subfigure (1) in Fig.6(a), the central glass bead is fully coated and its neighbouring pores are occupied by the dyed liquid. In this subfigure, the saturated cluster is defined as the combination of the central glass bead and all its neighbouring pore spaces. (II) For saturated clusters with two glass beads inside, as shown in the subfigure (2) in Fig.6(a), the volume of liquid in the inter-particle pore depends on the distance between two neighbouring beads, which can be quantified by the projection area of saturated clusters. (III) As for saturated clusters containing three beads, we assume that the pore structures in subfigure (3) in Fig.6(a) are the same. However, owing to the fact that the capillary effect in granular media can be influenced by multiple factors, including pore size, flow rate and the gravitational force, the projection area of these two saturated clusters are different. These factors can be quantified by the generalized Bond number 𝐵𝑜 * . In what follows, saturated clusters are extracted for different values of S r and 𝐵𝑜 * . By collecting all the data at each stage

Fig. 6

 6 Fig. 6 (a) The size distribution of saturated clusters for different liquid saturation levels with 𝐵𝑜 * =-0.003. (b) Saturated cluster distribution at the same degree of saturation (𝑆 l =0.5) with 𝐵𝑜 * =-0.003 (Dots: the experimental data, and Line: lognormal distribution by maximum likelihood estimation). (c) Scale parameter 𝜇 and (d) shape parameter 𝜎 as a function of generalised Bond number, 𝐵𝑜 * , for different saturation levels.

.

  3. The resulting log-normal distribution can then be plotted by substituting the specific parameters acquired from the MLE method, as inFig. 6(b). This figure illustrates the saturated cluster probability distribution as a function of area for different values of the generalized Bond number 𝐵𝑜 * at a specific saturation levels S r , and good agreement can be found between the MLE-fitted distribution and experimental data. Due to the removal of saturated regions that are smaller than a projection area of a typical glass bead (π mm 2 in this study), the corresponding fitting curve in Fig.6(b) starts at the point: (π, 0).After collecting parameters 𝜇 and 𝜎 for all series in the experiments, we applied Pearson's correlation coefficient method to quantify the relative linear relationship between them. To compute the correlation coefficient for both 𝜎 and 𝜇 as a function of 𝐵𝑜 * , we substitute the acquired values of these two parameters with respect to 𝐵𝑜 * into The correlated values for both parameters are shown inFigs. 6(c) and 6(d). In Fig.6(c), the correlation coefficient of 𝜇 and 𝐵𝑜 * is -0.8039, implying that all data points lay close to a line for which 𝜇 decreases as 𝐵𝑜 * increases. On the contrary, the value of correlation coefficient of 𝜎 and 𝐵𝑜 * is 0.6608, indicating that 𝜎 increases as 𝐵𝑜 * increases, as shown in Fig.6(d).

Table 1 .

 1 Evaluation of mean front velocity and generalized bond number (𝐵𝑜 * ) used for experiments.

	Tilt Angle (°)			0				5		10		15
	Drainage Rate (mL/min)	1	2	5	10	20	1	5	20	5	20	10
	Mean Front
	Velocity	0.116	0.231	0.579	1.157	2.315	0.116	0.579 2.315 0.579 2.315 1.157
	(mm/s)
	𝐵𝑜 * (-)	-1.556 ×10 -4

Table 2

 2 Pore-scale size parameters and wetted probabilities at different individual face area of Voronoi tessellation for 𝐵𝑜 * = -0.003, 𝑆 h =0.5. plotting their probability density distribution, it is possible to determine appropriate distribution parameters that represent the data set, and relate these to 𝐵𝑜 * .

	𝐵𝑜 * = -0.003

and

Acknowledgements

The authors are grateful for the financial support from Australian Research Council (Projects DE130101639 and DP170102886) and The University of Sydney SOAR Fellowship.

Appendix A: Voronoi and Delaunay tessellation analysis

The procedures of using Voronoi and Delaunay tessellation in our statistical analysis are shown in detail below:

1. We first compute Delaunay triangles and Voronoi cells associated with the network using the image processing method mentioned in Section 3. After numbering and collecting the vertex coordinates of every cell in both mesh regions, we are then able to create region-member functions associated to each cell, for instance, the 𝑖-th region-member function for the 𝑖-th cell.

2. Then, import all pictures representing different saturation levels and drainage conditions into our calculation. These pictures are then be performed with morphological image processes and transformed into a black and white image, while the black and white pixels stand for dry and wetted regions, respectively. The black islands inside the white area are the wetted clusters trapped in the air during drainage. To allocate our focus on the saturated clusters behind the liquid front, we removed the largest continuous black island in each black and white picture, leaving other relatively small wetted clusters inside the pictures.

3. With the help of the built-in function (PixelValuePositions) in Mathematica, we then extracted all the black pixel positions in every picture. By repeatedly applying different black points into the 𝑖th regionmember function, we returned to the results of true or false for whether these points are region members of the 𝑖th cell. This step gave us the results of the amount of points in the 𝑖th cell and we were then able to know the amount of black points in every cell by repeatedly employing this step on all rest cells in an individual picture. As a result, we were able to return to a list of number of black points in each cell aligned with the cell number. 4. By comparing the ratio of black pixel number to white pixel number, we classified the cells into two groups: dry and wet. After filtering the uncorrelated information from the aforementioned steps, we finally came to a list giving the area and wet (1) / dry (0) properties for every cell in an individual picture.

To enlarge our dataset, we also did five experiments under same conditions and merged all five lists to perform further statistical analysis.

Appendix B: Goodness-of-fit test results

Table . B1 shows two widely used tests, the Pearson's χ * and Anderson-Darling methods to test the validity of three assumed distribution models (Lognormal distribution, Weibull distribution and Gamma distribution). We applied a particular dataset (𝑆 h =0.5; 𝐵𝑜 * = -0.003) to the two test methods and listed the resulting values in Table . B1. These results are used to verify the validity of lognormal distribution which has been selected on the basis of other prior considerations.