Brownian bridge expansions for L\'evy area approximations and particular values of the Riemann zeta function

We study approximations for the L\'evy area of Brownian motion which are based on the Fourier series expansion and a polynomial expansion of the associated Brownian bridge. Comparing the asymptotic convergence rates of the L\'evy area approximations, we see that the approximation resulting from the polynomial expansion of the Brownian bridge is more accurate than the Kloeden-Platen-Wright approximation, whilst still only using independent normal random vectors. We then link the asymptotic convergence rates of these approximations to the limiting fluctuations for the corresponding series expansions of the Brownian bridge. Moreover, and of interest in its own right, the analysis we use to identify the fluctuation processes for the Karhunen-Lo\`eve and Fourier series expansions of the Brownian bridge is extended to give a stand-alone derivation of the values of the Riemann zeta function at even positive integers.


Introduction
One of the well-known applications for expansions of the Brownian bridge is the strong or L 2 (P) approximation of stochastic integrals. Most notably, the second iterated integrals of Brownian motion are required by high order strong numerical methods for general stochastic differential equations (SDEs), as discussed in [4,22,33]. Due to integration by parts, such integrals can be expressed in terms of the increment and Lévy area of Brownian motion. The approximation of multidimensional Lévy area is well-studied, see [5,8,11,13,14,23,25,32,35], with the majority of the algorithms proposed being based on a Fourier series expansion or the standard piecewise linear approximation of Brownian motion. Some alternatives include [5,11,25] which consider methods associated with a polynomial expansion of the Brownian bridge.
Since the advent of Multilevel Monte Carlo (MLMC), introduced by Giles in [16] and subsequently developed in [2,6,7,15,17], Lévy area approximation has become less prominent in the literature. In particular, the antithetic MLMC method introduced by Giles and Szpruch in [17] achieves the optimal complexity for the weak approximation of multidimensional SDEs without the need to generate Brownian Lévy area. That said, there are concrete applications where the simulation of Lévy area is beneficial, such as for sampling from non-log-concave distributions using Itô diffusions. For these sampling problems, high order strong convergence properties of the SDE solver lead to faster mixing properties of the resulting Markov chain Monte Carlo (MCMC) algorithm, see [26].
In this paper, we compare the approximations of Lévy area based on the Fourier series expansion and on a polynomial expansion of the Brownian bridge. We particularly observe their convergence rates and link those to the fluctuation processes associated with the different expansions of the Brownian bridge. The fluctuation process for the polynomial expansion is studied in [19], and our study of the fluctuation process for the Fourier series expansion allows us, at the same time, to determine the fluctuation process for the Karhunen-Loève expansion of the Brownian bridge. As an attractive side result, we extend the required analysis to obtain a stand-alone derivation of the values of the Riemann zeta function at even positive integers. Throughout, we denote the positive integers by N and the non-negative integers by N 0 .
Let us start by considering a Brownian bridge (B t ) t∈ [0,1] in R with B 0 = B 1 = 0. This is the unique continuous-time Gaussian process with mean zero and whose covariance function K B is given by, for s, t ∈ [0, 1], We are concerned with the following three expansions of the Brownian bridge. The Karhunen-Loève expansion of the Brownian bridge, see Loève [27, p. 144], is of the form, for t ∈ [0, 1], 2 sin(kπt) kπ 1 0 cos(kπr) dB r .
The Fourier series expansion of the Brownian bridge, see p. 198] or Kahane [21,Sect. 16.3], yields, for t ∈ [0, 1], where, for k ∈ N 0 , A polynomial expansion of the Brownian bridge in terms of the shifted Legendre polynomials Q k on the interval [0, 1] of degree k, see [12,19], is given by, for t ∈ [0, 1], where, for k ∈ N, These expansions are summarised in Table 1 in Appendix A and they are discussed in more detail in Section 2. For an implementation of the corresponding approximations for Brownian motion as Chebfun examples into MATLAB, see Filip, Javeed and Trefethen [9] as well as Trefethen [34].
We remark that the polynomial expansion (1.5) can be viewed as a Karhunen-Loève expansion of the Brownian bridge with respect to the weight function w on (0, 1) given by w(t) = 1 t(1−t) . This approach is employed in [12] to derive the expansion along with the standard optimality property of Karhunen-Loève expansions. In this setting, the polynomial approximation of (B t ) t∈[0,1] is optimal among truncated series expansions in a weighted L 2 (P) sense corresponding to the non-constant weight function w. To avoid confusion, we still adopt the convention throughout to reserve the term Karhunen-Loève expansion for (1.2), whereas (1.5) will be referred to as the polynomial expansion.
Before we investigate the approximations of Lévy area based on the different expansions of the Brownian bridge, we first analyse the fluctuations associated with the expansions. The fluctuation process for the polynomial expansion is studied and characterised in [19], and these results are (a k cos(2kπt) + b k sin(2kπt)) . 1] is the natural scaling to use because increasing N by one results in the subtraction of two additional Gaussian random variables. We use E to denote the expectation with respect to Wiener measure P.

