Camille Carvalho

Shilpa Khatri

Arnold D Kim

Asymptotic analysis for close evaluation of layer potentials

Keywords: Boundary integral equations, Laplace's equation, Layer potentials, Nearly singular integrals, Close evaluations

We study the evaluation of layer potentials close to the domain boundary. Accurate evaluation of layer potentials near boundaries is needed in many applications, including fluid-structure interactions and near-field scattering in nano-optics. When numerically evaluating layer potentials, it is natural to use the same quadrature rule as the one used in the Nyström method to solve the underlying boundary integral equation. However, this method is problematic for evaluation points close to boundaries. For a fixed number of quadrature points, N , this method incurs O(1) errors in a boundary layer of thickness O(1/N ). Using an asymptotic expansion for the kernel of the layer potential, we remove this O(1) error. We demonstrate the effectiveness of this method for interior and exterior problems for Laplace's equation in two dimensions.

Introduction

Boundary integral equation methods are useful for solving boundary value problems for linear, elliptic, partial differential equations [START_REF] Guenther | Partial Differential Equations of Mathematical Physics and Integral Equations[END_REF][START_REF] Kress | Linear Integral Equations[END_REF]. Rather than solving the partial differential equation directly, one represents the solution as a layer potential, an integral operator applied to a density. The density is the solution of an integral equation on the boundary of the domain that includes the prescribed boundary data. This formulation offers several advantages for the numerical solution of boundary value problems. First, the dominant computational cost is from the integral equation on the boundary whose dimension is lower than that of the domain. Second, this boundary integral equation can be solved to very high order using Nyström methods [START_REF] Atkinson | The Numerical Solution of Integral Equations of the Second Kind[END_REF][START_REF] Delves | Computational Methods for Integral Equations[END_REF]. Finally, the solution, given by this layer potential, can be evaluated anywhere in the domain without restriction to a particular mesh. For these reasons, boundary integral equations have found broad applications, including in fluid mechanics and electromagnetics.

One particular challenge in using boundary integral equation methods is the so-called close evaluation problem [START_REF] Barnett | Evaluation of layer potentials close to the boundary for Laplace and Helmholtz problems on analytic planar domains[END_REF][START_REF] Helsing | On the evaluation of layer potentials close to their sources[END_REF]. Since a layer potential is an integral over the boundary, it is natural to evaluate it numerically using the same quadrature rule used in the Nyström method to solve the boundary integral equation. In that case, we say that the layer potential is evaluated using its native quadrature rule. Numerical evaluation of the layer potential using its native quadrature rule inherits the high order accuracy associated with solving the boundary integral equation, except for points close to the boundary. For these close evaluation points, the native quadrature rule produces an O(1) error. This O(1) error is due to the sharply peaked kernel of the layer potential leading to its nearly singular behavior.

There are several problems that require accurate layer potential evaluations close to the boundary of the domain. For example, modeling of micro-organisms swimming, suspensions of droplets, and blood cells in Stokes flow use boundary integral methods [START_REF] Smith | A boundary element regularized Stokeslet method applied to cilia-and flagella-driven flow[END_REF][START_REF] Barnett | Spectrally accurate quadratures for evaluation of layer potentials close to the boundary for the 2D Stokes and Laplace equations[END_REF][START_REF] Marple | A fast algorithm for simulating multiphase flows through periodic geometries of arbitrary shape[END_REF][START_REF] Keaveny | Applying a second-kind boundary integral equation for surface tractions in Stokes flow[END_REF]. The key to these problems is the accurate computation of velocity fields or forces close to the boundary as these quantities provide the physical mechanisms leading to locomotion and other phenomena of interest. Another example is in the field of plasmonics [START_REF] Maier | Plasmonics: Fundamentals and Applications[END_REF], where one seeks to gain control of light at the sub-wavelength scales for applications such as nano-antennas [START_REF] Akselrod | Probing the mechanisms of large Purcell enhancement in plasmonic nanoantennas[END_REF][START_REF] Novotny | Antennas for light[END_REF] and sensors [START_REF] Mayer | A label-free immunoassay based upon localized surface plasmon resonance of gold nanorods[END_REF][START_REF] Sannomiya | In situ sensing of single binding events by localized surface plasmon resonance[END_REF]. Surface plasmons are subwavelength fields localized at interfaces between nano-scale metal obstacles and their surrounding dielectric background medium. Thus, these problems require accurate computation of electromagnetic fields near interfaces. These problems and others motivate the need to address the close evaluation problem.

The close evaluation problem for layer potentials has been studied previously for Laplace's equation. For example, Beale and Lai [START_REF] Beale | A method for computing nearly singular integrals[END_REF] have studied this problem in two dimensions by first regularizing the nearly singular kernel and then adding corrections for both the discretization and the regularization. The result of this approach is a uniform error in space. This method has been extended to threedimensional problems [START_REF] Beale | A simple method for computing singular or nearly singular integrals on closed surfaces[END_REF]. Helsing and Ojala [START_REF] Helsing | On the evaluation of layer potentials close to their sources[END_REF] have studied the Dirichlet problem in two dimensions by combining a globally compensated quadrature rule along with interpolation to achieve very accurate results over all regions of the domain. Barnett [4] has also studied this problem in two dimensions. In that work, Barnett has established a bound for the error associated with the periodic trapezoid rule. We make use of this result in our work below. To address the O(1) error in the close evaluation problem, Barnett has used surrogate local expansions with centers placed near, but not on, the boundary. This new method, called quadrature by expansion (QBX), leads to very accurate evaluations of the layer potential close to the boundary. Further error analysis of this method and extensions to evaluations on the boundary for the Helmholtz equation is presented in Klöckner et al [START_REF] Klöckner | Quadrature by expansion: A new method for the evaluation of layer potentials[END_REF]. For the special case of rectangular domains, Fikioris et al [START_REF] Fikioris | Strongly and uniformly convergent Green's function expansions[END_REF][START_REF] Fikioris | Strongly convergent Green's function expansions for rectangularly shielded microstrip lines[END_REF] have addressed the close evaluation problem by removing problematic terms from the explicit eigenfunction expansion of the fundamental solution.

Here, we develop a new method to address the close evaluation problem. We first determine the asymptotic behavior of the sharply peaked kernel and then use that to approximate the layer potential. Doing so relieves the quadrature rule from having to integrate over a sharp peak. Instead, the quadrature rule is used to correct the error made by this approximation. This new method is accurate, efficient, and easy to implement.

In this paper, we study the close evaluation problem in two dimensions for Laplace's equation. We use a Nyström method based on the periodic trapezoid rule to solve the boundary integral equation. We study the double-layer potential for the interior Dirichlet problem and the single-layer potential for the exterior Neumann problem. For both of these problems, we introduce an asymptotic expansion for the sharply peaked kernel of the layer potential for close evaluation points, which is the main cause for error. Using the Fourier series of this asymptotic expansion, we compute its contribution to the layer potential with spectral accuracy. Through several examples, we show that this asymptotic method effectively reduces errors in the close evaluation of layer potentials.

The remainder of this paper is as follows. In Section 2 we study in detail the illustrative example of the interior Dirichlet problem for a circular disk. For this problem, we obtain an explicit error when using the periodic trapezoid rule to evaluate the double-layer potential. This error motivates the use of an asymptotic expansion for the sharply peaked kernel in the general method we develop in Section 3 to evaluate the double-layer potential for the interior Dirichlet problem. In Section 4, we extend this method to the single-layer potential for the exterior Neumann problem. We discuss the general implementation of these methods in Section 5. Section 6 gives our conclusions.

Illustrative example: interior Dirichlet problem for a circular disk

We first study the close evaluation problem for ∆u = 0 in D = {r < a}, (2.1a)

u = f on ∂D = {r = a}. (2.1b)
For this problem, we compute an explicit error when applying an N -point periodic trapezoid rule (PTR N ). This error reveals the key factors leading to the large errors observed in the close evaluation problem. Moreover, this analysis provides the motivation for the general asymptotic method.

It is well understood that the solution of (2.1) is given by Poisson's formula [START_REF] Strauss | Partial Differential Equations[END_REF]. Here, we seek the solution as the double-layer potential [START_REF] Kress | Linear Integral Equations[END_REF],

u(x) = 1 2π |y|=a n y • (x -y) |x -y| 2 µ(y)dσ y . (2.2)
Here, x ∈ D denotes the evaluation point, y ∈ ∂D denotes the variable of integration, n y denotes the unit, outward normal at y, and dσ y denotes a differential boundary element. The density, µ(y), satisfies the following boundary integral equation,

- 1 2 µ(y) - 1 4πa |y |=a µ(y )dσ y = f (y), (2.3)
from which we determine that

µ(y) = 1 2πa |y |=a f (y )dσ y -2f (y). (2.4)
To numerically evaluate (2.2), we substitute the parameterization, x = (r cos t * , r sin t * ), and y = (a cos t, a sin t), with r < a, and t * , t ∈ [0, 2π], and obtain

u(r, t * ) = 1 2π 2π 0 ar cos(t -t * ) -a 2 a 2 + r 2 -2ar cos(t -t * ) µ(t)dt = 1 2π 2π 0 K(t -t * )µ(t)dt. (2.5)
Using PTR N with points t j = 2πj/N for j = 1, • • • , N , to approximate (2.5), we obtain

u(r, t * ) ≈ U N (r, t * ) = 1 N N j=1 K 2πj N -t * µ 2πj N . (2.6)
The error, U N (r, t * )-u(r, t * ), is not uniform in D. In particular, U N is very accurate for evaluation points far away from the boundary. On the other hand, it is O(1) when r ∼ a. The reason for this large error is due to the kernel, K(t -t * ). In Fig. 1, we show plots of K as a function of t -t * . The left plot of Fig. 1 shows that K becomes sharply peaked about t = t * when r ∼ a. We do not evaluate the double-layer potential on the boundary. Nonetheless, because K becomes sharply peaked as r → a, we say that the double-layer potential is nearly singular. The right plot of Fig. 1 shows that the piecewise linear approximation of K used by PTR 128 will grossly overestimate the magnitude of the double-layer potential for r = 0.99a. It is this error that leads to the O(1) errors produced by PTR N for close evaluation points. Barnett [4] has shown that this error is O(1) for a -r = O(1/N ). In light of this result, we say that the error exhibits a boundary layer of thickness O(1/N ) in which it undergoes rapid growth.

In the limit as N → ∞, PTR N converges because the boundary layer vanishes at a rate of O(1/N ). However, that is not the limit we consider here. Rather, we study the limit as the evaluation point approaches the boundary with N fixed. For that case, PTR N is unable to accurately capture the sharp peak of the kernel about t = t * that forms as r → a.

Using the error associated with PTR N [START_REF] Davis | On the numerical integration of periodic analytic functions[END_REF], we find U N defined in (2.6) satisfies

U N (r, t * ) = u(r, t * ) + ∞ l=-∞ l =0 p[lN ], (2.7)
where

p[k] = 1 2π 2π 0 K(t -t * )µ(t)e -ikt dt. (2.8)
The error in (2.7) is aliasing of high frequencies. To determine p[k] explicitly, we use the Fourier series representation of the kernel,

K(t -t * ) = ar cos(t -t * ) -a 2 a 2 + r 2 -2ar cos(t -t * ) = - 1 2 - 1 2 ∞ m=-∞ r a |m| e im(t-t * ) , (2.9)
and of the density

µ(t) = ∞ n=-∞ μ[n]e int , (2.10)
to find