The scaling by
The fluctuation processes (F N,2 t ) t∈[0,1] for the Fourier expansion converge in finite dimensional distributions as N → ∞ to the collection (F 2 t ) t∈[0,1] of zero-mean Gaussian random variables whose covariance structure is given by, for s, t ∈ [0, 1], The difference between the fluctuation result for the Karhunen-Loève expansion and the fluctuation result for the polynomial expansion, see [19,Theorem 1.6] or Section 2.3, is that there the variances of the independent Gaussian random variables follow the semicircle 1 π t(1 − t) whereas here they are constant on (0, 1), see Figure 1. The limit fluctuations for the Fourier series expansion further exhibit endpoints which are correlated.
As pointed out in [19], the reason for considering convergence in finite dimensional distributions for the fluctuation processes is that the limit fluctuations neither have a realisation as processes in C([0, 1], R), nor are they equivalent to measurable processes.
We prove Theorem 1.1 by studying the covariance functions of the Gaussian processes (F N,1 1] given in Lemma 2.2 and Lemma 2.3 in the limit N → ∞. The key ingredient is the following limit theorem for sine functions, which we see concerns the pointwise convergence for the covariance function of (F N,1 The above result serves as one of four base cases in the analysis performed in [18] of the asymptotic error arising when approximating the Green's function of a Sturm-Liouville problem through a truncation of its eigenfunction expansion. The work [18] offers a unifying view for Theorem 1.  The proof of Theorem 1.2 is split into an on-diagonal and an off-diagonal argument. We start by proving the convergence on the diagonal away from its endpoints by establishing locally uniform convergence, which ensures continuity of the limit function, and by using a moment argument to identify the limit. As a consequence of the on-diagonal convergence, we obtain the next corollary which then implies the off-diagonal convergence in Theorem 1.2.
Moreover, and of interest in its own right, the moment analysis we use to prove the on-diagonal convergence in Theorem 1.2 leads to a stand-alone derivation of the result that the values of the Riemann zeta function ζ : C \ {1} → C at even positive integers can be expressed in terms of the Bernoulli numbers B 2n as, for n ∈ N, see Borevich and Shafarevich [3]. In particular, the identity that is, the resolution to the Basel problem posed by Mengoli [28] is a consequence of our analysis and not a prerequisite for it.
We turn our attention to studying approximations of second iterated integrals of Brownian motion.
For an illustration of Lévy area for a two-dimensional Brownian motion, see Figure 2.
Remark 1.5. Given the increment W t − W s and the Lévy area A s,t , we can recover the second iterated integrals of Brownian motion using integration by parts as, for i, j ∈ {1, . . . , d} with i = j, ሺ2ሻ ሺ1ሻ Figure 2. Lévy area is the chordal area between independent Brownian motions.
We consider the sequences {a k } k∈N0 , {b k } k∈N and {c k } k∈N of Gaussian random vectors, where the coordinate random variables a  [22] and Milstein [31]. Further approximating terms so that only independent random coefficients are used yields the Kloeden-Platen-Wright approximation in [23,30,35]. Similarly, using the random coefficients from the polynomial expansion (1.5), we obtain the Lévy area approximation first proposed by Kuznetsov in [24]. These Lévy area approximations are summarised in Table 2 in Appendix A and have the following asymptotic convergence rates. Theorem 1.6 (Asymptotic convergence rates of Lévy area approximations). For n ∈ N, we set N = 2n and define approximations A n , A n and s A 2n of the Lévy area A 0,1 by, for i, j ∈ {1, . . . , d}, Then A n , A n and s A 2n are antisymmetric d × d matrices and, for i = j and as N → ∞, we have The asymptotic convergence rates in Theorem 1.6 are phrased in terms of N since the number of Gaussian random vectors required to define the above Lévy area approximations is N or N − 1, respectively. Of course, it is straightforward to define the polynomial approximation s A n for n ∈ N, see Theorem 5.4.
Intriguingly, the convergence rates for the approximations resulting from the Fourier series and the polynomial expansion correspond exactly with the areas under the limit variance function for each fluctuation process, which are We provide heuristics demonstrating how this correspondence arises at the end of Section 5.
By adding an additional Gaussian random matrix that matches the covariance of the tail sum, it is possible to derive high order Lévy area approximations with O(N −1 ) convergence in L 2 (P).
Wiktorsson [35] proposed this approach using the Kloeden-Platen-Wright approximation (1.11) and this was recently improved by Mrongowius and Rößler in [32] who use the approximation (1.10) obtained from the Fourier series expansion (1.3).
We expect that an O(N −1 ) polynomial-based approximation is possible using the same techniques. While this approximation should be slightly less accurate than the Fourier approach, we expect it to be easier to implement due to both the independence of the coefficients {c k } k∈N and the covariance of the tail sum having a closed-form expression, see Theorem 5.4. Moreover, this type of method has already been studied in [5,10,11] with Brownian Lévy area being approximated by where the antisymmetric d × d matrix λ 0,1 is normally distributed and designed so that Û A 0,1 has the same covariance structure as the Brownian Lévy area A 0,1 . Davie [5] as well as Flint and Lyons [10] generate each (i, j)-entry of λ 0,1 independently as λ (i,j) 0,1 ∼ N 0, 1 12 for i < j . In [11], it is shown that the covariance structure of A 0,1 can be explicitly computed conditional on both W 1 and c 1 . By matching the conditional covariance structure of A 0,1 , the work [11] obtains the approximation where the entries λ (i,j) 0,1 i < j are still generated independently, but only after c 1 has been generated.
By rescaling (1.13) to approximate Lévy area on k N , k+1 N and summing over k ∈ {0, . . . , N −1}, we obtain a fine discretisation of A 0,1 involving 2N Gaussian random vectors and N random matrices. In [5,10,11], the Lévy area of Brownian motion and this approximation are probabilistically coupled in such a way that L 2 (P) convergence rates of O(N −1 ) can be established. Furthermore, the efficient Lévy area approximation (1.13) can be used directly in numerical methods for SDEs, which then achieve L 2 (P) convergence of O(N −1 ) under certain conditions on the SDE vector fields, see [5,10]. We leave such high order polynomial-based approximations of Lévy area as a topic for future work.
The paper is organised as follows. In Section 2, we provide an overview of the three expansions we consider for the Brownian bridge, and we characterise the associated fluctuation processes 1] . Before discussing their behaviour in the limit N → ∞, we initiate the moment analysis used to prove the on-diagonal part of Theorem 1.2 and we extend the analysis to determine the values of the Riemann zeta function at even positive integers in Section 3. The proof of Theorem 1.2 follows in Section 4, where we complete the moment analysis and establish a locally uniform convergence to identify the limit on the diagonal, before we deduce Corollary 1.3, which then allows us to obtain the off-diagonal convergence in Theorem 1.2. We close Section 4 by proving Theorem 1.1. In Section 5, we compare the asymptotic convergence rates of the different approximations of Lévy area, which results in a proof of Theorem 1.6.

Series expansions for the Brownian bridge
We discuss the Karhunen-Loève expansion as well as the Fourier expansion of the Brownian bridge more closely, and we derive expressions for the covariance functions of their Gaussian fluctuation processes.
In our analysis, we frequently use a type of Itô isometry for Itô integrals with respect to a Brownian bridge, and we include its statement and proof for completeness.
has the same law as the Brownian bridge (B t ) t∈[0,1] . In particular, the random variable Using a similar expression for 1 0 g(t) dB t and applying the usual Itô isometry, we deduce that as claimed.
2.1. The Karhunen-Loève expansion. Mercer's theorem, see [29], states that for a continuous symmetric non-negative definite kernel K : which consists of eigenfunctions of the Hilbert-Schmidt integral operator associated with K and whose eigenvalues {λ k } k∈N are non-negative and such that, for s, t ∈ [0, 1], we have the representation which converges absolutely and uniformly on The Karhunen-Loève expansion of the Brownian bridge is then given by  1] for N ∈ N is a zero-mean Gaussian process with covariance function Proof. From the definition (1.7), we see that (F N,1 t ) t∈[0,1] is a zero-mean Gaussian process. Hence, it suffices to determine its covariance function. By Lemma 2.1, we have, for k, l ∈ N, Therefore, from (1.1) and (1.7), we obtain that, for all s, t ∈ [0, 1], as claimed.
Consequently, Theorem 1.2 is a statement about the pointwise convergence of the function N C N 1 in the limit N → ∞.
For our stand-alone derivation of the values of the Riemann zeta function at even positive integers in Section 3, it is further important to note that since, by Mercer's theorem, the representation Applying Lemma 2.1, we see that and, for k, l ∈ N, Since the random coefficients are Gaussian random variables with mean zero, by (2.3) and (2.4), this implies that, for k ∈ N, For the remaining covariances of these random coefficients, we obtain that, for k, l ∈ N, Using the covariance structure of the random coefficients, we determine the covariance functions of the fluctuation processes (F N,2 t ) t∈[0,1] defined in (1.8) for the Fourier series expansion. 1] for N ∈ N is a Gaussian process with mean zero and whose covariance function is Proof. Repeatedly applying Lemma 2.1, we compute that, for t ∈ [0, 1], as well as, for k ∈ N, From (2.5) and (2.8), it follows that, for s, t ∈ [0, 1], whereas (2.7) and (2.9) imply that as well as It remains to observe that, by (2.6) and (2.7), cos(2kπs) cos(2kπt) + sin(2kπs) sin(2kπt) 2k 2 π 2 .
By combining Corollary 1.3, the resolution (1.9) to the Basel problem and the representation (2.1), we can determine the pointwise limit of 2N C N 2 as N → ∞. We leave further considerations until Section 4.2 to demonstrate that the identity (1.9) is really a consequence of our analysis.
2.3. The polynomial expansion. As pointed out in the introduction and as discussed in detail in [12], the polynomial expansion of the Brownian bridge is a type of Karhunen-Loève expansion in the weighted L 2 (P) space with weight function w on (0, 1) defined by w(t) = 1 t(1−t) . An alternative derivation of the polynomial expansion is given in [19] by considering iterated Kolmogorov diffusions. The iterated Kolmogorov diffusion of step N ∈ N pairs a one-dimensional Brownian motion (W t ) t∈[0,1] with its first N − 1 iterated time integrals, that is, it is the stochastic process in R N of the form The shifted Legendre polynomial Q k of degree k ∈ N on the interval [0, 1] is defined in terms of the standard Legendre polynomial P k of degree k on [−1, 1] by, for t ∈ [0, 1], It is then shown that the first component of an iterated Kolmogorov diffusion of step N ∈ N conditioned to return to 0 ∈ R N in time 1 has the same law as the stochastic process The polynomial expansion (1.5) is an immediate consequence of the result [19, Theorem 1.4] which states that these first components of the conditioned iterated Kolmogorov diffusions converge weakly as N → ∞ to the zero process.
As for the Karhunen-Loève expansion discussed above, the sequence {c k } k∈N of random coefficients defined by (1.6) is again formed by independent Gaussian random variables. To see this, we first recall the following identities for Legendre polynomials [1, (12.23), (12.31), (12.32)] which in terms of the shifted Legendre polynomials read as, for k ∈ N, In particular, it follows that, for all k ∈ N, which, by Lemma 2.1, implies that, for k, l ∈ N, Since the random coefficients are Gaussian with mean zero, this establishes their independence.
The fluctuation processes (F N,3 t ) t∈[0,1] for the polynomial expansion defined by are studied in [19]. According to [19,Theorem 1.6], they converge in finite dimensional distributions as N → ∞ to the collection (F 3 t ) t∈[0,1] of independent Gaussian random variables with mean zero and variance that is, the variance function of the limit fluctuations is given by a scaled semicircle.