K(t -t * )µ(t) = - 1 2 ∞ n=-∞ μ[n]e int - 1 2 ∞ m=-∞ r a |m| e -imt * ∞ n=-∞ μ[n]e i(m+n)t . (2.11)
Substituting (2.11) into (2.8), and rearranging terms, we find that

p[k] = -μ[k] - 1 2 ∞ m=1 r a m e imt * μ[k + m] + e -imt * μ[k -m] .
(2.12)

We now obtain the error in using PTR N to evaluate the double-layer potential by substituting (2.12) into (2.7) which yields

E N (r, t * ) = U N (r, t * )-u(r, t * ) = ∞ l=1 -μ * [lN ] - 1 2 ∞ m=1 r a m μ * [lN -m]e imt * + μ * [lN + m]e -imt * + ∞ l=1 -μ[lN ] - 1 2 ∞ m=1 r a m μ[lN + m]e imt * + μ[lN -m]e -imt * . (2.13)
Here, we have assumed

µ is real, so μ[-k] = μ * [k]
, where [•] * denotes complex conjugation. Suppose we have chosen N to be large enough so that μ[lN ] 1 for l > 0. For that case, only terms in (2.13) proportional to μ[lN -m] will substantially contribute to the error. By neglecting the other terms, we obtain

E N (r, t * ) ∼ - ∞ l=1 ∞ m=1 r a m [Re{μ[lN -m]} cos(mt * ) + Im{μ[lN -m]} sin(mt * )] .
(2.14) Equation (2.14) is the asymptotic error made by PTR N . The key point is that when N is fixed, this error is not uniform for r ∈ [0, a). When r a, we see that the error is much smaller than when r ∼ a. Consider the specific case in which µ = 1, so that μ[0] = 1 and μ[k] = 0 for all k = 0. For that case, (2.14) simplifies to In practice, we would like to set N based on the resolution required to solve the boundary integral equation. It is neither desirable nor practical to increase N just to reduce aliasing in the evaluation of the double-layer potential. In light of this, we make the following observations.