Particular values of the Riemann zeta function
We demonstrate how to use the Karhunen-Loève expansion of the Brownian bridge or, more precisely, the series representation arising from Mercer's theorem for the covariance function of the Brownian bridge to determine the values of the Riemann zeta function at even positive integers. The analysis further feeds directly into Section 4.1 where we characterise the limit fluctuations for the Karhunen-Loève expansion.
The crucial ingredient is the observation (2.2) from Section 2, which implies that, for all n ∈ N 0 , .
For completeness, we recall that the Riemann zeta function ζ : C \ {1} → C analytically continues the sum of the Dirichlet series When discussing its values at even positive integers, we encounter the Bernoulli numbers. The Bernoulli numbers B n , for n ∈ N, are signed rational numbers defined by an exponential generating function via, for t ∈ (−2π, 2π), In particular, choosing m = 1 yields 1 + 2B 1 = 0, which shows that Moreover, since the function defined by, for t ∈ (−2π, 2π), is an even function, we obtain B 2n+1 = 0 for all n ∈ N, see [3, Theorem 5.8.2]. It follows from (3.2) that the Bernoulli numbers B 2n indexed by even positive integers are uniquely characterised by the recurrence relations These recurrence relations are our tool for identifying the Bernoulli numbers when determining the values of the Riemann zeta function at even positive integers.
The starting point for our analysis is (3.1), and we first illustrate how it allows us to compute ζ(2). Taking n = 0 in (3.1), multiplying through by π 2 , and using that 1 0 (sin(kπt)) 2 dt = 1 2 for k ∈ N, we deduce that We observe that this is exactly the identity obtained by applying the general result For working out the values for the remaining even positive integers, we iterate over the degree of the moment in (3.1). While for the remainder of this section it suffices to only consider the even moments, we derive the following recurrence relation and the explicit expression both for the even and for the odd moments as these are needed in Section 4.1. For k ∈ N and n ∈ N 0 , we set e k,n = 1 0 2 (sin(kπt)) 2 t n dt . Proof. For k ∈ N, the values for e k,0 and e k,1 can be verified directly. For n ∈ N with n ≥ 2, we integrate by parts twice to obtain e k,n = 1 0 2 (sin(kπt)) 2 t n dt as claimed.
Iteratively applying the recurrence relation, we find the following explicit expression, which despite its involvedness is exactly what we need. Proof. We proceed by induction over m. Since e k,0 = 1 and e k,1 = 1 2 for all k ∈ N, the expressions are true for m = 0 with the sums being understood as empty sums in this case. Assuming that the result is true for some fixed m ∈ N 0 , we use Lemma 3.1 to deduce that e k,2m+2 = 1 2m + 3 − (2m + 2)(2m + 1) 4k 2 π 2 e k,2m as well as which settles the induction step.

From Lemma 3.2, it follows that
.
Since ∞ k=1 k −2n converges for all n ∈ N, we can rearrange sums to obtain which in terms of the Riemann zeta function and after reindexing the sum rewrites as Multiplying through by (2m + 1)(2m + 2)(2m + 3) shows that, for all m ∈ N 0 , m+1 n=1 Ç 2m + 3 2n Comparing the last expression with the characterisation (3.3) of the Bernoulli numbers B 2n indexed by even positive integers implies that that is, we have established that, for all n ∈ N,

Fluctuations for the trigonometric expansions of the Brownian bridge
We first prove Theorem 1.2 and Corollary 1.3 which we use to determine the pointwise limits for the covariance functions of the fluctuation processes for the Karhunen-Loève expansion and of the fluctuation processes for the Fourier series expansion, and then we deduce Theorem 1.1.

Fluctuations for the Karhunen-Loève expansion.
For the moment analysis initiated in the previous section to allow us to identify the limit of N C N 1 as N → ∞ on the diagonal away from its endpoints, we apply the Arzelà-Ascoli theorem to guarantee continuity of the limit away from the endpoints. To this end, we first need to establish the uniform boundedness of two families of functions. Recall that the functions C In particular, for all N ∈ N and all t ∈ [0, 1], we have We further observe that It follows that, for all N ∈ N and all t ∈ [0, 1], , which is illustrated in Figure 3 and which establishes the claimed uniform boundedness.  Proof. According to Lemma 2.2, we have, for all t ∈ [0, 1], which implies that The desired result then follows by showing that, for ε > 0 fixed, the family is uniformly bounded, as illustrated in Figure 4. Employing a usual approach, we use the Dirichlet kernel, for N ∈ N, Integration by parts yields ds .
By the first mean value theorem for definite integrals, it follows that for t ∈ (0, π] fixed, there exists ξ ∈ [t, π], whereas for t ∈ [π, 2π) fixed, there exists ξ ∈ [π, t], such that Since cos N + 1 2 ξ is bounded above by one independently of ξ and as t 2 ∈ (0, π) for t ∈ (0, 2π) implies that 0 < sin t 2 ≤ 1, we conclude that, for all N ∈ N and for all t ∈ (0, 2π), which, for t ∈ [ε, 2π − ε], is uniformly bounded by 1/ sin ε 2 .  We can now prove the convergence in Theorem 1.2 on the diagonal away from the endpoints, which consists of a moment analysis to identify the moments of the limit function as well as an application of the Arzelà-Ascoli theorem to show that the limit function is continuous away from the endpoints. Alternatively, one could prove Corollary 1.3 directly with a similar approach as in the proof of Lemma 4.2, but integrating the Dirichlet kernel twice, and then deduce Theorem 1.2. However, as the moment analysis was already set up in Section 3 to determine the values of the Riemann zeta function at even positive integers, we demonstrate how to proceed with this approach.
By Lemma 4.1 and Lemma 4.2, the Arzelà-Ascoli theorem can be applied locally to any subsequence of {N C N 1 } N ∈N . Repeatedly using the Arzelà-Ascoli theorem and a diagonal argument, we deduce that there exists a subsequence of {N C N 1 } N ∈N which converges pointwise to a continuous limit function on the interval (0, 1). To prove that the full sequence converges pointwise and to identify the limit function, we proceed with the moment analysis initiated in Section 3. Applying Lemma 3.2, we see that, for m ∈ N 0 , The bound (4.1) together with For n ∈ N, we further have This shows that, for all n ∈ N 0 , If the sequence {N C N 1 } N ∈N failed to converge pointwise, we could use the Arzelà-Ascoli theorem and a diagonal argument to construct a second subsequence of {N C N 1 } N ∈N converging pointwise but to a different continuous limit function on (0, 1) compared to the first subsequence. Since this contradicts the convergence of moments, the claimed result follows.
We included the on-diagonal convergence in Theorem 1.2 as a separate statement to demonstrate that Corollary 1.3 is a consequence of Proposition 4.4, which is then used to prove the off-diagonal convergence in Theorem 1.2.
Proof of Corollary 1.3. Using the identity that, for k ∈ N, From (4.5) and Proposition 4.4, it follows that, for all t ∈ (0, 1), as claimed.
Since 0 < t − s < t + s < 2 for s, t ∈ (0, 1) with s < t, the convergence away from the diagonal is a consequence of Corollary 1.3.
Note that Theorem 1.2 states, for s, t ∈ [0, 1], if s = t and t ∈ (0, 1) 0 otherwise , which is the key ingredient for obtaining the characterisation of the limit fluctuations for the Karhunen-Loève expansion given in Theorem 1.1. We provide the full proof of Theorem 1.1 below after having determined the limit of 2N C N 2 as N → ∞.