E N (r, t * ) ∼ r a 2N -r a N cos(N t * ) 1 + r a 2N -2 r a N cos(N t * ) . (2.15) According to (2.15), |E N (a(1 -ε), t * )| = O((1 -ε) N ) = O(e -εN
• Equation (2.13) gives the error incurred by PTR N to approximate the double-layer potential. This error is due to aliasing. Equation (2.14) gives the asymptotic approximation of this error when the N -point grid sufficiently samples µ.

• The aliasing error is not uniform with respect to r. For the case in which µ = 1, the asymptotic error simplifies to (2.15). From this result, we find that the error grows rapidly and becomes O(1) in a boundary layer of thickness O(1/N ) near the boundary. This boundary layer is shown in Fig. 2.

• For points within the boundary layer, the sharply peaked kernel causes aliasing due to insufficient resolution. Fig. 1 shows how the sharp peak of the kernel when r/a = 0.99 is under-resolved on the grid for PTR 128 .

Alternatively, by substituting (2.9) and (2.10) into (2.5), we obtain

u(r, t * ) = ∞ n=-∞ K * [n]μ[n]e -int * ≈ N/2-1 n=-N/2 K * [n]μ[n]e -int * . (2.16) Since the coefficients, μ[n] for n = -N/2, • • • , N/2 -1
, can be computed readily using the Fast Fourier Transform, we introduce the truncated sum as an approximation in (2.16). The decay of μ[n] controls the error of this approximation. Therefore, choosing N to accurately solve the boundary integral equation yields a spectrally accurate approximation of the double-layer potential.

For general problems, we do not know K[n] explicitly as we do here. Instead, we compute an asymptotic expansion of the sharply peaked kernel. We determine the Fourier coefficients of this asymptotic expansion explicitly. Hence, we evaluate its contribution to the double-layer potential using an approximation like the one in (2.16). By removing the kernel's sharp peak in this way, we are left with a smooth function to integrate using the PTR N . We present this method to evaluate the double-layer potential for the interior Dirichlet problem for Laplace's equation in Section 3, and the single-layer potential for the exterior Neumann problem for Laplace's equation in Section 4.

Double-layer potential for the interior Dirichlet problem

u = f on ∂D, (3.1b)
with f an analytic function. We seek u as the double-layer potential [START_REF] Kress | Linear Integral Equations[END_REF],

u(x) = 1 2π ∂D K(x, y)µ(y)dσ y , x ∈ D, (3.2)
with K(x, y) = n y • x -y |x -y| 2 . (3.3)
The density, µ, satisfies the boundary integral equation,

- 1 2 µ(y) + 1 2π ∂D K(y, y )µ(y )dσ y = f (y), y ∈ ∂D. (3.4)
In what follows, we assume that we have solved (3.4) using PTR N .

To evaluate (3.2) when x is close to the boundary, we set

x = y * - ε |κ * | n * y , (3.5)
where y * is the closest point to x on the boundary, n * y is the unit, outward normal at y * , κ * is the signed curvature at y * , and ε > 0 is a small parameter. Fig. 3 gives a sketch of these quantities. Substituting (3.5) into (3.3) yields We have written K in (3.6) to reveal its inherent dependence on the stretched variable, y = y * + εY/|κ * |.

K y * - ε |κ * | n * y , y = |κ * | ε n y • |κ * |(y * -y)/ε -n y • n * y |κ * (y * -y)/ε| 2 -2n * y • |κ * |(y * -y)/ε + 1 . ( 3

Matched asymptotic expansion of the kernel

We determine the matched asymptotic expansion of (3.6) [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory[END_REF]. Consider first the outer expansion in which y * and y are held fixed and ε → 0 + , so that |Y| → ∞. To leading order, we find that

K out ∼ - |κ * | ε n y • Y |Y| 2 = n y • (y * -y) |y * -y| 2 . (3.7)
The error of (3.7) is O(ε). Since this outer expansion is the kernel in (3.4), we find that

1 2π ∂D K out (y * -y)µ(y)dσ y = f (y * ) + 1 2 µ(y * ). (3.8)
The inner expansion is (3.6) written in terms of the stretched variable, Y. We seek an explicit parameterization of this inner expansion using y(t) = (y 1 (t), y 2 (t)), with t ∈ [0, 2π]. It follows that dσ y = |y (t)|dt, the unit tangent is τ y (t) = (y 1 (t), y 2 (t))/|y (t)|, the outward unit normal is n y (t) = (y 2 (t), -y 1 (t))/|y (t)|, and the signed curvature is κ(t) = (y 1 (t)y 2 (t) -y 1 (t)y 2 (t))/|y (t)| 3 . Let y * = y(t * ) and κ * = κ(t * ) with t * ∈ [0, 2π]. We introduce the stretched parameter, t = t * + εT , and find by expanding about ε = 0 that

y(t * + εT ) = y(t * ) + εT |y (t * )|τ y (t * ) - 1 2 ε 2 T 2 κ * |y (t * )| 2 n y (t * ) -(τ y (t * ) • y (t * ))τ y (t * ) + O(ε 3 ). (3.9) It follows that n y (t * + εT ) • |κ * |(y(t * ) -y(t * + εT )) = - 1 2 ε 2 T 2 γ * + O(ε 3 ), (3.10) n y (t * ) • |κ * |(y(t * ) -y(t * + εT )) = 1 2 ε 2 T 2 γ * + O(ε 3 ), ( 3.11
)

n y (t * + εT ) • n y (t * ) = 1 - 1 2 ε 2 T 2 |γ * | 2 + O(ε 3 ), ( 3.12)
and

|κ(t * )[y(t * ) -y(t * + εT )]| 2 = ε 2 T 2 |γ * | + O(ε 3 ), ( 3
K in (T ; ε) = |κ(t * )| -ε -1 2 ε 2 T 2 γ * + O(ε 3 ) ε 2 T 2 |γ * | + ε 2 + O(ε 3 ) . (3.14)
Next, we substitute ε 2 T 2 ∼ 2-2 cos(t-t * ) into (3.14) and determine that the leading order asymptotic behavior of K in is given by

K in (t -t * ; ε) ∼ |κ(t * )| -(γ * + ε) + γ * cos(t -t * ) (2|γ * | + ε 2 ) -2|γ * | cos(t -t * ) . (3.15)
The error of (3.15) is at most O(ε).

To form the leading order matched asymptotic expansion, we establish asymptotic matching in the overlap region of the outer and inner expansions. We first evaluate (3.7) in the limit as y → y * and find that

K out → - κ * 2 , y → y * . (3.16)
Next, we evaluate (3.15) in the limit as ε → 0 + and find that

K in → - κ * 2 , ε → 0 + . (3.17)
Thus, the overlapping value is -κ * /2. It follows that the matched asymptotic expansion for the kernel of the double-layer potential is given by

K = K out + K in + κ * 2 + O(ε), ε → 0 + , (3.18)
with K out given in (3.7) and K in given in (3.15). The error of this matched asymptotic expansion has O(ε) error because the error of K out given by (3.7) is O(ε), and K in given by (3.15) is at most O(ε).

For example, we plot K, K in , and K out in Fig. 4 as a function of t -t * with t * = π and ε = 0.1 for the boundary curve r(t) = 1 + 0.3 cos 5t. The right plot shows the L ∞ -error made by (3.18) as a function of ε. The solid curve is the linear fit through these data and has slope 1.2034 consistent with the O(ε) error.

Fourier coefficients of K in

The inner expansion given by (3.15) accurately captures the sharp peak of the kernel at t = t * in the limit as ε → 0 + as shown in Fig. 4. To avoid using PTR N to integrate over this sharp peak, we seek the Fourier coefficients,

Kin [n] = 1 2π 2π 0 K in (t; ε)e -int dt, (3.19)
so that we may use an approximation similar to that given in (2.16). To do so, we rewrite (3.15) as

K in (t -t * ; ε) = - |κ(t * )| C 0 1 2 A 0 -A 1 cos(t -t * ) 1 + C 1 cos(t -t * ) , (3.20)
with [START_REF] Kress | Linear Integral Equations[END_REF] gives K in as a rational function of trigonometric polynomials which have been studied by Geer [START_REF] Geer | Rational trigonometric approximations using fourier series partial sums[END_REF] in the context of constructing Fourier-Padé approximations. Since |C 1 | < 1, we have where ρ = 1 -C 2 1 -1 /C 1 . We find that we can improve on our approximation by considering the specific case in which the boundary is a circle of radius a. For that case, K out = -κ * /2 which cancels with the asymptotic matching term in (3.18). If we set

A 0 = 2(γ * + ε), A 1 = γ * , C 0 = 2|γ * | + ε 2 , and C 1 = -2|γ * |/C 0 . Equation (3.
Kin [n] =    1+ρ 2 1-ρ 2 A0 2 + A 1 ρ , n = 0 1+ρ 2 1-ρ 2 A0ρ |n| 4 + A 1 (ρ |n|-1 + ρ |n|+1 ) , n = 0 , (3.21)
A 0 = 2(γ * + ε -ε|γ * |), (3.22a)
A 1 = γ * -ε|γ * |, (3.22b)
C 0 = 2(|γ * | -εγ * ) + ε 2 , (3.22c
)

C 1 = -2(|γ * | -εγ * )/C 0 , (3.22d)
instead of the coefficients defined above, we find that (3.20) gives the exact evaluation of the kernel at r = a(1 -ε). For this reason, we use (3.22) in (3.20) and (3.21) in practice. These coefficients just include the O(ε 3 T 2 ) terms in the asymptotic expansion of K in for a general boundary.

To compute the contribution by K in to the double-layer potential, we use the approximation

1 2π 2π 0 K in (t -t * ; ε)µ(t)|y (t)|dt ≈ N/2-1 n=-N/2 Kin * [n]μ y [n]e -int * , (3.23

Evaluating the double-layer potential

The new method developed here for close evaluation of the double-layer potential uses (3.18). For convenience, let us introduce the residual kernel,

K = K -K out -K in - κ * 2 . ( 3
u y * - ε |κ * | n * y ≈ 1 2π 2π 0 K(t -t * ; ε) + κ * 2 µ(t)|y (t)|dt+f (y(t * ))+ 1 2 µ(t * )+ N/2-1 n=-N/2
Kin * [n]μ y [n]e -int * .

(3.27) Applying PTR N with t j = 2πj/N to the remaining integral in (3.27), we arrive at

u y * - ε |κ * | n * y ≈ 1 N N j=1 K(t j -t * ; ε) + κ * 2 µ(t j )|y (t j )|+f (y(t * ))+ 1 2 µ(t * )+ N/2-1 n=-N/2
Kin * [n]μ y [n]e -int * .

(3.28) Equation (3.28) gives our method for computing the double-layer potential for close evaluation points. It avoids aliasing incurred by the sharp peak of K in by using (3.23). Integration of K out is replaced by f (y * ) + 1 2 µ(y * ), which comes from evaluating boundary integral equation (3.4) at y * . PTR N is now used only to integrate the term with the kernel, K + κ * /2. This term is important for taking into account additional, non-local contributions to the double-layer potential, which may be significant.

Numerical results

We present results of this method for evaluating the double-layer potential by computing the harmonic function, u(x) = -1 2π log |x -x 0 | with x 0 = (1.85, 1.65). We compute the solution interior to the boundary curve, r(t) = 1.55 + 0.4 cos 5t. The Dirichlet data in (3.1b) is determined by evaluating the harmonic function on the boundary. Boundary integral equation (3.4) is solved using PTR N and we use the resulting density, µ(t j ) with t j = 2πj/N for j = 1, • • • , N in the double-layer potential. We evaluate the double-layer potential using two methods: (1) PTR N and (2) asymptotic PTR N , the new method given in (3.28). We present the solution on a body-fitted grid in which evaluation points are found by moving along the normal into the domain from boundary grid points. The solution is evaluated on a grid of 200 equispaced points along each normal starting at the boundary until we reach a distance 1/κ max from the boundary, where κ max = max 0≤t * <2π |κ(t * )|. This grid captures the boundary layer, but does not coincide exactly with it since the boundary layer depends on the local curvature. For regions of high curvature, this body-fitted grid extends beyond the boundary layer.

In Fig. 5 we show the errors (log 10 -scale) in computing the double-layer potential using PTR 128 and asymptotic PTR 128 . These results show that asymptotic PTR N produces errors that are several orders of magnitude smaller than those of the PTR N . To give an indication of this improvement, the L ∞ error is 8.03 for PTR 128 and 1.85 × 10 -4 for asymptotic PTR 128 . To examine this error in more detail, we show in Fig. 6 a plot of the error in computing the double-layer potential evaluated at the points indicated in Fig. 5 (t * = 0, t * = π/2, and t * = π) as a function of ε. These three cases are plotted over different ranges of ε corresponding to 0 < ε < κ(t * )/κ max . We observe that asymptotic PTR N does significantly better than PTR N for small ε as expected. It reduces the O(1) error by at least four orders of magnitude. We find similar results over all values of t * .

We can further improve the new method using the identity for the double-layer potential [START_REF] Kress | Linear Integral Equations[END_REF] (see (5.1)) to rewrite (3.2) as follows:

u(x) = 1 2π ∂D K(x, y)(µ(y) -µ(y * ))dσ y -µ(y * ), x ∈ ∂D. ( 3
with g denoting an analytic function satisfying ∂D g(y)dσ y = 0. (

We seek v as the single-layer potential [START_REF] Kress | Linear Integral Equations[END_REF],

v(x) = 1 2π ∂D S(x, y)ϕ(y)dσ y , x ∈ R 2 \ D, (4.3)
with S(x, y) = -log |x -y|.

The density, ϕ(y), satisfies the boundary integral equation,

- 1 2 ϕ(y) + 1 2π ∂D ∂S(y, y ) ∂n y ϕ(y )dσ y = g(y), y ∈ ∂D. (4.5)
To study the close evaluation of (4.3), we now set

x = y * + ε |κ * | n * y . (4.6)
Substituting (4.6) into (4.4), we obtain

S y * + ε |κ * | n * y , y = -log ε + log |κ * | - 1 2 log |κ * (y * -y)/ε| 2 + 2n * y • |κ * |(y * -y)/ε + 1 . (4.7)
Just as we have done for K in (3.6), we have written (4.7) to show the underlying dependence on the stretched variable, y = y * + εY/|κ * |.

The outer expansion of (4.7) is S out ∼ -log |y * -y|. This outer expansion is singular. In contrast to the double-layer potential, this outer expansion does not correspond to the kernel boundary integral equation (4.5). One could use a high order quadrature rule that explicitly takes into account this singularity [START_REF] Avram | Quadrature methods for periodic singular and weakly singular Fredholm integral equations[END_REF][START_REF] Kress | Boundary integral equations in time-harmonic acoustic scattering[END_REF]. However, we choose to not use one here because we find in the numerical examples below that our method significantly reduces the dominant error. To compute the inner expansion, S in , we introduce the stretched parameter, t = t * + εT into the same parameterization of the boundary used for the double-layer potential. Making use of (3.10) and (3.12), we find by expanding as

ε → 0 + that S in (T ; ε) = log |κ * | - 1 2 log ε 2 T 2 |γ * | + ε 2 + O(ε 3 ) . (4.8)
Substituting ε 2 T 2 ∼ 2 -2 cos(t -t * ), we find that to leading order,

S in (t -t * ; ε) ∼ log |κ * | - 1 2 log (2|γ * | + ε 2 ) -2|γ * | cos(t -t * ) . (4.9)

Fourier coefficients of S in

Using the modified coefficients introduced in (3.22), we write (4.9) as

S in (t -t * ; ε) ∼ log |κ * | - 1 2 log C 0 - 1 2 log[1 + C 1 cos(t -t * )]. (4.10)
We now seek to compute

Ŝin [n] = δ n,0 log |κ * | - 1 2 log C 0 - 1 2π 2π 0 1 2 log[1 + C 1 cos(t -t * )]e -int dt, (4.11)
with δ n,0 denoting the Kronecker delta. To compute the integral in (4.11), we start with

d dt log[1 + C 1 cos(t -t * )] = C 1 sin(t -t * ) 1 + C 1 cos(t -t * ) . (4.12)
The right-hand side of (4.12) is another example of a rational trigonometric function studied by Geer [START_REF] Geer | Rational trigonometric approximations using fourier series partial sums[END_REF]. It can be readily shown that

1 2π 2π 0 sin(t -t * ) 1 + C 1 cos(t -t * ) e int dt = sgn(n) i 2 1 + ρ 2 1 -ρ 2 ρ |n|-1 -ρ |n|+1 , (4.13)
where ρ = 1 -C 2 1 /C 1 . It follows from term-by-term integration of the Fourier series with these coefficients that

Ŝin [n] =        log |κ * | - 1 2 log C 0 - C 1 2 1 + ρ 2 1 -ρ 2 1 ρ -ρ log(1 -ρ) n = 0, C 1 4|n| 1 + ρ 2 1 -ρ 2 ρ |n|-1 -ρ |n|+1 n = 0. (4.14)

Evaluating the single-layer potential

Given the inner expansion computed above, our method for evaluating the single-layer potential is to compute an approximation of

v y * + ε |κ * | n * y = 1 2π 2π 0 S(t -t * )ϕ(t)|y (t)|dt + 1 2π 2π 0 S in (t -t * ; ε)ϕ(t)|y (t)|dt, (4.15)
with S = S -S in . We use PTR N to evaluate the first integral with kernel, S, and a truncated convolution sum to evaluate the second integral with kernel, S in , In Fig. 7 we show the absolute error (log 10 -scale) in computing the single-layer potential using the PTR 128 and using asymptotic PTR 128 . The single-layer potential kernel is not as sharply peaked as the double-layer potential kernel, so the error in evaluating the single-layer potential is less than the error when evaluating the double-layer potential. Even so, we still observe a boundary layer of thickness O(1/N ) in which the error is O(1) due to aliasing when using PTR N . Asymptotic PTR N effectively reduces the error in the boundary layer. To give an indication of this improvement, the L ∞ error is 0.113 for PTR 128 and 5.39 × 10 -5 for asymptotic PTR 128 . In Fig. 8, we plot the error in computing the single-layer potential evaluated at the points indicated in Fig. 7 (t * = 0, t * = π/2 and t * = π) as a function of ε with 0 < ε < κ(t * )/κ max . These plots show that the asymptotic method reduces the error by at 3 orders of magnitude for small ε. We find similar results over all values of t * . For the case in which t * = π, the error of asymptotic PTR N becomes larger than that for PTR N for 0.36 < ε < 0.52. For this particular boundary curve, κ max is attained at t * = π. Hence, the body-fitted grid at t * = π plots the single-layer potential over 0 < ε < 1. For this range of ε, we consider points outside the boundary layer where PTR N is competitive with, and may even become more accurate than asymptotic PTR N . In fact, that is what is observed in Fig. 8 for t * = π.

v y * + ε |κ * | n * y ≈ 1 N N j=1

General implementation

In the results presented above, we are using a body-fitted grid, in which evaluation points are found by moving along the normal into the domain from boundary integration points. Often the solution is needed at more generally defined points. Asymptotic PTR N tacitly assumes in (3.5) for interior problems, or (4.6) for exterior problems, that y * , is the unique minimum distance from the boundary to the evaluation point, x. In the examples, the closest point on the boundary, y * coincides with a PTR N grid point from which we extended in the normal direction. We discuss here the more general problem. Suppose we have an evaluation point in the domain, x. Then the first problem to address is whether x is, in fact, close enough to the boundary to warrant special attention and requires use of asymptotic PTR N . To solve this problem, we make use of the identity for the double-layer potential [START_REF] Kress | Linear Integral Equations[END_REF],

1 2π ∂D n y • (x -y) |x -y| 2 dσ y =      -1 x ∈ D -1 2 x ∈ ∂D 0 x ∈ R 2 \ D .
(5.1) Evaluating (5.1) using PTR N suffers from the same aliasing problem that the more general layer potential evaluations do. Thus, we use it to determine if x lies within the boundary layer. To do this, we set a user-defined threshold for the error. If the error in evaluating (5.1) using PTR N is less than this threshold value, we keep the result computed using PTR N . Otherwise, we use the appropriate asymptotic approximation.

To use these asymptotic approximations, we must determine the parameter, t * , where t * = min 0≤t<2π |x -y(t)| 2 . For a general boundary, this problem may not have a unique solution. In practice, we find a unique minimizer for evaluation points that are identified to be in the boundary layer using PTR N evaluation of (5.1). Once t * is determined, all other quantities required for the asymptotic approximations follow. Finally, we evaluate the solution of the boundary value problem at hand at any point x using either PTR N or asymptotic PTR N .

We present results of this generalized method for evaluation of the double-layer potential and single-layer potential for the same problems presented in Section 3.4 and 4.3, respectively. We use a threshold of 1 × 10 -8 , as described above, to determine when the evaluation point is inside the boundary layer and asymptotic PTR N method is to be used. In Fig. 9 we show the error in computing the double-layer potential using PTR 256 and asymptotic PTR 256 when solving on a Cartesian grid with meshsize h = 0.005 within the boundary curve. Similarly, we present the evaluation single-layer potential in Fig. 10. Here, we are computing on a Cartesian grid with meshsize h = 0.005 exterior to the boundary curve. These results show, similar to the results while considering a body-fitted grid, that the error made by asymptotic PTR N is several orders of magnitude smaller than those made by PTR N . However, there are more variations in these errors because y * does not always coincide with a quadrature point. We choose to use 256 quadrature points here as this is what is actually needed to solve the boundary integral equations for the densities such that µ(t)|y (t)| and ϕ(t)|y (t)| are sufficiently resolved. We were able to use less points for the body-fitted grid as this restriction is not as strict when y * is a quadrature point.

Conclusions

We have presented a new method to address the close evaluation problem. When solving boundary value problems using boundary integrals equation methods, the solution is evaluated at desired points within the domain by numerically evaluating layer potentials. Using the same quadrature rule that is used to solve the integral equation for this evaluation achieves high order accuracy everywhere in the domain, except close to the boundary where an O(1) error is incurred. The new method developed here takes advantage of the knowledge of the sharply peaked kernel of layer potentials close to the boundary to reduce this error by several orders of magnitude. We have used asymptotic methods to analyze the kernels of the double-and single-layer potentials to relieve the numerical having to integrate over this sharply peaked kernel. The resulting method is straightforward to implement. We have presented results for both the interior Dirichlet problem and exterior Neumann problem for Laplace's equation and show a reduction in error of four to five orders of magnitude in the solution evaluation close to the boundary. Furthermore, we have presented how to generalize this method to solve within the whole domain, including how to determine when to use the new asymptotic method.

This asymptotic method has been recently applied to acoustic scattering problems by sound-soft obstacles [START_REF] Carvalho | Local analysis of near fields in acoustic scattering[END_REF]. For those problems, the sharp peaks in the kernels for the single-and double-layer potentials, both which are needed, have the same character as those for Laplace's equation discussed here. Hence, only small modifications are needed to apply these methods to wave propagation problems. We are currently extending this asymptotic method to three-dimensional problems. Furthermore, we are working on different applications which include extending the method presented here to the Stokes equations and studying surface plasmons.

Figure 2 :

 2 Figure 2: Contour plot of log 10 |E N (r, t)| where E N is given by (2.15) with a = 1 and N = 128.

  ), and |E N (r, t * )| → 1/2 as r → a. These results show that E N has a boundary layer of thickness O(1/N ) where it exponentially increases to values that are O(1). Fig.2shows a plot of (2.15) over the entire circular disk and a close-up near the boundary. These plots show the boundary layer about r = a where the error attains O(1) values.

  Consider a simply connected, open set denoted by D ⊂ R 2 , with analytic boundary ∂D. Let D = D ∪ ∂D. The function u ∈ C 2 (D) ∩ C 1 (D) satisfies ∆u = 0 in D, (3.1a)

Figure 3 :

 3 Figure 3: Sketch of the quantities introduced in (3.5) to study evaluation points close to the boundary.

  .13) with γ * = sgn[κ * ]|κ * y (t * )| 2 . Here, sgn[x] = x/|x| for x = 0 and sgn[x] = 0 for x = 0. Substituting (3.10) -(3.13) into (3.6), we find that

Figure 4 :

 4 Figure 4: [Left] Plot of the kernel, K, given in (3.6) (solid curve) and the leading order behavior of its inner expansion, K in given in (3.15) (dashed curve) and its outer expansion, K out given in (3.7) (dotted curve) as a function of t-t * with t * = π and ε = 0.1 for the boundary curve, r(t) = 1+0.3 cos 5t. [Right] L ∞ -error of the matched asymptotic expansion given in (3.18) evaluated at t * = π for ε = 0.0001, 0.001, 0.01, and 0.1. These computed errors are plotted as circles. The solid curve gives the result of fitting this data to the function, Cε p . This fit produced p = 1.2034 indicating the O(ε) error of the matched asymptotic expansion.

Figure 5 :Figure 6 :

 56 Figure 5: [Left] Plot of absolute error (log 10 -scale) in computing the double-layer potential using PTR 128 for the boundary r(t) = 1.55 + 0.4 cos 5t for the Dirichlet data, f (y) = 1 2π log |y -x 0 | with x 0 = (1.85, 1.65). [Right] Plot of absolute error (log 10 -scale) in computing the double-layer potential using the asymptotic PTR 128 given in (3.28) for the same problem. The "×" symbols on the boundary indicates the points corresponding to t * = 0, t * = π/2, and t * = π.

4. 3

 3 Numerical examplesWe present results for the evaluation of the single-layer potential by computing the harmonic function, v(x) = (x -x 0 )/|x -x 0 | 2 , with x 0 = (0.1, 0.4). We compute the solution exterior to the the boundary curve, r(t) = 1.55 + 0.4 cos 5t. The Neumann data in (4.1b) is determined by computing the normal derivative of the harmonic function on the boundary. Boundary integral equation (4.5) is solved using PTR N and we use the resulting density ϕ(t j ) with t j = 2πj/N for j = 1, • • • , N in the single-layer potential. We evaluate the single-layer potential using two methods: (1) PTR N and (2) asymptotic PTR N , the new method given in (4.16). We modify the body-fitted grid described above for the evaluation of the double-layer potential to evaluate exterior points. The solution is evaluated on a grid of 200 equispaced points along each normal starting at the boundary until we reach a distance 1/κ max from the boundary.

Figure 7 :

 7 Figure 7: [Left] Plot of the absolute error (log 10 -scale) in computing the single-layer potential using PTR 128 for the boundary r(t) = 1.55 + 0.4 cos 5t for the Neumann data, g(y) = ∂v ∂n with v(x) = (x -x 0 )/|x -x 0 | 2 , x 0 = (0.1, 0.4). [Right] Plot of the absolute error (log 10 -scale) in computing the single-layer potential using asymptotic PTR 128 given in (4.16) for the same problem. The "×" symbols on the boundary indicates the points corresponding to t * = 0, t * = π/2 and t * = π.

Figure 8 :

 8 Figure8: Plot of the absolute error as a function of ε made in evaluating the single-layer potential for different t * using PTR 128 (solid curve) and using asymptotic PTR 128 (dashed curve).

Figure 9 :

 9 Figure 9: [Left] Plot of the absolute error (log 10 -scale) in computing the double-layer potential using PTR 256 for the boundary r(t) = 1.55 + 0.4 cos 5t for the Dirichlet data, f (y) = 1 2π log |y -x 0 | with x 0 = (1.85, 1.65). We evaluate the solution inside the domain on a regular grid. [Right] Plot of the absolute error (log 10 -scale) in computing the double-layer potential using asymptotic PTR 256 given in (3.28) for the same problem. The asymptotic method is used in a boundary layer determined by a threshold on the error from evaluating (5.1).

Figure 10 :

 10 Figure 10: [Left] Plot of the absolute error (log 10 -scale) in computing the single-layer potential using PTR 256 for the boundary r(t) = 1.55 + 0.4 cos 5t for the Neumann data, g(y) = ∂v ∂n with v(x) = (x -x 0 )/|x -x 0 | 2 , x 0 = (0.1, 0.4). [Right] Plot of log 10 of the absolute error (log 10 -scale) in computing the single-layer potential using the asymptotic PTR 256 given in (4.16) for the same problem. The asymptotic method is used in a boundary layer determined by a threshold on the error from evaluating (5.1).

  .25) K = O(ε) and more importantly, it does not have a sharp peak about t = t * . We rewrite the double-

	layer potential as
	u y * -	ε |κ * |	n * y	=	1 2π	0	2π	K(t -t * ; ε)µ(t)|y (t)|dt +	1 2π	0	2π	K out (t -t * ; ε)µ(t)|y (t)|dt
								+	1 2π	0	2π	K in (t -t * ; ε)µ(t)|y (t)|dt +	κ * 4π	0	2π	µ(t)|y (t)|dt. (3.26)
	Substituting (3.8) and (3.23) into (3.26), we obtain

  .29) In(3.29), the integrand is now smoother as it vanishes at the point y = y * , and the error using PTR N drastically decreases. Applying asymptotic PTR N to (3.29), we obtain For the example problem discussed above, the L ∞ error is 7.13 × 10 -5 for PTR 128 applied to (3.29), and 2.54 × 10 -5 for asymptotic PTR 128 given by(3.30). This additional improvement (3.29) is only valid for the double-layer potential and does not generalize.

	u y * -	ε |κ * |	n * y	≈	1 N	N j=1	K(t j -t * ; ε) +	κ * 2	n=-N/2 (µ(t j )-µ(t N/2-1	Kin * [n]μ * y [n]e -int *
											(3.30)
	with						μ * y [n] =	1 2π	0	2π	(µ(t) -µ(t * ))|y (t)|e -int dt.	(3.31)

* ))|y (t j )|+f (y(t * ))+

4

  Single-layer potential for the exterior Neumann problem

	We now consider the exterior Neumann problem,
	∆v = 0 in R 2 \ D,	(4.1a)
	∂v ∂n	= g on ∂D,

Acknowledgments

The authors thank François Blanchette and Boaz Ilan for their thoughtful discussions leading up to this work. S. Khatri was supported in part by the National Science Foundation (PHY-1505061). A. D. Kim acknowledges support by Air Force Office of Scientific Research (FA9550-17-1-0238) and the National Science Foundation.


Mayela Authors

Barbara Cepeda-Vázquez
email: mayela.cepedavazquez@agroparistech.fr

Nicolas Rega

Valérie Descharles
email: nicolas.descharles@agroparistech.fr

Camel

How ingredients influence furan and aroma generation in sponge cake

Keywords: formulation, furanic compounds, glucose, egg white, heat-processed foods, quantification, headspace trap, multivariate analysis Highlights

A wide range of compounds can be formed during thermal processing of food, some of which are relevant for aroma (e.g. furfural), while others are of great health-concern (e.g. furan). This paper presents the study of formulation as affecting the simultaneous generation of furan and furfural, along with other aroma quality markers, in sponge cake by means of headspace trap / GC-MS. Ingredients were screened according to their category (fat, salt, sugar, egg-based). Glucose-containing formulation resulted in the highest content of furan and furfural (12.5±0.5 ng g -1 and 9.2±0.2 μg g -1 , dry basis, respectively), while their lowest amount was found in the egg-white recipe (3.1±0.1 ng g -1 for furan and 0.287±0.078 µg g -1 for furfural, dry basis). The latter also related negatively to all studied compounds. This work will be useful for developing novel strategies to deliver safe foods with appealing organoleptic attributes.

INTRODUCTION

Heat-treated foods, such as baked goods, are widely consumed among population, and their investigation is especially important from a food quality and safety point of view.

During thermal processing food constituents may induce chemical changes through reaction pathways (such as caramelization, Maillard reaction and lipid oxidation), leading to the formation of a great number of compounds, some of which are of health concern, whereas others play a role in aroma and color. A current challenge is thus, the development of safe yet appealing food for which food reactivity must be taken into account.

In recent years, much attention has been paid to furan, a heat-generated compound of health-concern due to its possible carcinogenicity, occurrence and consumption by population (International Agency for Research on Cancer, 1995; Joint FAO/WHO Expert Committee on Food Additives, 2011). It is formed through reaction pathways that are common to other compounds having a role in sensory quality [START_REF] Belitz | Food Chemistry[END_REF][START_REF] Crews | A review of the occurrence, formation and analysis of furan in heat-processed foods[END_REF][START_REF] Kroh | Caramelisation in food and beverages[END_REF][START_REF] Martins | A review of Maillard reaction in food and implications to kinetic modelling[END_REF][START_REF] Perez Locas | Origin and mechanistic pathways of formation of the parent furan -a food toxicant[END_REF][START_REF] Stadler | Heat-generated toxicants in foods: acrylamide, MCPD esters and furan[END_REF], such as furfural, a naturally occurring odorant in many food items [START_REF] Burdock | Fenaroli's handbook of flavor ingredients[END_REF][START_REF] Stofberg | Comsumption ratio and food predominance of flavoring materials[END_REF]. Getting insight into their simultaneous formation and their concentration in food is therefore important. Furan has been extensively quantified in a broad range of products, including baked goods [START_REF] Crews | A review of the occurrence, formation and analysis of furan in heat-processed foods[END_REF]; Joint FAO/WHO Expert Committee on Food Additives, 2011; [START_REF] Zoller | Furan in food: Headspace method and product survey[END_REF]. Furfural has also been widely investigated in food, including qualitative, semi-quantitative and ultimately quantitative studies in bakery products (Ait Ameur, [START_REF] Ait Ameur | The fate of furfurals and other volatile markers during the baking process of a model cookie[END_REF][START_REF] Maire | Lipid oxidation in baked products : Impact of formula and process on the generation of volatile compounds[END_REF][START_REF] Petisca | Assessment of hydroxymethylfurfural and furfural in commercial bakery products[END_REF][START_REF] Pico | Wheat bread aroma compounds in crumb and crust : A review[END_REF][START_REF] Poinot | Influence of formulation and process on the aromatic profile and physical characteristics of bread[END_REF][START_REF] Pozo-Bayón | Influence of Eggs on the Aroma Composition of a Sponge Cake and on the Aroma Release in Model Studies on Flavored Sponge Cakes[END_REF][START_REF] Rega | On-line dynamic HS-SPME for monitoring endogenous aroma compounds released during the baking of a model cake[END_REF]. Up to date, however, furan and furfural simultaneous quantitative determination in food are scarce [START_REF] Chaichi | Furanic compounds and furfural in different coffee products by headspace liquid-phase micro-extraction followed by gas chromatography -mass spectrometry : survey and effect of brewing procedures[END_REF][START_REF] Pérez-Palacios | Impact of cooking and handling conditions on furanic compounds in breaded fish products[END_REF]Catarina Petisca, Pérez-Palacios, Pinho, & Ferreira, 2014b), and no studies are reported in baked goods.

Very recently, we optimized headspace trap (HS trap) extraction for furan and furfural simultaneous quantification in sponge cake by gas chromatography coupled to mass spectrometry (GC-MS), which proved to be sensitive, precise and competitive in comparison to other techniques [START_REF] Cepeda-Vázquez | Multivariate optimization of headspace trap for furan and furfural simultaneous determination in sponge cake[END_REF]. This method has been successfully applied to sponge cakes having different properties (fat content and humidity) and can thus be used for formulation studies in similar matrices.

Sponge cake is indeed an interesting bakery product to study the relationship between quality and reactivity, since it is made-up with several ingredients (i.e. sugar, salt, egg, fat and flour), leading to a complex chemical composition, and contains a number of compounds that may also have a role in aroma quality [START_REF] Maire | Lipid oxidation in baked products : Impact of formula and process on the generation of volatile compounds[END_REF][START_REF] Pozo-Bayón | Influence of Eggs on the Aroma Composition of a Sponge Cake and on the Aroma Release in Model Studies on Flavored Sponge Cakes[END_REF][START_REF] Rega | On-line dynamic HS-SPME for monitoring endogenous aroma compounds released during the baking of a model cake[END_REF]. Moreover, studies in food are all-important because of the diverse nature of ingredients, their possible synergistic effect and impact on overall quality. This paper presents the study of ingredients in sponge cake as influencing the generation of furan, furfural and other relevant aroma compounds (such as aldehydes and pyrazines) by means of HS trap / GC-MS.

Chemicals and standards

Furan (purity: 99%), furfural (purity: 99%), 2-methylbutanal (purity: 95%), 2,5dimethylpyrazine (purity: 99.9%), benzaldehyde (purity: 99%), hexanal (purity: 98%), 2pentylfuran (purity: 97%), 1-hydroxy-2-propanone (purity: 90%) and d4-furan (purity: 98%) were supplied by Sigma-Aldrich. 3-Methylbutanal (purity: 98%) was obtained from Merck, 2,3,5-trimethylpyrazine was supplied by Acros (purity: 99%) and d4-furfural (purity: 99.7%) was obtained from CIL Cluzeau, CDN Isotopes (Pointe-Claire, Canada). Methanol (purity: 99.9%) was purchased from Carlo Erba (Val-de-Reuil, France). Ultrapure water (> 5 MΩ.cm) was produced using a Millipore Elix 3 System from Millipore SAS (Molsheim, France). Heptane and tetrahydrofuran were provided by Carlo Erba (Val de Reuil, France) and standard (±)-α-tocopherol by Sigma (Steinhein, Germany).

Ingredients

Flour was kindly provided by Grands Moulins de Paris (Ivry-sur-Seine, France). Caster sucrose was obtained from Tereos (Lille, France) and glucose (dextrose monohydrate) was kindly provided by Roquette Frères (Lestrem, France). Non hydrogenated palm oil was obtained from Palva Celys (Rezé, France) and sunflower oil by Lesieur (Asnières-sur-Seine, France). Pasteurized whole egg and egg white were supplied by Agrodoubs (Flagey, France) and Ovoteam (Massy, France). NaCl was supplied by Salinor (Clichy, France).

Sucrose, glucose and salt were stored at room temperature, whole egg and sunflower oil were kept at 4°C and flour, egg white and palm oil were stored at -20°C. Flour and egg white were thawed at 4°C and palm oil in a water bath at 50°C prior to use.

A batch of purified sunflower oil was also obtained by removing tocopherols (naturally occurring antioxidants). Purification was carried out in three successive vacuum filtrations (150 mbar and twice 120 mbar) of 200 g of sunflower oil through 70 g of dry alumina (adapted from [START_REF] Yoshida | Participation of Free Fatty Acids in the Oxidation of Purified Soybean Oil During Microwave Heating[END_REF]). This purified sunflower oil was put into amber glass bottles filled with nitrogen and stored at 4°C until use.

Absence of tocopherols in purified sunflower oil was confirmed by high performance liquid chromatography (HPLC), for which 1 g of oil was diluted in heptane for a final concentration of 0.1 g mL -1 . Sample was vortexed for 60 s at maximum speed and put into an ultrasound bath for 30 s for proper homogenization, followed by filtration through a 0.22 nm nylon filter membrane (Interchim, Montluçon, France). Analyses were carried out under isocratic conditions in an Alliance Waters 2695 system equipped with a photodiode array detector (298 nm) and a Lichrospher 100 DIOL column (5 μm, 250 x 4 mm, AIT France, Houilles, France). Mobile phase was a 96.15/3.85 (v/v) heptane/tetrahydrofuran mixture (1 mL min -1 ) and standard tocopherol concentration was 100 mg L -1 (adapted from [START_REF] Afnor | Corps gras d'origines animale et végétale -Détermination des teneurs en tocophérols et en tocotriénols par chromatographie en phase liquide à haute performance[END_REF]).

Sponge cake preparation

Formulations for the screening of ingredients

Reference sponge cake (f1) was formulated according to [START_REF] Cepeda-Vázquez | Multivariate optimization of headspace trap for furan and furfural simultaneous determination in sponge cake[END_REF].

Five other formulations (f2-f6) were initially prepared in order to perform the screening of ingredients in relation to their category (fat, salt, sugar, egg-based) and their ability to act as precursors of different reaction pathways, namely lipid oxidation, caramelization, Maillard reaction, Strecker's degradation and other chemical interactions (see Table 1).

Concerning the fat category, three types of formulas were designed based on polyunsaturated fatty acids (PUFAs) content (mainly linolenic and linoleic acids, well-known furan precursors) [START_REF] Crews | A review of the occurrence, formation and analysis of furan in heat-processed foods[END_REF][START_REF] Stadler | Heat-generated toxicants in foods: acrylamide, MCPD esters and furan[END_REF], and thus their likeliness towards enhancing lipid oxidation. Palm oil (f1) is mainly composed of saturated and monounsaturated fat with around 9% of linoleic acid and a minor fraction of linolenic acid. Sunflower oil (f3) is particularly rich in linoleic acid (> 65%) (U.S. Department of Agriculture (USDA), 2016). Since these two oils naturally contain tocopherols (i.e. endogenous antioxidants), an additional formulation (f7) was prepared with purified sunflower oil to further test the lipid oxidation pathway in the absence of antioxidants.

Batter preparation

Batter of reference sponge cake (500 g) was prepared according to Cepeda-Vázquez et al.

(2017) by mixing and beating whole eggs, sucrose and salt using a kitchen appliance (Kenwood Chef, KM300, UK) for 10 min at maximum speed. Flour was gently added during 1.5 min at minimum speed and the whole was further beaten for 0.5 min. Palm oil was gently added during 15 s at minimum speed and the batter beaten for 1 min at the same speed. Additional sponge cakes were prepared by replacing or removing whole egg, sucrose, salt and palm oil from reference formulation according to Table 1. Formulation containing egg white was further battered for 120 s for proper incorporation.

Baking trials

For each formulation, aliquots of 20 g of batter were poured into 21 aluminum molds (8.0 cm x 4.5 cm x 3.5 cm) and baked at 170°C for 25 min in an instrumented oven (Bongard, Wolfisheim, France) specifically designed to ensure thermal homogeneity [START_REF] Fehaili | An instrumented oven for the monitoring of thermal reactions during the baking of sponge cake[END_REF].

Immediately after baking, sponge cakes were put into hermetically sealed glass jars and frozen at -20°C. Twenty sponge cakes were used for composite sampling for chemical analyses, and one for color measurements.

Independent baking trials (n=3) were performed to check the reproducibility of sample preparation on the reference formula (f1) as detailed in section 2.8.

Composite sampling for chemical analysis

For each baking trial, sponge cakes were ground by means of a Grindomix GM200 knife mill, equipped with a stainless steel bowl and titan knifes (Retsch GmbH, Haan, Germany) using the following conditions: five cakes per grinding operation -3,000 rpm for 10 s, 6,000 rpm for 20 s and 3,000 rpm for up to 50 s. All batches were then mixed together at 2,000 rpm for 5 s and the resulting composite sample was stored at -20°C until analysis. The same procedure was adopted for all formulas.

Physical and physicochemical properties

Sponge cake formulations (f1 to f6) were characterized by measuring the properties described below.

Dry matter determination

Moisture content was determined gravimetrically on 1 g of composite samples (n=3) after oven drying at 105°C (Memmert, Schwabach, Germany) to nearest milligram (AE 200,d: 0.1mg, Mettler, Switzerland) (adapted from [START_REF] Fehaili | An instrumented oven for the monitoring of thermal reactions during the baking of sponge cake[END_REF]).

pH measurement

Aliquots of 1 g of ground samples were mixed with 25 mL of ultrapure water (n=3), stirred for 1 min and further stirred during pH measurement of the resulting mixture (PH3, Hach Lange, Barcelona, Spain) (adapted from [START_REF] Ait Ameur | The fate of furfurals and other volatile markers during the baking process of a model cookie[END_REF]).

Color determination

Color measurements were performed on the upper face of a whole cake at 4 points using a 6834 spectro-guide sphere gloss colorimeter (BYK-Gardner, Germany) and expressed in the CIELab color space (adapted from [START_REF] Poinot | Influence of formulation and process on the aromatic profile and physical characteristics of bread[END_REF]).

Selection of quality markers

Furan and furfural were selected as main furanic quality markers for food safety and aroma, respectively. Other typical odorant compounds found in baked goods were selected for being representative of different chemical classes (furans, aldehydes, pyrazines, etc.) and reaction pathways (mainly caramelization, Maillard reaction, lipid oxidation and Strecker's degradation) (Table 2).

HS trap / GC-MS procedure

Furan and furfural were quantified by stable isotope dilution assays according to [START_REF] Cepeda-Vázquez | Multivariate optimization of headspace trap for furan and furfural simultaneous determination in sponge cake[END_REF], and concentrations were estimated. The remaining aroma compounds were studied using semi-quantitative analysis based on chromatographic peak areas and relative abundances were assessed.

Quantitative and semi-quantitative analyses (n=3) were performed simultaneously by means of a TurboMatrix Headspace Sampler HS 40 Trap (Perkin Elmer, Llantrisant, UK) equipped with an air monitoring trap containing a two-sorbent bed (Carbotrap and Carbosieve SIII, Perkin Elmer), and a Trace GC Ultra gas chromatography system coupled to an ISQ single quadrupole mass spectrometer (Thermo Scientific, Rodano, Italy).

Standard solutions for furan and furfural quantification

Stock and standard solutions were prepared in a room at 18°C, while HS samples were prepared in a separate room to avoid any cross-contamination from the surrounding air.

Stock solutions of d4-furan, d4-furfural, furan and furfural in methanol were prepared separately at a concentration of 2.5 g L -1 , put into hermetic amber glass bottles and stored at -20°C for up to one month. A standard solution of d4-furan at a concentration of 25 ng μL - 1 was prepared daily by diluting the stock solution in ultrapure water. A mixed standard solution containing 0.5 and 15 ng μL -1 of d4-furan and d4-furfural, respectively, was prepared daily by diluting the standard solution of d4-furan and the stock solution of d4furfural in ultrapure water and kept at 4°C until use. All glassware was previously baked out at 55°C for at least 24 h to avoid any contamination prior to analysis (Cepeda-Vázquez et al., 2017).

HS sample preparation

For vial preparation, sponge cake composite samples were kept in an ice bath to avoid volatile loss. Aliquots of ground sponge cake (0.588 g, dry basis) and ultrapure water (9.412 g) were weighed into 20 mL vials for a 16 water / sample amount ratio (dry basis) and 10 g of total amount (water + sample amount, dry basis) [START_REF] Cepeda-Vázquez | Multivariate optimization of headspace trap for furan and furfural simultaneous determination in sponge cake[END_REF]. All amounts were exactly measured in mass. A volume of 40 μL of the d4-furan + d4-furfural standard solution was added for a final set of 20 ng of d4-furan and 600 ng of d4-furfural per vial, respectively. Vials were capped with aluminum polytetrafluoroethylene coated silicone septa at each step of preparation and sealed immediately after labeled standards addition to avoid loss of volatiles. Vials were then vortexed for 5 s at maximum speed. All HS vials were previously baked out in a muffle furnace at 350°C for at least 1 h to avoid any HS glassware contamination before analysis [START_REF] Cepeda-Vázquez | Multivariate optimization of headspace trap for furan and furfural simultaneous determination in sponge cake[END_REF].

HS trap conditions

Conditions were set as follows, according to the multivariate optimization of HS trap extraction for furan and furfural simultaneous quantitative analysis [START_REF] Cepeda-Vázquez | Multivariate optimization of headspace trap for furan and furfural simultaneous determination in sponge cake[END_REF]: thermostatting temperature (65°C) and time (15 min), number of pressurization cycles (4) and dry purge time (0.9 min). For all assays, shaker mode was on during vial thermostatting (automatically varying through a frequency range). Needle temperature was set to 75°C and vials were pressurized to 40 psi during 0.5 min, followed by a decay time of 2 min per pressurization cycle. After HS extraction, the trap was heated from 28 to 250°C at a 40°C s -1 rate and a pressure of 23 psi, in order to desorb and release the analytes into the transfer line for 0.2 min. Trap hold time was set to 15 min between GC runs and a sample blank (empty vial) was run every three injections to confirm complete desorption. Vial and water blanks were run daily to check for the absence of contaminants.

GC-MS analysis

Analysis conditions were adapted from Cepeda-Vázquez et al. (2017). A 100% polyethylene glycol column (VF-WAXms) was used for separation (60 m x 0.25 mm x 0.5 μm; Agilent, Netherlands). Helium was used as carrier gas at a constant flow rate of 1.2 mL min -1 . Injection was done in split mode (split ratio of 1:8). The GC oven was programmed as follows: initial temperature 40°C (held for 4 min), then raised at 10°C min -1 until 150°C and then raised at 20°C min -1 until reaching a final temperature of 240°C (held for 8 min). MS transfer line and ion source temperatures were set to 250°C and 200°C, respectively. The ionization mode was electron impact (EI), 70 eV. Data acquisition was done in nine segments, including a full scan mode from 25 to 150 m/z at 0.2 s per scan and the corresponding selected ion monitoring (SIM) mode with a dwell time of 0.1 s each, except for furfural and d4-furfural 0.05 s each (Table 2). Chromatographic peak areas for each compound were calculated by extracting the quantifier ions from the SIM mode acquisition data using Quan Browser, Xcalibur 2.1.0 SP1, build 1160 (Thermo Fisher Scientific Inc., USA). Native molecules (i.e. generated in sponge cake during baking) were confirmed by means of the Wiley 8 and NIST 08 mass spectra libraries, calculation of normal alkane retention index (RI) and comparison to NIST Chemistry WebBook SRD 69 (2017) indices, and analysis of pure standards.

Sample preparation variability

In order to assess variability during batter preparation, baking and grinding, independent assays were run for reference formulation (three independent baking trials). Sponge cakes were sampled, frozen and ground as previously described to obtain a composite sample per baking trial. Each composite sample was then analyzed by HS trap / GC -MS in triplicate, both for furan and furfural quantitation, as well as aroma markers semiquantitation as detailed in section 2.7.

Data treatment

One-way ANOVA, Tukey-Kramer's honest significant difference (HSD) and Dunnett's tests were performed on dry matter, pH, color measurements as well as furan and furfural quantitative results. Reference formulation (f1) was defined as the control (Dunnett's test).

For the markers here considered, statistical differences among formulations were assessed by means of ANOVA and Tukey-Kramer's HSD test on raw chromatographic peak areas of quantifier ions, while principal component analysis (PCA) was carried out on the meancentered areas. All analyses were performed using JMP 10.0.0 (SAS Institute Inc.).

RESULTS AND DISCUSSION

Variability in sample preparation

Low variability was observed for furan and furfural quantitation on single sponge cake preparation experiments (≤ 3.5% for both compounds, n=3), as previously reported for the method [START_REF] Cepeda-Vázquez | Multivariate optimization of headspace trap for furan and furfural simultaneous determination in sponge cake[END_REF]. Experiment reproducibility assessed through three independent trials, including batter preparation, baking, grinding and HS trap analysis, both for furan and furfural quantitation was also low: 6.4 and 6.8%, respectively. Higher variability was observed for semi-quantitative analyses: 2-methylbutanal and 3methylbutanal (10%), hexanal (6%), 2-pentylfuran (11%), 1-hydroxy-2-propanone (11%), 2,5-dimethylpyrazine and 2,6-dimethylpyrazine (15%), 2,3,5-trimethylpyrazine (15%), benzaldehyde (8%). Yet, moderate variability was observed overall, showing good reproducibility throughout experiments.

Sample physical and physicochemical characterization

Low moisture content, neutral pH and mid-browning level (L*) were observed in general (Table 3). However, according to the ANOVA and Tukey-Kramer's test, significant differences were found among sponge cake formulations within each physical and physicochemical variable and were also grouped differently across them. Surprisingly, major differences were only found for glucose (f5) and egg white formulation (f6), exhibiting maximum and minimum values for pH, L* and b*. Indeed, glucose recipe showed the highest browning level (L*< 50) and the lowest pH (acid), as opposed to egg white formulation which exhibited the lowest browning level (L*> 50) and the highest pH (alkaline).

Low pH in formulation with glucose might be due to the formation of acidic compounds and correlated to advanced browning, as previously reported in cookies [START_REF] Ait Ameur | The fate of furfurals and other volatile markers during the baking process of a model cookie[END_REF]. Indeed, while a release of H + occurs during caramelization [START_REF] Kroh | Caramelisation in food and beverages[END_REF], organic acids are also formed during the Maillard reaction [START_REF] Martins | A review of Maillard reaction in food and implications to kinetic modelling[END_REF], which may lead to pH fall. Egg white's pH (7.6 to 7.9), which rises with temperature [START_REF] Belitz | Food Chemistry[END_REF], could explain higher pH in formulation f6. Low dry matter content in this ingredient (12.1%) [START_REF] Belitz | Food Chemistry[END_REF] may also explain low dry matter in this formulation, and in turn result in a lower browning level, given its higher water content.

Furan and furfural content in sponge cake

Furan content in sponge cake remains in the range of 3.1-12.5 ng g -1 , dry basis (Figure 1) which can be considered as a low to mid-level when compared to different types of bread (≤ 30 ng g -1 ), and similar or low when compared to other products, such as baby food containing fruit (≤ 16 ng g -1 ), cereal-based snacks (≤ 143 ng g -1 ) or coffee brew (≤ 199 ng g - 1 ) [START_REF] Zoller | Furan in food: Headspace method and product survey[END_REF]. Furfural content in sponge cake (0.29-0.83 μg g -1 , dry basis), except in the case of glucose recipe (9.2 μg g -1 , dry basis), is similar to levels observed for cake/pastry (0.26-2.6 μg g -1 ) (Petisca et al., 2014a).

Furan and aroma generation in sponge cake

A scheme based on current literature on the main reaction pathways (i.e. caramelization, Maillard reaction, Strecker's degradation and lipid oxidation) leading to the selected compounds (Table 2) and considering main precursors in sponge cake ingredients is here put forward (Figure 2). It will be discussed in the following sections to better understand how formulation impacts the generation of key quality markers in sponge cake.

Furan and furfural generation

Furfural can only be generated through caramelization or Maillard reaction [START_REF] Belitz | Food Chemistry[END_REF][START_REF] Kroh | Caramelisation in food and beverages[END_REF][START_REF] Martins | A review of Maillard reaction in food and implications to kinetic modelling[END_REF][START_REF] Srivastava | Kinetic study of furan and furfural generation during baking of cake models[END_REF][START_REF] Whitfield | Volatiles from interactions of Maillard reactions and lipids[END_REF][START_REF] Yaylayan | Origin of Carbohydrate Degradation Products in L -Alanine/D-[13C] Glucose Model Systems[END_REF], while furan can be additionally formed by PUFAs oxidation, as well as amino acid thermal degradation [START_REF] Crews | A review of the occurrence, formation and analysis of furan in heat-processed foods[END_REF][START_REF] Perez Locas | Origin and mechanistic pathways of formation of the parent furan -a food toxicant[END_REF][START_REF] Srivastava | Kinetic study of furan and furfural generation during baking of cake models[END_REF][START_REF] Stadler | Heat-generated toxicants in foods: acrylamide, MCPD esters and furan[END_REF] (Figure 2). Yet, our results show that furfural generation in sponge cake seems to be favored over furan's, since furfural content (ppb to ppm) was found to be 65 to 730-fold that of furan (ppb), depending on the formulation under study. Both furan and furfural content were mainly impacted by two categories of ingredients, namely sugar and egg-based (Figure 1).

The glucose-containing recipe (f5) exhibited their highest concentration, suggesting that the extent of non-enzymatic browning reactions (caramelization, Maillard reaction) in this formulation is greater than in the reference cake (f1) containing sucrose. This result is not surprising since glucose (reducing sugar) can directly undergo enolization and further degradation (through caramelization or Maillard reaction) (Figure 2) [START_REF] Kroh | Caramelisation in food and beverages[END_REF][START_REF] Martins | A review of Maillard reaction in food and implications to kinetic modelling[END_REF], while sucrose (non-reducing sugar) first has to undergo hydrolysis before releasing its reducing monomers (glucose and fructose) [START_REF] Quintas | Multiresponse modelling of the caramelisation reaction[END_REF]. Similar results were previously reported in sponge cake where sucrose replacement by reducing sugars led to higher amounts of furanic compounds (i.e. furfural and HMF (5-hydroxymethylfurfural)) [START_REF] Zhang | Effects of sugars in batter formula and baking conditions on 5-hydroxymethylfurfural and furfural formation in sponge cake models[END_REF]. Our results also show a higher molar yield regarding furan generation in the case of the glucose-containing formulation (f5) as compared to the sucrose recipe (f1) (18.0 and 13.0 ng g -1 mol -1 , dry basis, respectively). Opposite results were reported in model solutions (Perez Locas & [START_REF] Perez Locas | Origin and mechanistic pathways of formation of the parent furan -a food toxicant[END_REF], suggesting that sucrose undergoes hydrolysis less effectively in sponge cake than in simple models. Indeed, sucrose degradation into reducing sugars (fructose and glucose) may be limited in sponge cake [START_REF] Zhang | Effects of sugars in batter formula and baking conditions on 5-hydroxymethylfurfural and furfural formation in sponge cake models[END_REF]. This can be easily explained: model systems (single or binary mixtures with a high water content, heated at elevated temperatures, i.e. 250°C) do not reproduce what really occurs during baking of sponge cake (complex chemical composition and structure, where gradual water evaporation and drying of the surface lead to crust formation, whose temperature tends to that of the oven, i.e. 170°C) [START_REF] Fehaili | An instrumented oven for the monitoring of thermal reactions during the baking of sponge cake[END_REF]. This underlines the importance of performing studies in complex matrices.

Interestingly, the increase in furfural content in our recipe (f5) with respect to our reference cake (f1) was much greater than that in furan (around 17-fold and 3-fold, respectively). This suggests that, although glucose is a precursor of both compounds, the pathway leading to furfural (i.e. 1,2-enolization) is favored over the one majorly contributing to furan (i.e. 2,3enolization) [START_REF] Belitz | Food Chemistry[END_REF][START_REF] Kroh | Caramelisation in food and beverages[END_REF][START_REF] Martins | A review of Maillard reaction in food and implications to kinetic modelling[END_REF][START_REF] Perez Locas | Origin and mechanistic pathways of formation of the parent furan -a food toxicant[END_REF]. Indeed, in a recent kinetic study in cake models containing either glucose-only or glucose and leucine, furfural generation was also favored over furan's [START_REF] Srivastava | Kinetic study of furan and furfural generation during baking of cake models[END_REF].

The egg white formulation (f6) showed the lowest furan and furfural concentrations. Indeed, whole egg is richer in free amino acids (U.S. Department of Agriculture (USDA), 2016), which might enhance the Maillard reaction pathway or might even form furan through amino acid thermal degradation (Perez Locas & Yaylayan, 2004) (Figure 2). Besides, while PUFAs are present in egg yolk, phospholipids also are and might exert a protective effect by forming antioxidant species when reacting with carbonyl compounds [START_REF] Maire | Lipid oxidation in baked products : Impact of formula and process on the generation of volatile compounds[END_REF], hindering lipid oxidation as a result.

Interestingly, all formulations (f1-f3) designed to study furan formation via lipid oxidation had low and very similar furan contents, suggesting that furan formation through PUFA oxidation is not particularly favored. Likewise, in other products such as baby food model systems containing starch, sugars and soybean oil, a higher furan content was only obtained when oxidizing oil in extreme conditions (14 days at 60°C) prior to use [START_REF] Owczarek-Fendor | Furan formation in starch-based model systems containing carbohydrates in combination with proteins, ascorbic acid and lipids[END_REF]. Also, naturally occurring tocopherols in our fat ingredients (sunflower and palm oil) could exert an antioxidant effect as reported by [START_REF] Maire | Lipid oxidation in baked products : Impact of formula and process on the generation of volatile compounds[END_REF], hence preventing further furan generation through PUFAs oxidation. However, this hypothesis is here discarded since formulations with raw (f3) or purified sunflower oil (f7) resulted in no differences in furan content (7.1 ± 0.3 and 7.1 ± 0.2 ng g -1 , dry basis, respectively, P: 0.337), even if (f7) did show a higher amount of other lipid oxidation markers (i.e. hexanal and 2-pentylfuran).

Salt removal from recipe (f4) did not seem to exert any effect on furan formation (P: 0.0779), while it did on furfural, since a lower amount was found with respect to (f1) (P: 0.0005).

Previous studies on salt effect are scarce and results rather unclear. Undetectable furan levels were reported both in liquid models containing glucose or sucrose added with NaCl [START_REF] Fan | Factors affecting thermally induced furan formation[END_REF]. In the case of furfural, salt addition showed different effects depending on the complexity of the matrix: no effect on furfural generation in a glucosecontaining solution, and different effects on furfural content in cookies (increasing in the presence of leavening agents and decreasing in their absence) [START_REF] Kocadağlı | Effects of Sodium Chloride, Potassium Chloride, and Calcium Chloride on the Formation of α-Dicarbonyl Compounds and Furfurals and the Development of Browning in Cookies during Baking[END_REF]. Indeed, the complexity of the food matrix is to be considered (e.g. use of leavening agents).

To conclude, furan generation in sponge cake would seem to proceed through common pathways to furfural (i.e. caramelization and Maillard reaction). Their generation is enhanced in the presence of glucose (reducing sugar) and hindered in the absence of egg yolk (probably due to differences in Maillard reaction progression, given the restriction on free amino acids).

Impact of formulation on selected quality markers

PCA was carried out on all considered markers in order to determine correlations among them and relate them to changes in formulation. According to this analysis, components 1, 2 and 3 account for 95.5% of variability (Figure 3). Three highly correlated groups of compounds are observed in the PC1 vs PC2 loading plot: 1) furan, furfural and hexanal, 2) Strecker aldehydes, 2-pentylfuran and 1-hydroxypropanone, and 3) pyrazines. However, in the PC1 vs PC3 loading plot, hexanal and 1-hydroxypropanone exhibit a particular behavior: hexanal is opposite to furan and furfural, and 1-hydroxypropanone to Strecker aldehydes, when considering PC3. These results, along with possible reaction pathways (Figure 2), give rise to several remarks. First of all, furan strongly correlates to furfural and overall behaves differently from all the other volatile markers. 2-Pentylfuran and hexanal are both well-known lipid oxidation products [START_REF] Belitz | Food Chemistry[END_REF][START_REF] Cerny | Thermal Generation of Aroma-Active Volatiles in Food[END_REF][START_REF] Perez Locas | Origin and mechanistic pathways of formation of the parent furan -a food toxicant[END_REF]) generated from linoleic and arachidonic acid via 2,4-decadienal [START_REF] Grein | 4-Hydroxy-2-alkenals and other products formed by water-mediated oxidative decomposition of α,β-unsaturated aldehydes[END_REF][START_REF] Whitfield | Volatiles from interactions of Maillard reactions and lipids[END_REF]; they strongly correlate when considering PC2, but not regarding PC1. This suggests that additional pathways could be active in the formation of 2-pentylfuran in sponge cake, either through cyclization of other alkenals [START_REF] Adams | Amino Acid Catalysis of 2-Alkylfuran Formation from Lipid Oxidation-Derived α,β-Unsaturated Aldehydes[END_REF][START_REF] Hidalgo | Strecker Type Degradation of Phenylalanine by 4-Hydroxy-2-nonenal in Model Systems[END_REF] or via a lipid oxidation and amino group interaction [START_REF] Adams | Amino Acid Catalysis of 2-Alkylfuran Formation from Lipid Oxidation-Derived α,β-Unsaturated Aldehydes[END_REF][START_REF] Whitfield | Volatiles from interactions of Maillard reactions and lipids[END_REF]. It should also be noted that hexanal may be formed prior to sponge cake baking, that is during the beating of batter, as previously shown by [START_REF] Maire | Lipid oxidation in baked products : Impact of formula and process on the generation of volatile compounds[END_REF].

Pyrazines are produced through a well-known pathway involving aminoketone degradation [START_REF] Cerny | Thermal Generation of Aroma-Active Volatiles in Food[END_REF][START_REF] Whitfield | Volatiles from interactions of Maillard reactions and lipids[END_REF]. This very specific origin is well represented on the PCA loading plots where they are closely correlated to one another (Figure 2). Benzaldehyde would be produced either by Strecker's degradation or lipid oxidation, as previously reported by [START_REF] Whitfield | Volatiles from interactions of Maillard reactions and lipids[END_REF], explaining its high correlation both to other Strecker aldehydes and 2-pentylfuran.

1-Hydroxy-2-propanone is an especially interesting marker. It could be produced either through sugar degradation [START_REF] Yaylayan | Origin of Carbohydrate Degradation Products in L -Alanine/D-[13C] Glucose Model Systems[END_REF] or Maillard reaction [START_REF] Belitz | Food Chemistry[END_REF][START_REF] Martins | A review of Maillard reaction in food and implications to kinetic modelling[END_REF]) via 1-deoxy-2,3-hexodiulose, a common intermediate to furan, but also from aminoketones produced during Strecker's pathway [START_REF] Whitfield | Volatiles from interactions of Maillard reactions and lipids[END_REF].

Moreover, as shown in simple models, 1-hydroxy-2-pronanone could yield pyrazines in the presence of NH3 along with aliphatic aldehydes (e.g. hexanal) and these compounds may be formed during cake preparation (through Strecker degradation from cysteine and lipid oxidation, respectively [START_REF] Whitfield | Volatiles from interactions of Maillard reactions and lipids[END_REF]). All of the above would explain 1-hydroxy-2propanone's particular location in the PCA loading plots (Figure 3), correlated to Strecker aldehydes, positively influencing PC1 (like pyrazines, unlike furan) and PC2 (like furan, unlike pyrazines).

The effect of sponge cake formulation is clearly visible in the PC1 vs PC2 plane (85%). Indeed, four clusters of recipes are observed in the score plot: glucose (f5, top left), egg white (f6, bottom left), with no added fat nor salt (f2, f4, center) and reference and sunflower oil (f1, f3, right center). This indicates that in the presence of glucose, the peak intensities of furan, furfural and hexanal are higher, and those of the remaining compounds are lower.

However, given hexanal's particular behavior in the PC1 vs PC3 plane, it can only be concluded that the presence of glucose might induce a quicker progression of caramelization and Maillard reaction in sponge cake, yielding higher amounts of furan and furfural and limiting the formation of compounds related to amino acid degradation (mainly Strecker-related and pyrazines). In the absence of egg yolk (f6), the peak intensities of nearly all compounds are lower. In this case, the restriction on some precursors (e.g. amino acids) might be responsible for limiting the progression of certain pathways, in particular those related to amino acid degradation. Indeed, our work confirms previous findings where a higher response of 2-methylbutanal, 3-methylbutanal, benzaldehyde and 2-pentylfuran were observed in a whole egg sponge cake recipe, compared to an egg white-only formulation [START_REF] Maire | Lipid oxidation in baked products : Impact of formula and process on the generation of volatile compounds[END_REF].

All considered markers were significantly influenced by formulation as shown by the ANOVA procedure, confirming PCA results (see supplementary table). Overall, sponge cake reference formulation containing sucrose and whole egg has a greater variety in aroma compounds such as Strecker aldehydes and pyrazines; glucose recipe increases furan and furfural response, most likely due to this sugar's reducing nature; and egg white formulation causes a decrease in the peak intensity of aroma compounds, probably related to amino acid-related pathways.

Clearly, an appropriate selection of key markers, a suitable experimental design and the use of multivariate analysis tools make it possible to have better understanding of the link between formulation and reaction pathways in real food, as in this case with sponge cake.

CONCLUSIONS

Furan and furfural content, as well as the generation of key aroma markers from different reaction pathways taking place during baking, are clearly influenced by sponge cake formulation. Sugar and egg-based are undoubtedly the most interesting categories of ingredients involved in their formation in this product.

Furan content is mainly affected by the presence of glucose (increase) and the absence of egg yolk (decrease). Reaction pathways active in its formation in sponge cake are most likely caramelization, Maillard reaction or even amino acid degradation, but not lipid oxidation. Furfural is highly correlated to furan, but its higher content in all recipes suggests that its generation is favored over furan's.

Key aroma compounds (Strecker aldehydes, pyrazines and 2-pentylfuran) are positively related to formulation containing whole egg, probably due to a higher contribution of free amino acids in egg yolk.

This work represents a step forward in delivering baked goods having optimal food safety and aroma properties. It contributes to the understanding of the two-way relationship among key quality markers and formulation. Our findings will be useful for developing furan mitigation strategies based on formulation and also for understanding how changes in recipes (e.g. sucrose replacement, low fat, fat free or salt free recipes) might affect the generation of key compounds in these bakery products. [START_REF] Perez Locas | Origin and mechanistic pathways of formation of the parent furan -a food toxicant[END_REF][START_REF] Srivastava | Kinetic study of furan and furfural generation during baking of cake models[END_REF]d Stadler, 2012;e Belitz, Grosch, & Schieberle, 2009;f Cerny, 2010;g Whitfield, 1992;h Adams, Bouckaert, Lancker, Meulenaer, & Kimpe, 2011;i Grein, Huffer, Scheller, & Schreier, 1993;[START_REF] Hidalgo | Strecker Type Degradation of Phenylalanine by 4-Hydroxy-2-nonenal in Model Systems[END_REF]k Martins, Jongen, & van Boekel, 2001;l Yaylayan & Keyhani, 2000;m Kroh, 1994;n [START_REF] Crews | A review of the occurrence, formation and analysis of furan in heat-processed foods[END_REF][START_REF] Perez Locas | Origin and mechanistic pathways of formation of the parent furan -a food toxicant[END_REF][START_REF] Srivastava | Kinetic study of furan and furfural generation during baking of cake models[END_REF][START_REF] Stadler | Heat-generated toxicants in foods: acrylamide, MCPD esters and furan[END_REF][START_REF] Belitz | Food Chemistry[END_REF][START_REF] Cerny | Thermal Generation of Aroma-Active Volatiles in Food[END_REF][START_REF] Whitfield | Volatiles from interactions of Maillard reactions and lipids[END_REF][START_REF] Adams | Amino Acid Catalysis of 2-Alkylfuran Formation from Lipid Oxidation-Derived α,β-Unsaturated Aldehydes[END_REF][START_REF] Grein | 4-Hydroxy-2-alkenals and other products formed by water-mediated oxidative decomposition of α,β-unsaturated aldehydes[END_REF][START_REF] Hidalgo | Strecker Type Degradation of Phenylalanine by 4-Hydroxy-2-nonenal in Model Systems[END_REF][START_REF] Martins | A review of Maillard reaction in food and implications to kinetic modelling[END_REF][START_REF] Yaylayan | Origin of Carbohydrate Degradation Products in L -Alanine/D-[13C] Glucose Model Systems[END_REF][START_REF] Kroh | Caramelisation in food and beverages[END_REF] Figure 3 Principal component analysis on selected volatile compounds in sponge cake formulations. f1: reference; f2: no added fat; f3: sunflower oil; f4: no salt; f5: glucose; f6: egg white; 2MB: 2-methylbutanal; 3MB: 3-methylbutanal; 2PF: 2-pentylfuran; 1H2P: 1-hydroxy-2-propanone; 2,5 DMP: 2,5-dimethylpyrazine; 2,6 DMP: 2,6dimethylpyrazine; 2,3,5 TMP: 2,3,5-trimethylpyrazine Supplementary

  oil; SO: sunflower oil; pSO: purified sunflower oil; S: sucrose; G: glucose; E: whole egg; EW: egg white; PUFAs: polyunsaturated fatty acids

Figure 1 Figure 2

 12 Figure 1 Furan (ng/g) and furfural (µg/g) content in sponge cake according to formulation. (a) furan and furfural content in formulations 1 to 6; (b) formulations having the lowest furfural content; f1: reference; f2: no added fat; f3: sunflower oil; f4: no salt; f5: glucose; f6: egg white; values are the means and 95% confidence intervals (n=3); mean values with different letters are significantly different (P: <0.05)

Table 1 Sponge cake formulation according to ingredient categories Ingredients (g)

 1

	Formulation description	Testing for	Flour	Fat	Salt	Sugar	Egg-based	Total
				Type 55 PO	SO	pSO NaCl	S	G	E	EW
			effect of fat type and amount on lipid
		reference (palm oil,	oxidation; effect of salt and sugar type
	f1	salt, sucrose, whole	on browning reactions; effect of egg-
		egg)	based ingredients on aroma
			generation

Table 2 Selected volatile markers in sponge cake, associated pathways, sensory attributes and identification method

 2

	Compound	Class	Reaction pathway Flavor/odor	CAS No.	RIref.*	RIexp.	Selected ions (m/z)**	Identification
	furan	furan	MR, C, LO, AA	-	110-00-9	813 p	811	39, 68	MS, RI, STD
			degradation a, b, c, d
	2-methylbutanal	aldehyde	S e, f, g	almond, cocoa, fermented,	96-17-3	930 p	928	41, 44, 57	MS, RI, STD
				hazelnut, malty n, o
	3-methylbutanal	aldehyde	S e, f, g	malty, roasty cucumber-like o	590-86-3	934 p	933	41, 44, 57	MS, RI, STD
	hexanal	aldehyde	LO e, f, g	apple, fat, fresh, green, oil,	66-25-1	1103 p	1099	44, 56	MS, RI, STD
				grassy, tallow n, o
	2-pentylfuran	furan	LO, LO & AA	butter, floral, fruit, green bean,	3777-69-3	1240 p	1248	81, 138	MS, RI, STD
			interaction b, e, g, h, i, j	mushroom, raw nuts n, o
	1-hydroxy-2-propanone hydroxyketone MR, C, S e, g, k, l	pungent, sweet caramellic,	116-09-6	1319 p	1345	43, 74	MS, RI, STD
				ethereal o
	2,5-dimethylpyrazine	pyrazine	S interaction g	cocoa, roast beef, roasted nut,	123-32-0	1357 p	1357	42, 108	MS, RI, STD
				crust-like, popcorn n, o
	2,6-dimethylpyrazine	pyrazine	S interaction g	cocoa, coffee, green, roast beef,	108-50-9	1362 p	1362	42, 108	MS, RI
				roasted nut, roasted n, o
	2,3,5-trimethylpyrazine pyrazine	S interaction f, g	cocoa, earthy, must, potato,	14667-55-1	1429 p	1439	42, 122	MS, RI, STD
				roast n, o
	furfural	furan	MR, C e, g, k, l, m	almond, baked potatoes, bread,	98-01-1	1503 p	1507	39, 95, 96	MS, RI, STD
				burnt, spice, soil, roasted,
				sweet, toasted n, o
	benzaldehyde	aldehyde	S, LO g	bitter almond, burnt sugar,	100-52-7	1533 p	1587	77, 106	MS, RI, STD
				cherry, malt, roasted pepper,
				caramel n, o
	Internal standards
	d4-furan				6142-90-1	-	-	42, 72	-
	d4-furfural				1219803-80-1	-	-42, 98, 100	-
	MR: Maillard reaction; C: caramelization; LO: lipid oxidation; AA: amino acid; S: Strecker's degradation; RI: normal alkane retention index; MS: mass spectra; STD: pure standard
	a Crews & Castle, 2007; b Perez Locas &

  Flavor & Extract Manufacturers Association (FEMA), 2017; o Pico, Bernal, & Gómez, 2015; p National Institute of Standards and Technology (NIST), 2017 *Retention indices on a ZB-WAX-like column; **Quantifier ions are in bold

Table 3 Sponge cake physical and physicochemical characterization

 3 Values are the means from 3 replicates for dry matter and pH, and 4 replicates for color measurements. Values followed by different letters within columns are significantly different (P: <0.05)

	Formulation	dry matter (%)	pH	L*	a*	b*
	f1	92.8 a	7.0 d 53.2 c 12.8 a 15.5 c
	f2	90.0 bc	7.1 d 56.3 c 13.7 a 18.6 bc
	f3	89.5 c	7.2 c 55.5 c 14.7 a 19.0 bc
	f4	93.6 a	7.7 b 62.6 b 12.9 a 21.7 b
	f5	91.3 b	5.5 e 43.9 d	6.2 b	4.6 d
	f6	87.3 d	8.0 a 70.1 a 12.9 a 26.5 a

Table Peak areas of the selected markers in sponge cake

  Values are the means of raw chromatographic peak areas from 3 replicates. Values followed by different letters within columns are significantly different (P<0.05)

	Formulation	furan	2MB+3MB	hexanal	2PF	1H2P	2,5 DMP+2,6 DMP 2,3,5 TMP	furfural	benzaldehyde
	1	1019293 b	117334872 a 1036261 d 6123387 b	793669 a	2165852 a 401325 a	983668 bc	1737559 a
	2	897436 cd	80511591 c 3833174 c 2383420 d	621995 bcd	1791143 b 280398 bc	1332845 b	1364217 c
	3	969202 bc	91631525 b 1073709 d 7135143 a	746201 ab	1924555 b 308601 b	927058 bc	1516750 bc
	4	869030 d	95836968 b 4831596 b 6144870 b	476103 d	1458376 c 245609 c	598459 c	1632867 ab
	5	2291855 a	110936634 a 5566375 a 4429624 c	641014 abc	35616 e	21244 e	16721826 a	1029226 d
	6	506805 e	9527479 d	426097 e	737079 e	488099 cd	1146420 d	98061 d	468726 c	384190 e
	2MB: 2-methylbutanal; 3MB: 3-methylbutanal; 2PF: 2-pentylfuran; 1H2P: 1-hydroxy-2-propanone; 2,5 DMP: 2,5-dimethylpyrazine; 2,6 DMP: 2,6-dimethylpyrazine; 2,3,5 TMP:
	2,3,5-trimethylpyrazine

ACKNOWLEDGEMENTS

The authors wish to thank Paola Soto for technical support in oil purification and running HPLC analysis, and also David Blumenthal and Pierre Giampaoli for fruitful discussions.

The authors have declared no conflict of interest.

FUNDING

This work was supported by the Mexican National Council for Science and Technology (Consejo Nacional de Ciencia y Tecnología -CONACYT) and the French Embassy in Mexico (joint grant: 314446).