4.2.
Fluctuations for the Fourier series expansion. Instead of setting up another moment analysis to study the pointwise limit of 2N C N 2 as N → ∞, we simplify the expression for C N 2 from Lemma 2.3 and deduce the desired pointwise limit from Corollary 1.3. cos(2kπt) Applying the identity (4.6) as well as the representation (2.1) and using the value for ζ(2) derived in Section 3, we have ∞ k=1 cos(2kπt) Once again exploiting the value for ζ(2), we obtain Using the expression for C N 2 from Lemma 2.3, it follows that, for s, t ∈ [0, 1], This implies that if t − s is an integer then, as a result of the limit (4.5), This can be summarised as, for s, t ∈ [0, 1], We finally prove Theorem 1.1 by considering characteristic functions. and 2N C N 2 , respectively. By the pointwise convergences (4.7) and (4.8) of the covariance functions in the limit N → ∞, for any n ∈ N and any t 1 , . . . , t n ∈ [0, 1], the characteristic functions of the Gaussian random vectors (F N,i t1 , . . . , F N,i tn ), for i ∈ {1, 2}, converge pointwise as N → ∞ to the characteristic function of the Gaussian random vector (F i t1 , . . . , F i tn ). Therefore, the claimed convergences in finite dimensional distributions are consequences of Lévy's continuity theorem.

Approximations of Brownian Lévy area
In this section, we consider approximations of second iterated integrals of Brownian motion, which is a classical problem in the numerical analysis of stochastic differential equations (SDEs), see [22]. Due to their presence within stochastic Taylor expansions, increments and second iterated integrals of multidimensional Brownian motion are required by high order strong methods for general SDEs, such as stochastic Taylor [22] and Runge-Kutta [33] methods. Currently, the only methodology for exactly generating the increment and second iterated integral, or equivalently the Lévy area, given by Definition 1.4, of a d-dimensional Brownian motion is limited to the case when d = 2. This algorithm for the exact generation of Brownian increments and Lévy area is detailed in [13]. The approach adapts Marsaglia's "rectangle-wedge-tail" algorithm to the joint density function of W 1 , A (1,2) 0,1 , which is expressible as an integral, but can only be evaluated numerically. Due to the subtle relationships between different entries in A 0,1 , it has not been extended to d > 2.
Obtaining good approximations of Brownian Lévy area in an L 2 (P) sense is known to be difficult. For example, it was shown in [8] that any approximation of Lévy area which is measurable with respect to N Gaussian random variables, obtained from linear functionals of the Brownian path, cannot achieve strong convergence faster than O(N − 1 2 ). In particular, this result extends the classical theorem of Clark and Cameron [4] which establishes a best convergence rate of O(N − 1 2 ) for approximations of Lévy area based on only the Brownian increments {W (n+1)h − W nh } 0≤n≤N −1 . Therefore, approximations have been developed which fall outside of this paradigm, see [5,11,32,35]. In the analysis of these methodologies, the Lévy area of Brownian motion and its approximation are probabilistically coupled in such a way that L 2 (P) convergence rates of O(N −1 ) can be established.
We are interested in the approximations of Brownian Lévy area that can be obtained directly from the Fourier series expansion (1.3) and the polynomial expansion (1.5) of the Brownian bridge. For the remainder of the section, the Brownian motion (W t ) t∈[0,1] is assumed to be d-dimensional and (B t ) t∈[0,1] is its associated Brownian bridge.
We first recall the standard Fourier approach to the strong approximation of Brownian Lévy area.
Theorem 5.1 (Approximation of Brownian Lévy area via Fourier coefficients, see [22, p. 205] and [31, p. 99]). For n ∈ N, we define a random antisymmetric d×d matrix A n by, for i, j ∈ {1, . . . , d}, where the normal random vectors {a k } k∈N0 and {b k } k∈N are the coefficients from the Brownian bridge expansion (1.3), that is, the coordinates of each random vector are independent and defined according to (1.4). Then, for i, j ∈ {1, . . . , d} with i = j, we have Remark 5.2. Using the covariance structure given by (2.5), (2.6), (2.7) and the independence of the components of a Brownian bridge, it immediately follows that the coefficients {a k } k∈N0 and {b k } k∈N are jointly normal with a 0 ∼ N 0, 1 3 I d , a k , b k ∼ N 0, 1 2k 2 π 2 I d , cov(a 0 , a k ) = − 1 k 2 π 2 I d and cov(a l , b k ) = 0 for k ∈ N and l ∈ N 0 .
In practice, the above approximation may involve generating the N independent random vectors {a k } 1≤k≤N followed by the coefficient a 0 , which will not be independent, but can be expressed as a linear combination of {a k } 1≤k≤N along with an additional independent normal random vector. Without this additional normal random vector, we obtain the following discretisation of Lévy area.  [23,30,35]). For n ∈ N, we define a random antisymmetric d × d matrix A n by, for i, j ∈ {1, . . . , d}, where the sequences {a k } k∈N and {b k } k∈N of independent normal random vectors are the same as before. Then, for i, j ∈ {1, . . . , d} with i = j, we have Proof. As for Theorem 5.1, the result follows by direct calculation. The constant is larger because, for i ∈ {1, . . . , d} and k ∈ N, which yields the required result.
Finally, we give the approximation of Lévy area corresponding to the polynomial expansion (1.5). Although this series expansion of Brownian Lévy area was first proposed in [24], a straightforward derivation based on the polynomial expansion (1.5) was only established much later in [25]. However in [24,25], the optimal bound for the mean squared error of the approximation is not identified. We will present a similar derivation to [25], but with a simple formula for the mean squared error.
Theorem 5.4 (Polynomial approximation of Brownian Lévy area, see [24, p. 47] and [25]). For n ∈ N 0 , we define a random antisymmetric d × d matrix s A n by, for n ∈ N and i, j ∈ {1, . . . , d}, where the normal random vectors {c k } k∈N are the coefficients from the polynomial expansion (1.5), that is, the coordinates are independent and defined according to (1.6), and we set s A (i,j) 0 Then, for n ∈ N 0 and for i, j ∈ {1, . . . , d} with i = j, we have Remark 5.5. By applying Lemma 2.1, the orthogonality of shifted Legendre polynomials and the independence of the components of a Brownian bridge, we see that the coefficients {c k } k∈N are independent and distributed as c k ∼ N 0, 1 2k+1 I d for k ∈ N. Proof. It follows from the polynomial expansion (1.5) that, for i, j ∈ {1, . . . , d} with i = j, where the series converge in L 2 (P). To simplify (5.1), we use the identities in (2.11) for shifted Legendre polynomials as well as the orthogonality of shifted Legendre polynomials to obtain that, for k, l ∈ N, Evaluating the above integrals gives, for k, l ∈ N, In particular, for k, l ∈ N, this implies that Therefore, by the bounded convergence theorem in L 2 (P), we can simplify the expansion (5.1) to where, just as before, the series converges in L 2 (P).
where the second line follows by integration by parts. As and Q 1 (t) = 2t − 1, the above and (5.3) imply that, for i, j ∈ {1, . . . , d}, By the independence of the normal random vectors in the sequence {c k } k∈N , it is straightforward to compute the mean squared error in approximating A 0,1 and we obtain, for n ∈ N and for i, j ∈ {1, . . . , d} with i = j, In particular, the polynomial approximation of Brownian Lévy area is more accurate than the Kloeden-Platen-Wright approximation, both of which use only independent Gaussian vectors.
Remark 5.7. It was shown in [8] that 1 π 2 1 N is the optimal asymptotic rate of mean squared convergence for Lévy area approximations that are measurable with respect to N Gaussian random variables, obtained from linear functionals of the Brownian path.
As one would expect, all the Lévy area approximations converge in L 2 (P) with a rate of O(N − 1 2 ) and thus the main difference between their respective accuracies is in the leading error constant. More concretely, for sufficiently large N , the approximation based on the Fourier expansion of the Brownian bridge is roughly 11% more accurate in L 2 (P) than that of the polynomial approximation. On the other hand, the polynomial approximation is easier to implement in practice as all of the required coefficients are independent. Since it has the largest asymptotic error constant, the Kloeden-Platen-Wright approach gives the least accurate approximation for Brownian Lévy area.
We observe that the leading error constants for the Lévy area approximations resulting from the Fourier series and the polynomial expansion coincide with the average L 2 (P) error of their respective fluctuation processes, that is, applying Fubini's theorem followed by the limit theorems for the fluctuation processes which, for instance, for the polynomial approximation follows directly from (5.2) and Remark 5.5, then in terms of the fluctuation processes (F N t ) t∈[0,1] defined by the error of the Lévy area approximation can be expressed as Thus, by Itô's isometry and Fubini's theorem, the leading error constant in the mean squared error is indeed given by This connection could be interpreted as an asymptotic Itô isometry for Lévy area approximations.  Fourier series (Kloeden-Platen [22] and Milstein [31]) Fourier series (Kloeden-Platen-Wright [23] and Milstein [30]) r and Q k denoting the shifted Legendre polynomial of degree k Table 2. Table summarising the Lévy area expansions considered in this paper.