# Accurate and efficient structure factors in
ultrasoft pseudopotential and projector augmented wave DFT

Benjamin X. Shi<EMAIL_ADDRESS>Rebecca J. Nicholls  Jonathan R. Yates
Department of Materials, University of Oxford, Parks Road, Oxford, OX1 3PH,
United Kingdom

(August 28, 2024)

###### Abstract

Structure factors obtained from diffraction experiments are one of the most
important quantities for characterizing the electronic and structural
properties of materials. Methods for calculating this quantity from plane-wave
density functional theory (DFT) codes are typically prohibitively expensive to
perform, requiring the electron density to be constructed and evaluated on
dense real-space grids. Making use of the projector functions found in both
the Vanderbilt ultrasoft pseudopotential and projector augmented wave methods,
we implement an approach to calculate structure factors which avoids the use
of a dense grid by separating the rapidly changing contributions to the
electron density and treating them on logarithmic radial grids. This approach
is successfully validated against structure factors obtained from all-electron
DFT and experiments for three prototype systems, allowing structure factors to
be obtained at all-electron accuracy at a fraction of the cost of previous
approaches for plane-wave DFT.

††preprint: APS/123-QED

##  I  Introduction

The structure factor (SF) is a fundamentally important quantity in the
physical and biological sciences. Experimental SFs have become the principal
tool for characterizing the geometrical structure (e.g., atomic positions and
elements) of a range of crystalline matter, from simple materials [1] to
complex bio-molecules [2]. Beyond structural information, the (x-ray) SFs are
the Fourier coefficients of the electron density (ED) - $n(\mathbf{r})$ \- and
allow for its reconstruction through an inverse Fourier transform [3]. The ED
contains a wealth of information, as established by the Hohenberg-Kohn
theorems [4], which states that the ground-state electronic properties of a
system are a unique functional of the ground-state ED. Reconstructed EDs have
been used to investigate the properties of a range of materials, from the
anisotropic elastic constants of Al [5] to the electronic origins of high-
temperature cuprate superconductors [6, 7].

Accurate SFs calculated from first-principles methods, particularly density
functional theory (DFT) are necessary for the above applications of SFs. For
example, SFs computed from DFT have been used to augment experimental SFs to
allow full reconstruction of the ED. By themselves, SFs obtained from x-ray
[8, 7, 9], $\gamma$-ray [10, 11, 12] or electron diffraction [13, 14, 15, 16]
experiments can only provide a finite set of SFs, introducing Fourier series
truncation errors to the reconstructed ED [17] if used alone. Many of these
diffraction experiments are also incapable of obtaining the phase of the
complex SFs in non-centrosymmetric crystals – the phase problem [18]. DFT
helps to alleviate these problems as it can generate many SFs to augment those
that are not provided by experiments while also providing phase information.

Accurate SFs from DFT are also useful for assessing the quality of
experimental diffraction techniques, which can suffer from problems afflicting
their accuracy. For example, extinction effects [19, 20] and the source of
x-rays [21] can affect measurements in x-ray diffraction experiments. In
electron diffraction, material preparation [22] and instrument distortions
[23] can serve as potential error sources.

High-precision SFs from experiments can also help validate the approximations
used to make DFT calculations computationally tractable. Based upon the
Hohenberg-Kohn theorems, practical DFT calculations attempt to approximate the
exact energy functional, which is unknown, to the ED through density
functional approximations (DFAs). There is a whole ‘zoo’ [24] of available
DFAs, with no systematic manner to determine their accuracy. In recent years,
there has been evidence [25] to show that although modern DFAs give improved
energetic descriptions of (atomic) systems, the description of the ED is
worse. This deficiency arises because many modern DFAs have been constructed
through empirical fitting of reference energies, typically neglecting the ED
due to the lack of available references. Towards this end, high-precision SFs
obtained from diffraction experiments can be compared to those obtained from
various DFAs to assess their relative accuracy, and it has been successfully
performed for several materials [26, 27, 14, 28].

Beyond the need for accurate SFs from DFT, it is also highly desirable that
they can be obtained efficiently without incurring heavy computational
burden/time to allow for more complex systems to be tackled. Unfortunately,
current approaches to calculate the SF from DFT methods are not efficient,
requiring a high computational cost. There are an array of DFT methods, but
the key methods used for computing the SF are either all-electron (AE) or
plane-wave pseudopotential (PP) DFT methods. AE DFT approaches treat all the
electrons in the system explicitly, enabling highly accurate calculations.
This approach is the most common approach for computing SFs due to the ease at
obtaining the SF from the ED. However, it comes at significant computational
cost due to the $\sim O(N^{3})$ [29] scaling of DFT.

In PP DFT, the inclusion of a pseudopotential [30] means that only valence
electrons need to be treated, decreasing the computational cost significantly.
The resulting ED from the PP DFT SCF calculation cannot be directly used to
obtain the SF (as in AE DFT) as it includes only valence electrons and has
been smoothed – pseudized – near the core. In modern plane-wave DFT codes, the
Vanderbilt ultrasoft pseudopotential [31] or the projector augmented wave
(PAW) [32, 33] methods can be employed to restore the AE total ED from the
pseudized (PS) valence ED by adding a compensating augmentation charge and
including the (frozen-)core ED. In the typical plane-wave DFT codes used, the
AE total ED must be reconstructed on a real-space fast Fourier transform (FFT)
grid. This FFT grid has to be several orders of magnitude denser than the
default used for the PS valence ED to have sufficient spatial resolution to
accurately capture the rapid oscillations of the AE total ED near the nucleus.
As a result, the calculations require a large amount of memory and time, with
several studies explicitly highlighting the difficulty with converging the
total ED to a sufficient precision due to computational limitations [28, 34].

In this paper, we propose a highly efficient and accurate method of
calculating the AE SF from PP DFT that can be implemented in both the
Vanderbilt ultrasoft pseudopotential and PAW methods. It works by removing the
rapidly changing contributions (i.e. core and augmentation charges) from the
AE total ED and treating their contributions to the AE SF separately. With
this change, the FFT grid does not need to be increased beyond its default
size since it only represents the PS valence ED. The core and augmentation
contributions to the ED, which exist only within a small region around the
atom, will have their radial components treated on atom-centered logarithmic
radial support grids, while the angular components are treated analytically.
As the radial support grids are one-dimensional, they can be made very dense
near the origin where the ED contributions vary the most, allowing for high
accuracy at a low computational cost. This method is successfully validated
against AE DFT and experiment for a range of materials.

##  II  The AE Electron Density

This section details how the AE ED is obtained from PP DFT calculations using
the Vanderbilt ultrasoft pseudopotential and PAW methods. The theory is
applicable for both methods, so unless otherwise stated, PAW DFT will be used
to denote both methods hereafter.

In PAW DFT, the AE wave function $\psi_{n}(\mathbf{r})$ for each of the $n$
valence (Kohn-Sham) orbitals can be reconstructed from its corresponding
pseudized (PS) wave function $\tilde{\psi}_{n}(\mathbf{r})$ through a linear
transformation [32]:

$\psi_{n}(\mathbf{r})=\tilde{\psi}_{n}(\mathbf{r})+\sum_{\mathbf{R}ju}\left[\phi^{j}_{u}(\mathbf{r}_{\mathbf{R}}^{j})-\tilde{\phi}^{j}_{u}(\mathbf{r}_{\mathbf{R}}^{j})\right]\braket{\tilde{p}^{j}_{u}}{\tilde{\psi}_{n}}.$
(1)

Within this expression, $\phi^{j}_{u}(\mathbf{r}_{\mathbf{R}}^{j})$ and
$\tilde{\phi}^{j}_{u}(\mathbf{r}_{\mathbf{R}}^{j})$ are the AE and PS partial
waves respectively for each atom $j$ in the unit cell, with the projectors
$\bra{\tilde{p}^{j}_{u}}$ designed to be dual to the PS partial waves:
$\braket{\tilde{p}^{j}_{u}}{\tilde{\phi}^{j^{\prime}}_{u^{\prime}}}=\delta_{jj^{\prime}}\delta_{uu^{\prime}}$.
The AE partial waves are a set of wave functions obtained from the
corresponding reference atom [35], where $u$ is the composite index for the
angular momentum quantum numbers $l,m$ as well as an index $k$ to label
partial waves constructed at different reference energies [32]. The vectors

$\mathbf{r}_{\mathbf{R}}^{j}=\mathbf{r}-\mathbf{r}_{j}-\mathbf{R}$ (2)

are used to denote the spatial dependence of the partial waves and projectors
to emphasize that these functions are atom-centered. These equations assume a
periodic material, where $\mathbf{R}$ are the (infinite) set of lattice
vectors and $\mathbf{r}_{j}$ denotes the position of atom $j$ in the unit
cell. Computationally, the partial waves and projectors can be expressed as a
radial function (stored on logarithmic radial grids) multiplied by a spherical
harmonic:

$\phi_{u}(\mathbf{r})=R_{lk}(r)Y_{lm}(\hat{\mathbf{r}}),$ (3)

with the PS partial waves differing from the AE partial waves only within a
cutoff $r_{c}^{l}$, where its radial component has been pseudized.

The AE valence ED $n_{\text{val}}(\mathbf{r})$ can be given as a sum of two
contributions:

$n_{\textrm{val}}(\mathbf{r})=\tilde{n}_{\textrm{val}}(\mathbf{r})+n_{\textrm{aug}}(\mathbf{r}).$
(4)

The first term is the PS valence ED, which is the ED resulting from the PS
(valence) wave functions in Eq. 1:

$\tilde{n}_{\textrm{val}}(\mathbf{r})=\sum_{n}f_{n}|\tilde{\psi}_{n}(\mathbf{r})|^{2},$
(5)

where $f_{n}$ are the occupation numbers. The FFT grid is designed to store
this smooth function, with its default size sufficient to sample and represent
it fully. The second term is the augmentation charge, which restores the PS
valence ED to the AE valence ED, taking the form:

$n_{\text{aug}}(\mathbf{r})=\sum_{\mathbf{R}ju_{1}u_{2}}\rho_{u_{1}u_{2}}^{j}Q_{u_{1}u_{2}}^{j}(\mathbf{r}_{\mathbf{R}}^{j}).$
(6)

The augmentation functions $Q_{u_{1}u_{2}}^{j}(\mathbf{r}_{\mathbf{R}}^{j})$
are defined as

$Q_{u_{1}u_{2}}^{j}(\mathbf{r})=\phi^{j}_{u_{1}}(\mathbf{r})^{*}\phi^{j}_{u_{2}}(\mathbf{r})-\tilde{\phi}^{j}_{u_{1}}(\mathbf{r})^{*}\tilde{\phi}^{j}_{u_{2}}(\mathbf{r}),$
(7)

where $u_{1}$ and $u_{2}$ are two sets of $u$ indices, with
$\rho_{u_{1}u_{2}}^{j}$ giving the occupancy of each $u_{1},u_{2}$
augmentation function channel for atom $j$:

$\rho_{u_{1}u_{2}}^{j}=\sum_{n}f_{n}\braket{\tilde{\psi}_{n}}{\tilde{p}_{u_{1}}^{j}}\braket{\tilde{p}_{u_{2}}^{j}}{\tilde{\psi}_{n}}.$
(8)

The augmentation function $Q_{u_{1}u_{2}}^{j}(\mathbf{r})$ will be localized
around atom $j$ in an ‘augmentation’ region. The two most common types of
pseudopotentials are the norm-conserving PPs (NCPs) and ultrasoft PPs (USPs).
NCPs are constructed to preserve the norm of the wave (i.e. the integral of
the augmentation functions within the augmentation region is zero by
construction), while USPs relax this condition, requiring fewer plane-waves to
describe the PS valence wave functions at the cost of additional complexity.
It is standard practice when using NCPs to neglect any augmentation to the
charge density during the calculation of the ground-state. Following from Eqs.
3 and 7, $Q_{u_{1}u_{2}}^{j}(\mathbf{r})$ can be expressed as the product of a
radial function $\Delta R_{u_{1}u_{2}}^{j}(r)$ and two spherical harmonics:

$Q_{u_{1}u_{2}}^{j}(\mathbf{r})=\Delta
R_{u_{1}u_{2}}^{j}(r)Y^{*}_{l_{1}m_{1}}(\hat{\mathbf{r}})Y_{l_{2}m_{2}}(\hat{\mathbf{r}}).$
(9)

Within the PAW method, the radial functions in the augmentation functions are
treated on atom-centered logarithmic radial support grids [33], with the
spherical harmonics treated analytically, while in the Vanderbilt ultrasoft
pseudopotential approach, these functions are typically pseudized and placed
onto the FFT grid [36]. While this pseudization gives accurate total energies
[37], the use of a PS compared to the AE augmentation charge introduces
significant errors in the calculations of SF, as will be shown in Sec. IV.

The PAW method typically utilizes the frozen core approximation. Under this
scheme, the core ED is a superposition of the core EDs
$\rho_{\textrm{core}}^{j}$ obtained from isolated atoms:

$n_{\textrm{core}}(\mathbf{r})=\sum_{\mathbf{R}j}\rho^{j}_{\textrm{core}}(|\mathbf{r}_{\mathbf{R}}^{j}|).$
(10)

These atomic core EDs are spherically symmetric and obtained from solving the
radial Kohn-Sham Schrödinger equation [35] on the logarithmic radial support
grids.

Overall, the total AE ED will then be the sum of the PS valence, augmentation
charge and frozen-core EDs:

$n(\mathbf{r})=\tilde{n}_{\textrm{val}}(\mathbf{r})+n_{\textrm{aug}}(\mathbf{r})+n_{\textrm{core}}(\mathbf{r}).$
(11)

##  III  The AE Structure Factors

This section details how the AE SFs can be obtained from the AE ED in PAW DFT
codes. The AE SFs $F(\mathbf{H})$ are the Fourier coefficients of the AE
(total) ED:

$F(\mathbf{H})=\mathcal{F}[n(\mathbf{r})],$ (12)

where $\mathbf{H}=h\mathbf{a}^{*}+k\mathbf{b}^{*}+l\mathbf{c}^{*}$ is the
scattering vector corresponding to the $(hkl)$ plane and the crystallographic
convention for the primitive reciprocal lattice vectors
$\mathbf{a}^{*},\mathbf{b}^{*}$ and $\mathbf{c}^{*}$ has been used. Current
methods of obtaining the AE SFs from the AE ED in PAW DFT codes involves
constructing the AE ED onto a uniform grid before applying a fast Fourier
transform (FFT). This approach is highly inefficient because the FFT grid used
has to be several orders of magnitude denser than the default (designed for
the PS valence ED only) to accommodate the rapidly varying augmentation and
core charges near the nuclei. As shown in Fig. S1 of the Supplemental Material
[38], this leads to orders of magnitude increases in time and peak RAM. The
impact of the resulting additional computational burden has been noted in
previous studies [28].

In this section, we propose and derive a new approach to calculate AE SFs
efficiently and accurately in the PAW DFT method. It works by separating the
three individual contributions to the AE total ED from PAW DFT and treating
the two terms that require high spatial resolution near the nuclei:
$n_{\textrm{aug}}$ and $n_{\textrm{core}}$, on logarithmic radial support
grids and analytically for radial and angular components respectively. These
radial grids are one-dimensional with a high density of points placed near the
nuclei to achieve high accuracy efficiently.

From the linearity principle of the Fourier transform, the SF can be separated
into three contributions:

$F(\mathbf{H})=\tilde{F}_{\textrm{val}}(\mathbf{H})+F_{\textrm{core}}(\mathbf{H})+F_{\textrm{aug}}(\mathbf{H}),$
(13)

where $\tilde{F}_{\textrm{val}}(\mathbf{H}),F_{\textrm{core}}(\mathbf{H})$ and
$F_{\textrm{aug}}(\mathbf{H})$ are the Fourier transforms of
$\tilde{n}_{\textrm{val}}(\mathbf{r}),n_{\textrm{core}}(\mathbf{r})$ and
$n_{\textrm{aug}}(\mathbf{r})$ respectively. By nature,
$\tilde{n}_{\textrm{val}}(\mathbf{r})$ is constructed to be fully described on
the coarse (default) FFT grid size. Thus, its Fourier transform can be
computed efficiently using the FFT method, a fundamental component of all
plane-wave DFT codes.

The theory for treating the core and augmentation contributions to the SF in
the next two subsections relies on some of the methods developed for the
independent atom model (IAM) – commonly used in crystallography. The IAM is
constructed as a summation of isolated atomic densities $\rho^{j}(\mathbf{r})$
about their atomic positions:

$n_{\textrm{IAM}}(\mathbf{r})=\sum_{\mathbf{R}j}\rho^{j}(\mathbf{r}_{\mathbf{R}}^{j}),$
(14)

and it is possible to show that its SF takes the form [3]:

$F_{\textrm{IAM}}(\mathbf{H})=\sum_{j}f^{j}(\mathbf{H})\exp(i2\pi\mathbf{H}\cdot\mathbf{r}_{j}),$
(15)

where $f^{j}(\mathbf{H})$ is the atomic scattering factor, defined to be the
Fourier transform of the corresponding $\rho^{j}(\mathbf{r})$:

$f^{j}(\mathbf{H})=\int\rho^{j}(\mathbf{r})\exp(i2\pi\mathbf{H}\cdot\mathbf{r})\mathop{}\\!\mathrm{d}{\mathbf{r}}.$
(16)

###  III.1 Core contribution

Recalling Eq. 10, the core ED is a summation of atom-centered densities, much
like the IAM. Thus, using the SF expressions derived for the IAM, the SF
contribution from the core ED takes the form:

$F_{\textrm{core}}(\mathbf{H})=\sum_{j}f^{j}_{\textrm{core}}(\mathbf{H})\exp(i2\pi\mathbf{H}\cdot\mathbf{r}_{j}),$
(17)

where the atomic core scattering factor $f^{j}_{\textrm{core}}(\mathbf{H})$ is
the Fourier transform of the corresponding
$\rho^{j}_{\textrm{core}}(\mathbf{r})$:

$f_{\textrm{core}}^{j}(\mathbf{H})=\int\rho^{j}_{\textrm{core}}(|\mathbf{r}|)\exp(i2\pi\mathbf{H}\cdot\mathbf{r})\mathop{}\\!\mathrm{d}{\mathbf{r}}.$
(18)

As $\rho_{\textrm{core}}(|\mathbf{r}|)$ is a spherically symmetric function
about the origin, it is more convenient to use the spherical polar coordinate
system. Within this coordinate system, the plane-wave function can be expanded
into complex spherical harmonic functions [39]:

$\exp(i2\pi\mathbf{H}\cdot\mathbf{r})=4\pi\sum_{l=0}^{\infty}\sum_{m=-l}^{l}i^{l}j_{l}(2\pi
Hr)Y_{l,m}(\hat{\mathbf{H}})Y_{l,m}^{*}(\hat{\mathbf{r}}),$ (19)

where $j_{l}(r)$ is the spherical Bessel function of order $l$. Substituting
Eq. 19 into Eq. 18 would then yield:

$f_{\textrm{core}}^{j}(H)=\int\rho_{\textrm{core}}^{j}(r)4\pi r^{2}j_{0}(2\pi
Hr)\mathop{}\\!\mathrm{d}r,$ (20)

as only the $l=0$ terms persist due to the rotational invariance of
$\rho_{\textrm{core}}^{j}(r)$ (and in turn $f_{\textrm{core}}^{j}(H)$). This
equation can be evaluated on the native logarithmic radial grids for each atom
$j$ within the unit cell and as a result, the core contribution to the SF no
longer needs to be evaluated on the FFT grid.

###  III.2 Augmentation contribution

Like the core ED and IAM, the augmentation charge also consists of a summation
of atom-centered functions, so it can be written to take the form:

$F_{\textrm{aug}}(\mathbf{H})=\sum_{j}\exp(i2\pi\mathbf{H}\cdot\mathbf{r}_{j})\sum_{u_{1}u_{2}}f^{j}_{u_{1}u_{2}}(\mathbf{H}),$
(21)

where

$f^{j}_{u_{1}u_{2}}(\mathbf{H})=\rho_{u_{1}u_{2}}^{j}\int
Q_{u_{1}u_{2}}^{j}(\mathbf{r})\exp(i2\pi\mathbf{H}\cdot\mathbf{r})\mathop{}\\!\mathrm{d}\mathbf{r}.$
(22)

The augmentation charge differs from the core ED in that these atom-centered
functions $Q_{u_{1}u2}^{j}(\mathbf{r})$ are not spherically symmetric about
the origin, instead being expressed as a multiple of a radial function and two
spherical harmonics (see Eq. 9). The spherical harmonics form a complete set,
so their products can be expressed as an expansion of single spherical
harmonics:

$Y_{l_{1}m_{1}}^{\ast}(\hat{\mathbf{r}})Y_{l_{2}m_{2}}(\hat{\mathbf{r}})=\sum_{L}C_{l_{1}l_{2}L}^{m_{1}m_{2}}Y_{LM}(\hat{\mathbf{r}}),$
(23)

where $M=-m_{1}+m_{2}$, $|l_{1}-l_{2}|\leq L\leq l_{1}+l_{2}$ and
$C_{l_{1}l_{2}L}^{m_{1}m_{2}}$ are the Clebsch-Gordan coefficients. The
augmentation functions $Q_{u_{1}u2}^{j}(\mathbf{r})$ can then be rewritten as

$Q_{u_{1}u2}^{j}(\mathbf{r})=\Delta
R_{u_{1}u_{2}}^{j}(r)\sum_{L}C_{l_{1}l_{2}L}^{m_{1}m_{2}}Y_{LM}(\hat{\mathbf{r}}).$
(24)

When this new formulation of $Q_{u_{1}u_{2}}^{j}(\mathbf{r})$ is substituted
into Eq. 22, $f^{j}_{u_{1}u_{2}}(\mathbf{H})$ can be further simplified to

$\begin{split}f_{u_{1}u_{2}}^{j}(\mathbf{H})=\rho_{u_{1}u_{2}}^{j}\sum_{L}C_{l_{1}l_{2}L}^{m_{1}m_{2}}Y_{LM}(\hat{\mathbf{H}})4\pi
i^{L}\\\ \times\int j_{L}(2\pi Hr)\Delta
R_{u_{1}u_{2}}^{j}(r)r^{2}\mathop{}\\!\mathrm{d}r,\end{split}$ (25)

where the expansion of the plane-wave in spherical harmonics (Eq. 19) and
subsequently the orthonormal properties of complex spherical harmonics has
been used to arrive at the final expression. Within Eq. 25, the $\int
j_{L}(2\pi Hr)\Delta R_{u_{1}u_{2}}^{j}(r)r^{2}\mathop{}\\!\mathrm{d}r$
integral can be evaluated to high precision on the logarithmic radial support
grids, while the angular components are treated analytically using spherical
harmonics, removing the need to evaluate the augmentation contribution on the
FFT grid.

With the method outlined in this section, only the default FFT grid size has
to be used to evaluate converged AE SFs in PAW DFT since it only houses the PS
valence ED, with the other terms being evaluated efficiently on radial support
grids, resulting in orders of magnitude savings in computational cost (see
Fig. S1 of the Supplemental Material [38]).

##  IV  Validation against AE DFT

To ensure that the approach outlined in Sec. III is accurate, we have compared
the SFs produced via this approach in PAW DFT against AE DFT calculations for
the same DFA for three prototypical systems.

###  IV.1 Computational details

The approach was implemented in CASTEP [40], which is a plane-wave PP DFT code
that can use either NCP or USP for ground-state calculation and the PAW method
as a post-processing approach to calculate properties such as hyperfine
tensors, NMR properties [41] and now SFs. We refer to AE SFs obtained from
using this full approach (with augmentation charges treated on radial support
grids) as AE-USP or AE-NCP depending on the class of PP used. To investigate
the importance of using the full valence ED, results will also be presented
for what we shall refer to as the PS-USP and PS-NCP approaches. The PS-NCP
approach uses the pseudized valence ED and the ED from the core electrons,
without any augmentation charge. The PS-USP also includes the contribution
from the pseudized augmentation charge evaluated on the same FFT grid as the
valence ED. Details of the pseudopotentials used are given in Sec. II of the
Supplemental Material [38], with all input parameters and output files from
this work made available online (see data availability statement at end of
paper).

The AE DFT SF calculations are performed using the WIEN2k code [42], which
utilizes a APW+lo basis set. For both AE and PAW DFT calculations, the energy
cutoff ($R_{\textrm{MT}}K$ for AE DFT) and number of $k$-points were varied to
ensure that the SFs were numerically converged to 0.001 (e)lectrons, the
limits of experimental accuracy. SFs were compared for both the LDA and PBE
DFAs in three materials that highlight different bonding and unit cell
systems: diamond Si (covalent), hcp Mg (metallic) and rocksalt MgO (ionic). In
each system, the amplitudes of SFs with scattering length
$\mathbf{s}=|\mathbf{H}|/2$ below 1.5 Å-1 were analyzed, removing any SFs that
have equivalent amplitudes due to symmetry or under Friedel’s law [43]. All
forbidden reflections, under the spherical IAM, were also removed, except for
the (222) reflection in Si as it has a noticeable amplitude due to the
asphericity of the atoms [44]. The $R$-factor is used to assess the accuracy
of our calculations relative to a reference and is defined as:

$R^{\textrm{ref.}}=\dfrac{\sum_{\mathbf{H}}|F^{\textrm{ref.}}(\mathbf{H})-F^{\textrm{calc.}}(\mathbf{H})|}{\sum_{\mathbf{H}}|F^{\textrm{ref.}}(\mathbf{H})|},$
(26)

where the reference will be either WIEN2k calculations in this section or
experimental values in the following section.

###  IV.2 Results and discussion

Figure 1:  The difference between APW+lo WIEN2k SF calculations and those
obtained from CASTEP as a function of scattering vector $s$ for (a) Si, (b) Mg
and (c) MgO. In CASTEP, PAW DFT calculations were performed employing both
USPs and NCPs, with each either using an all-electron (AE) or pseudized (PS)
augmentation charge. The PBE functional was used.

In Fig. 1, the difference between AE DFT WIEN2k and PAW DFT CASTEP
calculations, using AE or PS augmentation charges for both USPs and NCPs, were
performed for scattering vectors up to 1.5 Å-1 for all three materials
systems. These plots were for the PBE DFA as we found that the LDA mirrored
the plots (see Fig. S2 of the Supplemental Material [38]), suggesting that the
differences arise due to the particular pseudization schemes for the
pseudopotentials and augmentation charge, which is independent of DFA. For the
two monoatomic systems, we see that in the $s$ range of 0.25-1.25 Å-1, there
is a systematic underestimation of the SFs for PAW DFT methods that utilize a
pseudized (PS) augmentation charge (e.g., PS-NCP and PS-USP). This
underestimation arises because these high Fourier components are removed when
the oscillating valence ED is made smooth. The magnitude of this
underestimation depends on the exact pseudization schemes employed for the
pseudopotentials and augmentation charge (for USP) around each atom. Beyond
1.1 Å-1, the SF contributions from the PS valence ED and PS augmentation
charge become less than 0.01 e, so that the SF behavior becomes a reflection
of the core ED used by the pseudopotential.

Table 1:  $R^{\textrm{WIEN2k}}$ (%) values for the PAW DFT method employing either the USP and NCP in combination with either AE or PS augmentation charges. |  | AE-USP | PS-USP | AE-NCP | PS-NCP
---|---|---|---|---|---
Si | LDA | 0.05 | 0.56 | 0.04 | 1.03
PBE | 0.04 | 0.57 | 0.04 | 0.99
Mg | LDA | 0.04 | 3.75 | 0.02 | 0.29
PBE | 0.03 | 3.66 | 0.02 | 0.28
MgO | LDA | 0.05 | 3.26 | 0.08 | 2.01
PBE | 0.05 | 3.19 | 0.08 | 1.98

Regardless of the pseudopotential, the deviations resulting from a PS
augmentation charge is completely removed when the AE augmentation charge is
used. The $R^{\textrm{WIEN2k}}$ (in Table 1) lowers by at least an order of
magnitude for all studied systems from a $R^{\textrm{WIEN2k}}$ of 0.28-3.75 %
for PS-USP and PS-NCP to a $R^{\textrm{WIEN2k}}$ of 0.02-0.08 % for AE-USP and
AE-NCP. The errors of AE-USP and AE-NCP with respect to the APW+lo WIEN2k SF
values are small enough that they are of the same order as the experimental
errors that can arise in high quality x-ray diffraction experiments. Thus, SFs
from AE-USP and AE-NCP (but not PS-USP or PS-NCP) can be used to compare
against experiments, as we have done in Sec. V.

MgO is the only system studied here which contains more than one type of
atomic species. For the PS-NCP calculation, there are two observed behaviors
in the SFs, with the $(hkl)$ reflections where $h,k,l$ are all odd
overestimating the SFs, while those with $h,k,l$ all even underestimating the
SFs. This contrasting behavior arises because these two sets arise from
different types of reflections. Using Eq. 20 for an IAM, the structure factors
where $h,k,l$ are all odd take the form: $4(f_{\textrm{Mg}}-f_{\textrm{O}})$,
arising from a scattering difference between its ions [3] while those that are
all even arises from a constructive summation:
$4(f_{\textrm{Mg}}+f_{\textrm{O}})$. Thus, the SFs with odd $h,k,l$ become
sensitive to the specific pseudopotentials of the two atoms, and may
overestimate if one pseudopotential underestimates its scattering factor more
strongly than the other.

##  V Comparison against experiment

Ultimately, the aim of computing AE SFs from DFT is to use it together with
experimental SFs, either to augment the (limited) experimental SFs for better
ED reconstructions or to validate the approximations in both experiments and
DFT. In this section, we will showcase the strength of our approach for the
latter application of AE SFs by comparing DFT with two different DFAs and
pseudopotentials against experiment, revealing key insights into the effects
of these approximations on the accuracy of the ED.

###  V.1 Computational details

Figure 2:  (a) The difference between experimental atomic scattering factors
and those obtained from AE-USP CASTEP as a function of scattering vector $s$
for Si for both the LDA and PBE density functional approximations (DFAs). (b)
The total difference between the atomic scattering factors of the PBE and LDA
DFAs and the separation of this difference into its core and valence
(including augmentation) contributions.

X-ray SFs obtained from diffraction experiments are influenced by the thermal
vibrations of the atoms within the crystal [43], tending to reduce the
intensity of the diffracted beams. Comparison of (static) SFs obtained from
DFT requires incorporation of these thermal effects to allow for direct
comparison, or alternatively, removal of the thermal effects from experimental
numbers.

To a good approximation, the SF is the Fourier transform of the thermally
averaged ED $\langle n(\mathbf{r})\rangle$. Within the IAM, thermal effects
can be incorporated into the static SF by multiplying each atomic scattering
factor $f_{j}(\mathbf{H})$ in Eq. 16 by an isotropic harmonic temperature
factor $T_{j}(\mathbf{H})$:

$F(\mathbf{H})=\sum_{j}\exp(i2\pi\mathbf{H}\cdot\mathbf{r}_{j})f_{j}(\mathbf{H})T_{j}(\mathbf{H}),$
(27)

where

$T_{j}(\mathbf{H})=\exp(-B_{j}|\mathbf{H}|^{2}/4).$ (28)

The Debye-Waller (DW) factor $B_{j}$ for atom $j$ can be obtained from either
fitting of the static SF to experimental SFs [45, 46] or using ab-initio
phonon dispersion curves [47]. The derivation of the temperature factor makes
two key assumptions: (i) the nuclei vibrate isotropically about their
equilibrium positions and (ii) the atomic densities follow the nuclear motion
perfectly. The second assumption requires that the crystal ED can be divided
into a summation of atom-centered density terms.

Within PAW DFT, the second assumption can be applied to the core and
augmentation charges, but not the PS valence ED, since it is ‘delocalized’ and
cannot be assigned to any one atom. For mono-atomic systems, this problem is
trivial because the temperature factor can be applied to the static structure
factor as a whole, but it cannot be performed for systems with more than one
atomic species. Likewise, this problem also manifests in the APW+lo AE DFT
approach of WIEN2k, where its MTs are localized but the interstitial regions
are not. Prior studies using the APW+lo codes have overcome this problem by
applying the average DW factor of the atomic species – the average method – to
the interstitial region [14, 48]; this approximation is sufficient because the
interstitial region makes up a small proportion of the total ED. The valence
electrons in PAW DFT make up a more significant contribution to the total ED
compared to the interstitial region in APW+lo AE DFT (see Table S1 of the
Supplemental Material [38]), such that it may be inaccurate to use an average
DW factor. Hence, in this study, we will instead use the Hirshfeld
partitioning [49] scheme (see Sec. S4 of the Supplemental Material [38]) to
divide the valence electron density into atom-centered densities. This
approach is found to give a small improvement over the average method in terms
of $R^{\textrm{EXP}}$ across all functionals and pseudopotentials tested for
MgO (see Table S2 of the Supplemental Material [38]). It is expected to give
even larger improvements in systems with a large disparity in valence
electrons or DW factors where it can account for the relative contributions of
the different species.

The SFs with thermal vibrations included through the above scheme was compared
to experimental SF for the same three systems as the previous section: Mg
[27], Si [26], and MgO [14]. Experimental details can be found in their
respective references. Thermal effects had already been removed from the
structure factors provided by the experimental study on Mg, so we did not
apply a Debye-Waller factor to this system. Compared to theoretical SFs
computed from DFT, only a small number of SFs (particularly for Mg and MgO)
were available from experiment and only these were compared in the subsequent
analysis. To be consistent with the prior experimental and theoretical
literature on Si, we will compare the ‘effective’ atomic scattering factor of
this system. Its relation to the structure factor is given by:

$f(\mathbf{H})=f(hkl)=\dfrac{F(hkl)}{8\cos(\frac{\pi}{4}(h+k+l))}.$ (29)

This equation was derived for the spherical IAM, which predicts the $(222)$
reflection to be forbidden. However, both DFT and experiment predict
noticeable intensity in this reflection due to the asphericity of the ED, so
it was included as well, with its atomic scattering factor defined as
$f(222)=F(222)/8$ to give the relative contribution from each atom of the
conventional unit cell [26].

###  V.2 Results and discussion

Table 2:  $R^{\textrm{EXP}}$ (%) values for the APW+lo AE DFT method and PAW DFT method employing either the USP and NCP in combination with either AE or PS augmentation charges. |  | APW+lo | AE-USP | PS-USP | AE-NCP | PS-NCP
---|---|---|---|---|---|---
Si | LDA | 0.24 | 0.20 | 0.71 | 0.25 | 1.52
PBE | 0.13 | 0.13 | 0.57 | 0.14 | 1.33
Mg | LDA | 0.48 | 0.45 | 2.17 | 0.53 | 0.86
PBE | 0.36 | 0.34 | 2.03 | 0.42 | 0.74
MgO | LDA | 0.34 | 0.32 | 1.43 | 0.32 | 0.58
PBE | 0.32 | 0.30 | 1.26 | 0.35 | 0.53

Table 2 evaluates $R^{\textrm{EXP}}$ for Si, Mg and MgO from APW+lo WIEN2k and
PAW CASTEP calculations. We have compared PAW results incorporating either AE
or PS augmentation charges. For both USPs and NCPs, the use of PS augmentation
charges leads to an $R^{\textrm{EXP}}$ which can be three to five times larger
than the AE DFT results. This large $R^{\textrm{EXP}}$ is particularly caused
by a systematic underestimation of high Fourier components (see Fig. S3 of the
Supplemental Material [38]). When the AE augmentation charge is used, the
resulting $R^{\textrm{EXP}}$ are all within 0.06 % of the AE DFT results, with
the USPs generally performing better than their NCP counterparts. In fact, the
USPs even appear to be better than AE DFT by up to 0.04 % in
$R^{\textrm{EXP}}$ across the three systems and two DFAs. However, this
improvement is not statistically significant because each experimentally
determined SF has an associated error in its value, which we found to
propagate to an error of around $\pm 0.05\%$ in the $R^{\textrm{EXP}}$ for Si.

In Fig. 2(a), the difference in atomic scattering factors for Si between
experiment and the AE-USP method from CASTEP is plotted for LDA and PBE. We
have focused on Si specifically because there are many available experimental
SFs which have been consolidated to high precision from multiple studies. In
general, both LDA and PBE underestimates the SFs in the range of $s=$ 0.4-0.8
Å-1. The improved $R^{\textrm{EXP}}$ for PBE arises within the 0.4-0.6 Å-1
region, where its underestimation of experimental SFs is less, agreeing with
the observations of Zuo et al [26].

Compared to APW+lo DFT, where a large portion of the valence electrons are
treated in the MT spheres, our approach makes it simple (see Eq. 13) to obtain
the core and valence (PS valence and AE augmentation charge) contributions to
the SF; the difference in these two contributions between PBE and LDA are
plotted for Si in Fig. 2(b). Below 0.4 Å-1, the valence electrons play a
significant contribution to the differences in the total scattering factors,
but this effect becomes small beyond 0.4 Å-1, which arises because the valence
electrons contribute less than 2% to the total scattering factor at that point
(see Fig. S4 of the Supplemental Material [38]). Thus, the marked improvement
in the 0.4-0.6 Å-1 region of PBE over LDA arises predominantly from a better
description of the core electrons. This core ED makes up the frozen-core in
PAW DFT, which, in turn, determines the corresponding pseudopotentials. As it
is difficult to construct pseudopotentials for more sophisticated DFAs such as
hybrid functionals [50], the GGA pseudopotential is often used instead. This
uncontrolled approximation will make it problematic to compare the ED or SFs
of these more advanced DFAs since they share the same (frozen-)core ED as the
DFA used to generate the pseudopotential, so any analysis will not reflect
improvements in the core ED.

##  VI  Conclusion

In conclusion, we have proposed an efficient approach of obtaining accurate
x-ray structure factors for the Vanderbilt ultrasoft pseudopotential and
projector augmented wave methods within DFT. Compared to prior approaches,
involving constructing the total ED on a uniform regular grid, this approach
circumvents such a need by evaluating the core and augmentation charges on
logarithmic radial support grids, significantly reducing the calculation cost
and time; thus, extending the range and size of systems that can be studied.
This approach was implemented in CASTEP and used to study three systems: Si,
Mg and MgO. Comparison of the SFs to all-electron DFT has shown that it is
capable of achieving all-electron accuracy if an all-electron augmentation
charge is used for both norm-conserving and ultrasoft pseudopotentials.
Further comparisons to experimental SFs have shown that our approach can give
a detailed comparison of the ED of different DFAs and where their deficiencies
might occur (i.e. in the core or the valence electrons).

## Data availability

The data that supports the findings of this study are available within the
article and its supplemental material. The input and output files associated
with this study and all analysis can be found on GitHub at
benshi97/Data_XRD_Structure_Factor and Binder.

## Acknowledgments

The authors acknowledge support from the UKCP Consortium, funded by the
Engineering and Physical Sciences Research Council (EPSRC) under grant number
EP/P022561/1. R.J.N. gratefully acknowledges financial support from the EPSRC,
Grant No. EP/L022907/1.

## References

  * Groom  _et al._ [2016] C. R. Groom, I. J. Bruno, M. P. Lightfoot, and S. C. Ward, Acta Crystallogr. B 72, 171 (2016).
  * Berman  _et al._ [2000] H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I. N. Shindyalov, and P. E. Bourne, Nucleic Acids Res. 28, 235 (2000).
  * Hammond [2015] C. Hammond,  _The Basics of Crystallography and Diffraction_ , fourth edition ed., International Union of Crystallography Texts on Crystallography (Oxford University Press, Oxford, New York, 2015).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Nakashima  _et al._ [2011] P. N. H. Nakashima, A. E. Smith, J. Etheridge, and B. C. Muddle, Science 331, 1583 (2011).
  * Zuo  _et al._ [1999] J. M. Zuo, M. Kim, M. O’Keeffe, and J. C. H. Spence, Nature 401, 49 (1999).
  * Forgan  _et al._ [2015] E. M. Forgan, E. Blackburn, A. T. Holmes, A. K. R. Briffa, J. Chang, L. Bouchenoire, S. D. Brown, R. Liang, D. Bonn, W. N. Hardy, N. B. Christensen, M. V. Zimmermann, M. Hücker, and S. M. Hayden, Nat. Commun. 6, 10064 (2015).
  * Mori  _et al._ [2002] H. Mori, S. Lee, A. Yamamoto, S. Tajima, and S. Sato, Phys. Rev. B 65, 092507 (2002).
  * Aubert  _et al._ [2011] E. Aubert, S. Lebègue, M. Marsman, T. T. T. Bui, C. Jelsch, S. Dahaoui, E. Espinosa, and J. G. Ángyán, J. Phys. Chem. A 115, 14484 (2011).
  * Jauch and Reehuis [2005] W. Jauch and M. Reehuis, Acta Crystallogr. A 61, 422 (2005).
  * Jauch and Reehuis [2009] W. Jauch and M. Reehuis, Phys. Rev. B 80, 125126 (2009).
  * Jauch and Reehuis [2011] W. Jauch and M. Reehuis, Phys. Rev. B 83, 115102 (2011).
  * Shibata  _et al._ [1999] S. Shibata, F. Hirota, and T. Shioda, J. Mol. Struct. 485–486, 1 (1999).
  * Zuo  _et al._ [1997a] J. M. Zuo, M. O’Keeffe, P. Rez, and J. C. H. Spence, Phys. Rev. Lett. 78, 4777 (1997a).
  * Ogata  _et al._ [2008] Y. Ogata, K. Tsuda, and M. Tanaka, Acta Crystallogr. A 64, 587 (2008).
  * Peng and Nakashima [2021] D. Peng and P. N. H. Nakashima, Phys. Rev. Lett. 126, 176402 (2021).
  * Wu  _et al._ [2004] L. Wu, Y. Zhu, T. Vogt, H. Su, J. W. Davenport, and J. Tafto, Phys. Rev. B 69, 064501 (2004).
  * Taylor [2003] G. Taylor, Acta Crystallogr. D 59, 1881 (2003).
  * Howard and Jones [1977] C. J. Howard and R. D. G. Jones, Acta Crystallogr. A 33, 776 (1977).
  * Schmidt  _et al._ [1985] M. C. Schmidt, R. Colella, and D. R. Yoder-Short, Acta Crystallogr. A 41, 171 (1985).
  * Schmøkel  _et al._ [2013] M. S. Schmøkel, L. Bjerg, F. K. Larsen, J. Overgaard, S. Cenedese, M. Christensen, G. K. H. Madsen, C. Gatti, E. Nishibori, K. Sugimoto, M. Takata, and B. B. Iversen, Acta Crystallogr. A 69, 570 (2013).
  * Peng and Nakashima [2017] D. Peng and P. N. H. Nakashima, J. Appl. Crystallogr. 50, 602 (2017).
  * Gruene and Mugnaioli [2021] T. Gruene and E. Mugnaioli, Chem. Rev. 121, 11823 (2021).
  * [24] M. Bursch, J.-M. Mewes, A. Hansen, and S. Grimme, Angew. Chem., Int. Ed. Engl. e202205735 (2022).
  * Medvedev  _et al._ [2017] M. G. Medvedev, I. S. Bushmarinov, J. Sun, J. P. Perdew, and K. A. Lyssenko, Science 355, 49 (2017).
  * Zuo  _et al._ [1997b] J. M. Zuo, P. Blaha, and K. Schwarz, J. Phys.: Condens. Matter 9, 7541 (1997b).
  * Friis  _et al._ [2003] J. Friis, G. K. H. Madsen, F. K. Larsen, B. Jiang, K. Marthinsen, and R. Holmestad, J. Chem. Phys. 119, 11359 (2003).
  * Sæterli  _et al._ [2011] R. Sæterli, E. Flage-Larsen, J. Friis, O. M. Løvvik, J. Pacaud, K. Marthinsen, and R. Holmestad, Ultramicroscopy 111, 847 (2011).
  * Aarons  _et al._ [2016] J. Aarons, M. Sarwar, D. Thompsett, and C.-K. Skylaris, J. Chem. Phys. 145, 220901 (2016).
  * Schwerdtfeger [2011] P. Schwerdtfeger, ChemPhysChem 12, 3143 (2011).
  * Vanderbilt [1990] D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
  * Blöchl [1994] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Flage-Larsen  _et al._ [2010] E. Flage-Larsen, O. M. Løvvik, Ø. Prytz, and J. Taftø, Comput. Mater. Sci. 47, 752 (2010).
  * Koelling and Harmon [1977] D. D. Koelling and B. N. Harmon, J. Phys. C: Solid State Phys. 10, 3107 (1977).
  * Laasonen  _et al._ [1993] K. Laasonen, A. Pasquarello, R. Car, C. Lee, and D. Vanderbilt, Phys. Rev. B 47, 10142 (1993).
  * Lejaeghere  _et al._ [2016] K. Lejaeghere, G. Bihlmayer, T. Björkman, P. Blaha, S. Blügel, V. Blum, D. Caliste, I. E. Castelli, S. J. Clark, A. D. Corso, S. de Gironcoli, T. Deutsch, J. K. Dewhurst, I. D. Marco, C. Draxl, M. Dułak, O. Eriksson, J. A. Flores-Livas, K. F. Garrity, L. Genovese, P. Giannozzi, M. Giantomassi, S. Goedecker, X. Gonze, O. Grånäs, E. K. U. Gross, A. Gulans, F. Gygi, D. R. Hamann, P. J. Hasnip, N. a. W. Holzwarth, D. Iuşan, D. B. Jochym, F. Jollet, D. Jones, G. Kresse, K. Koepernik, E. Küçükbenli, Y. O. Kvashnin, I. L. M. Locht, S. Lubeck, M. Marsman, N. Marzari, U. Nitzsche, L. Nordström, T. Ozaki, L. Paulatto, C. J. Pickard, W. Poelmans, M. I. J. Probert, K. Refson, M. Richter, G.-M. Rignanese, S. Saha, M. Scheffler, M. Schlipf, K. Schwarz, S. Sharma, F. Tavazza, P. Thunström, A. Tkatchenko, M. Torrent, D. Vanderbilt, M. J. van Setten, V. V. Speybroeck, J. M. Wills, J. R. Yates, G.-X. Zhang, and S. Cottenier, Science 351, aad3000 (2016).
  * [38] See Supplemental Material at [URL will be inserted by publisher] for computational details, pseudopotential information, a description of the Hirshfeld partitioning scheme for applying Debye-Waller factors and more detailed structure factor comparisons between different materials, DFT functionals and pseudopotential approximations .
  * Shmueli [2001] U. Shmueli, ed.,  _International Tables for Crystallography. Vol. B: Reciprocal Space_ , 2nd ed. (Kluwer Acad. Publ, Dordrecht, 2001).
  * Clark  _et al._ [2005] S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. I. J. Probert, K. Refson, and M. C. Payne, Z. Kristallogr. - Cryst. Mater. 220, 567 (2005).
  * Yates  _et al._ [2007] J. R. Yates, C. J. Pickard, and F. Mauri, Phys. Rev. B 76, 024401 (2007).
  * Blaha  _et al._ [2020] P. Blaha, K. Schwarz, F. Tran, R. Laskowski, G. K. H. Madsen, and L. D. Marks, J. Chem. Phys. 152, 074101 (2020).
  * Coppens [1997] P. Coppens,  _X-Ray Charge Densities and Chemical Bonding_ , International Union of Crystallography Texts on Crystallography No. 4 (International Union of Crystallography, New York, 1997).
  * Bragg [1920] W. H. Bragg, Proc. Phys. Soc. London 33, 304 (1920).
  * Lawrence [1973] J. L. Lawrence, Acta Crystallogr. A 29, 94 (1973).
  * Valvoda and Jeńný [1977] V. Valvoda and J. Jeńný, Phys. Status Solidi A 44, 557 (1977).
  * Malica and Dal Corso [2019] C. Malica and A. Dal Corso, Acta Crystallogr. A 75, 624 (2019).
  * Dudarev  _et al._ [2000] S. L. Dudarev, L.-M. Peng, S. Y. Savrasov, and J.-M. Zuo, Phys. Rev. B 61, 2506 (2000).
  * Hirshfeld [1977] F. L. Hirshfeld, Theoret. Chim. Acta 44, 129 (1977).
  * Yang  _et al._ [2018] J. Yang, L. Z. Tan, and A. M. Rappe, Phys. Rev. B 97, 085130 (2018).


11footnotetext: Univ Rennes, INRIA, IRMAR - UMR 6625, F-35000 Rennes, France.
<EMAIL_ADDRESS> Univ Rennes, INRIA, CNRS, IRMAR -
UMR 6625, F-35000 Rennes, France<EMAIL_ADDRESS>Univ Rennes, CNRS, IRMAR - UMR 6625, F-35000 Rennes, France.
<EMAIL_ADDRESS> Univ Rennes, CNRS, IRMAR - UMR
6625, F-35000 Rennes, France<EMAIL_ADDRESS>The work of I. Almuslimani is supported by the Swiss National Science
Foundation, project No: P2GEP2_195212.

#  Uniformly accurate schemes for drift–oscillatory stochastic differential
equations

Ibrahim Almuslimania,$\ast$, Philippe Chartierb, Mohammed Lemouc and Florian
Méhatsd

###### Abstract.

In this work, we adapt the micro-macro methodology to stochastic differential
equations for the purpose of numerically solving oscillatory evolution
equations. The models we consider are addressed in a wide spectrum of regimes
where oscillations may be slow or fast. We show that through an ad-hoc
transformation (the micro-macro decomposition), it is possible to retain the
usual orders of convergence of Euler-Maruyama method, that is to say, uniform
weak order one and uniform strong order one half. We also show that the same
orders of uniform accuracy can be achieved by a simple integral scheme. The
advantage of the micro-macro scheme is that, in contrast to the integral
scheme, it can be generalized to higher order methods.

Keywords:  highly-oscillatory, stochastic differential equations, micro-macro
decomposition, uniform accuracy.

AMS subject classification (2010):  65L20, 74Q10, 35K15.

##  1\. Introduction

In this paper, we aim at constructing uniformly accurate numerical schemes for
solving Itô stochastic differential equations (SDEs) with a (possibly highly)
oscillatory drift term of the form

$\displaystyle dX(t)=$ $\displaystyle
f_{t/\epsilon}(X(t))dt+\sigma(X(t))dW(t),\quad X(0)=X_{0}\in{\mathbb{R}}^{d},$
(1.1)

where $X(t)$ is a stochastic process with values in ${\mathbb{R}}^{d}$. A
standard assumption in the literature of oscillatory problems and averaging
theory is that the drift function
$(\theta,x)\in{\mathbb{T}}\times{\mathbb{R}}^{d}\mapsto f_{\theta}(x)$ is
assumed to be periodic111 Typically, the deterministic part of (1.1) is a
result of a change of variable applied to the ”non-filtered” equation of the
form $\dot{v}=-\epsilon^{-1}Av+g(v)$, where the matrix $A$ is such that
$\exp(A)=I$, which leads to the periodic vector field $f_{\theta}$ [9, Remark
1.1]. with respect to $\theta$ (we shall denote accordingly the torus by
${\mathbb{T}}={\mathbb{R}}/{\mathbb{Z}}\equiv[0,1]$), see e.g. [19, 9, 8]. For
simplicity and ease of the presentation, and without loss of generality, we
assume in the rest of the paper that
$(\theta,x)\in{\mathbb{T}}\times{\mathbb{R}}^{d}\mapsto f_{\theta}(x)$ is
1-periodic. The diffusion function $\sigma$ is defined as a smooth function
$x\in{\mathbb{R}}^{d}\mapsto\sigma(x)\in{\mathbb{R}}^{d\times m}$. Finally,
$W(t)=(W_{1}(t),...,W_{m}(t))^{T}$ is an array of $m$ independent one-
dimensional Weiner processes. Precise regularity assumptions on $f$ and
$\sigma$ are made in Assumption 1.2. Let us emphasize that $\epsilon$ is here
a parameter whose value can freely vary in the interval $]0,1]$ and that
equation (1.1) is not restricted to its asymptotic regime where $\epsilon$
tends to zero , i.e., the oscillations may be slow or fast without any
restriction on their frequency. This model may appear in many interesting
applications including the perturbation of deterministic oscillatory problems
by adding a random noise term for modeling purposes. For example, in Section
4, we consider a stochastic Hénon-Heiles model, one might also think of
stochastic Fermi-Pasta-Ulam problem and the stochastic nonlinear Schrödinger
equation [21].

Analogously to the case of deterministic differential equations, in the
highly-oscillatory (stiff) regime ($\epsilon\ll 1)$, standard numerical
methods for SDEs, such as Euler-Maruyama method, face a severe time step
restriction ($\Delta t=\mathcal{O}(\epsilon)$) when applied directly to (1.1).
This issue is well documented in the literature for deterministic ODEs [19],
and several classes of methods were introduced in order to deal with stiffness
[5, 8, 10, 11, 12]. However, none of the methods introduced therein qualify as
uniformly accurate methods as they do not produce numerical approximations
with an accuracy and at a cost both independent of the value of
$\epsilon\in]0,1]$. This motivated the introduction of a new methodology based
on averaging techniques and exposed in [9]: there, the authors elaborate a new
technique enabling standard numerical methods to retain their non-stiff order
with uniform accuracy for all $\epsilon\in]0,1]$. In the case of stiff SDEs,
as for ODEs, we can differentiate between two kinds of stiffness: stiff
dissipative SDEs and highly-oscillatory SDEs. While for stiff dissipative
SDEs, many interesting integrators were introduced in the last two decades
with nice stability and convergence properties [1, 2, 4, 3], the numerical
solution of highly-oscillatory (drift) SDEs has not received so much attention
except, for example, [21] in which the author derives weak second order
multirevolution composition methods which are accurate only for very small
values of $\epsilon$ and fail for $\epsilon\simeq 1$. This contribution is, up
to our knowledge, the first attempt to adapt the technique of micro-macro
decomposition to the SDE context in order to construct a _uniformly accurate_
integrator that works well for every $\epsilon\in]0,1]$.

More precisely, we focus , in the present paper, on constructing uniformly
accurate methods for highly-oscillatory SDEs in the spirit of the methodology
explained in [9]. We derive a micro-macro system by introducing a change of
variable that leads us to treat the average decay and the fast oscillations
separately. We show that applying Euler-Maruyama method to the micro-macro
system gives, under appropriate assumptions, an approximation of uniform weak
order 1 and strong order 1/2 for any value of $\epsilon\in(0,1]$ and with no
restriction on the time step. In more mathematical terms, we prove that, under
Assumption 1.2 below, the Euler-Maruyama scheme for solving the micro-macro
system derived from (1.1), provides approximations $X_{n}$ on a uniform grid
$\\{t_{n}\\}_{n=0}^{N}$ such that

$\displaystyle\forall
n=0,\ldots,N,\quad|\mathbb{E}(\phi(X(t_{n})))-\mathbb{E}(\phi(X_{n}))|\leq
C\,h\mbox{ and }\mathbb{E}(|X(t_{n})-X_{n}|^{2})^{\frac{1}{2}}\leq
Ch^{\frac{1}{2}},$

where $h=t_{n+1}-t_{n}$ (for simplicity we consider constant time step size)
and the constant $C$ is independent of $n$, $h$ and $\epsilon$. We prove as
well that the same result can be obtained using the integral scheme introduced
in Section 2. After Noticing the simplicity of implementation and the uniform
accuracy of the integral scheme, a natural question arises: why do we consider
the micro-macro decomposition? The answer is, interestingly, deriving the
micro-macro method gives insight about possible generalizations to higher-
order schemes inspired from deterministic averaging theory developed in [9].

The rest of the paper is organized as follows. In the next subsection, we fix
the main notations used throughout the paper and we make assumptions. In
Section 2 we introduce the integral scheme and we prove its uniform
convergence (weak and strong) with respect to the parameter $\epsilon$. In
Section 3 we derive the micro-macro scheme, and we show its uniform
convergence properties. Finally, in Section 4 we present some numerical
experiments that illustrate the efficiency of the schemes.

### Notations and Assumptions

Noteworthy, our results are obtained under quite standard assumptions that we
recall below. In particular, we shall constantly suppose that Assumption 1.2
below holds true in the sequel.

######  Definition 1.1.

Let $q\in\mathbb{N}$, we define $C_{poly}^{q}$ as the set of $q$-times
differentiable functions from ${\mathbb{R}}^{d}\to{\mathbb{R}}$ whose all
derivatives up to order $q$ have at most polynomial growth. We define as well
$C_{lin}^{q}$ in the same way but with linear growth instead of polynomial.

######  Assumption 1.2.

The functions $x\mapsto\sigma_{i,j}(x)$, $i=1,\ldots,d$, $j=1,\ldots,m$ are
Lipschitz functions of the set $C^{2}_{poly}$. The function $(\theta,x)\mapsto
f_{\theta}(x)\in{\mathbb{R}}^{d}$ is defined on
${\mathbb{T}}\times{\mathbb{R}}^{d}$ and it is uniformly Lipschitz continuous
with respect to $x$ (i.e with Lipschitz constant independent of $\theta$). In
addition, its components $(f_{\theta}(\cdot))_{i}$, $i=1,\ldots,d$, are
functions of $C_{poly}^{4}\cap C_{lin}^{2}$. Furthermore, $f$ is continuously
differentiable with respect to $\theta$ and
$x\mapsto\partial_{\theta}f_{\theta}(x)$ is uniformly Lipschitz continuous
with respect to the variable $x$.

######  Remark 1.3.

The Lipschitz continuity of $f_{\theta}$ and $\sigma$ ensures the existence
and uniqueness of the solution of (1.1) in $C([0,T];L^{2}(\Omega))$. It
ensures in addition a linear growth, that is to say

$\forall(\theta,x)\in{\mathbb{T}}\times{\mathbb{R}}^{d},\quad|f_{\theta}(x)|+|\sigma(x)|\leq
K(1+|x|),$

where $K$ is a constant independent of $\theta$.

######  Notations.

The derivative with respect to the space variable $x$ will be denoted by a
_prime_. For example, $f^{\prime}_{\theta}(x)=\partial_{x}f_{\theta}(x)$ (and
gradient or Jacobien when needed), and
$f^{\prime\prime}_{\theta}(x)=\partial^{2}_{x}f_{\theta}(x)$ (Hessian matrix
or second derivative tensor when needed). The letter $L$ will be used
throughout the paper as the maximum Lipschitz constant of the concerned
functions. The constants represented by the capital letter $C$ (indexed or
not) are generic constants. We denote by $|\cdot|$ the Euclidean norm of
${\mathbb{R}}^{d}$.

We will conduct all our analysis for $m=1$ in (1.1). This comes for the sake
of simplicity and ease of presentation, and to avoid any confusion due to the
heavy computations done in the proofs. However, we emphasize that our results
apply straightforwardly to SDEs of the form (1.1) with any number of Weiner
processes. Hence, $\sigma$ will denote a function from
${\mathbb{R}}^{d}\to{\mathbb{R}}^{d}$ and $W(t)$ will denote a standard one-
dimensional Weiner process.

##  2\. Integral scheme

Let $T>0$ be a given final time, $N\in\mathbb{N}\backslash\\{0\\}$, and
$h=T/N$, and let $t_{n}=nh,~{}n=0,\dots,N$. We consider the integral scheme

$X_{n+1}=X_{n}+h\left<f\right>(X_{n})+\epsilon\int_{t_{n}/\epsilon}^{t_{n+1}/\epsilon}(f_{\theta}(X_{n})-\left<f\right>(X_{n}))d\theta+\sigma(X_{n})\Delta
W_{n},$ (2.1)

where $\left<f\right>(x)$ is the value at $x$ of the average of $f$ defined as

$\left<f\right>(x)=\int_{\mathbb{T}}f_{\theta}(x)d\theta,$

and $\Delta W_{n}=W(t_{n+1})-W(t_{n})\sim\mathcal{N}(0,h)$ which can be
replaced by $\sqrt{h}\xi_{n}$ where $\xi_{n}\sim\mathcal{N}(0,1)$. Note that
the integral term in (2.1) can be rewritten as
$\epsilon(F_{t_{n+1}/\epsilon}(X_{n})-F_{t_{n}/\epsilon}(X_{n}))$, where

$F_{\theta}(x)=\int_{0}^{\theta}\left(f_{\tau}(x)-\left<f\right>(x)\right)d\tau.$

######  Theorem 2.1.

Consider the integral scheme (2.1), where $t_{n}=nh$ for $n=0,\ldots,N$, and
$N\in\mathbb{N}\backslash\\{0\\}$. Let $X(t)$ be the solution of (1.1). Then,
for all $\phi\in C_{poly}^{4}$, there exists $C>0$ independent of $h$, $n$ and
$\epsilon\in(0,1]$ such that,

$\displaystyle\forall n\in\\{0,\ldots,N\\},\quad$
$\displaystyle|\mathbb{E}(\phi(X(t_{n})))-\mathbb{E}(\phi(X_{n}))|\leq Ch,$
(2.2) $\displaystyle\forall n\in\\{0,\ldots,N\\},\quad$
$\displaystyle\mathbb{E}(|X(t_{n})-X_{n}|^{2})^{\frac{1}{2}}\leq
Ch^{\frac{1}{2}}.$ (2.3)

######  Lemma 2.2.

Under the assumptions of Theorem 2.1, one has for all $n=0,\dots,N-1$

$\left|\mathbb{E}\left(\int_{t_{n}}^{t_{n+1}}[f_{t/\epsilon}(X(t))-f_{t/\epsilon}(X(t_{n}))]dt\right)\right|\leq
Ch^{2}.$ (2.4)

###### Proof.

We first rewrite the integrand $f_{t/\epsilon}(X(t))-f_{t/\epsilon}(X(t_{n}))$
of (2.4) as

$\displaystyle\int_{t_{n}}^{t}d(f_{t/\epsilon}(X(s)))$
$\displaystyle=\int_{t_{n}}^{t}[f^{\prime}_{t/\epsilon}(X(s))f_{t/\epsilon}(X(s))$
$\displaystyle+\frac{1}{2}f^{\prime\prime}_{t/\epsilon}(X(s))(\sigma(X(s)),\sigma(X(s)))]ds$
$\displaystyle+\int_{t_{n}}^{t}f^{\prime}_{t/\epsilon}(X(s))\sigma(X(s))dW(s),$

so that

$\displaystyle\left|\mathbb{E}\left(\int_{t_{n}}^{t_{n+1}}[f_{t/\epsilon}(X(t))-f_{t/\epsilon}(X(t_{n}))]dt\right)\right|$
$\displaystyle\leq\left|\mathbb{E}\left(\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{t}\biggl{[}f^{\prime}_{t/\epsilon}(X(s))f_{t/\epsilon}(X(s))\right.\right.$
$\displaystyle\left.\left.+\frac{1}{2}f^{\prime\prime}_{t/\epsilon}(X(s))(\sigma(X(s)),\sigma(X(s)))\biggr{]}dsdt\right)\right|$
$\displaystyle+\left|\mathbb{E}\left(\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{t}f^{\prime}_{t/\epsilon}(X(s))\sigma(X(s))dW(s)dt\right)\right|,$

and using Fubini’s theorem, we get the following upper-bound of the left-hand
side of (2.4)

$\displaystyle\left|\int_{t_{n}}^{t_{n+1}}\mathbb{E}\left(\int_{t_{n}}^{t}f^{\prime}_{t/\epsilon}(X(s))\sigma(X(s))dW(s)\right)dt\right|+\mathcal{O}(h^{2}),$

where $\mathcal{O}(h^{2})$ is bounded thanks to the polynomial growth
assumptions and the boundedness of the moments. The expectation in previous
integral is an expectation of a stochastic integral, hence it is null since
$X(s)$ is independent of the increments of $W$ for times above $s$. ∎

######  Lemma 2.3.

Under the assumptions of Theorem 2.1, we have for all $n=0,\dots,N-1$

$\mathbb{E}(|X(t)-X(t_{n})|^{2}~{};~{}X(t_{n})=x)\leq C(x)(t-t_{n}).$

###### Proof.

From the integral form of equation (1.1)

$X(t)=X(t_{n})+\int_{t_{n}}^{t}f_{s/\epsilon}(X(s)ds+\int_{t_{n}}^{t}\sigma(X(s))dW(s).$

we may write

$\displaystyle\mathbb{E}(|X(t)-X(t_{n})|^{2}~{};~{}X(t_{n})=x)$
$\displaystyle=\mathbb{E}\left(\left|\int_{t_{n}}^{t}f_{s/\epsilon}(X(s))ds+\int_{t_{n}}^{t}\sigma(X(s))dW(s)\right|^{2}\right)$
$\displaystyle\leq
2(t-t_{n})\int_{t_{n}}^{t}\mathbb{E}(f_{s/\epsilon}(X(s))^{2})ds+2\int_{t_{n}}^{t}\mathbb{E}(\sigma(X(s))^{2})ds$
$\displaystyle\leq
2(t-t_{n})\int_{t_{n}}^{t}K(1+\mathbb{E}(|X(t_{n})|^{2}))ds+2\int_{t_{n}}^{t}K(1+\mathbb{E}(|X(t_{n})|^{2}))ds$
$\displaystyle\leq C(x)(t-t_{n}),$

where we have used, between the first and second lines, the Young’s
inequality, then the Cauchy-Schwarz inequality for the first integral, and Itô
isometry for the second integral. Finally, we have used the bound
$\mathbb{E}(|X(s)|^{2})\leq K(1+|x|^{2})$, borrowed from [13], between the
third and the fourth lines. ∎

###### Proof of Theorem 2.1.

We hereby follow the methodology introduced in [20], which consists in
bounding the moments (first step), then proving the weak convergence of the
scheme (second step) and finally establishing its strong convergence (third
step).

### Step $1$

In order to bound the moments of $X_{n}$ of arbitrary order, we shall resort
to [18, Lemma 2.2, p. 102], which requires the following estimates

$\displaystyle|\mathbb{E}(X_{n+1}-X_{n}|X_{n}=x)|$
$\displaystyle=\left|\mathbb{E}\left(\epsilon\int_{t_{n}/\epsilon}^{t_{n+1}/\epsilon}f_{\theta}(x)d\theta\right)\right|$
$\displaystyle\leq\epsilon\int_{t_{n}/\epsilon}^{t_{n+1}/\epsilon}K(1+|x|)d\theta=K(1+|x|)h,$

and

$\displaystyle|X_{n+1}-X_{n}|$
$\displaystyle=\left|\epsilon\int_{t_{n}/\epsilon}^{t_{n+1}/\epsilon}f_{\theta}(X_{n})d\theta+\sigma(X_{n})\Delta
W_{n}\right|$ $\displaystyle\leq K(1+|X_{n}|)h+K(1+|X_{n}|)\frac{|\Delta
W_{n}|}{\sqrt{h}}\sqrt{h}\leq M_{n}(1+|X_{n}|)\sqrt{h},$

where $M_{n}=K\left(\sqrt{h}+\frac{|\Delta W_{n}|}{\sqrt{h}}\right)$ which is
of bounded moments since $\frac{\Delta W_{n}}{\sqrt{h}}\sim\mathcal{N}(0,1)$.

### Step $2$

Let $\phi\in C^{4}_{poly}$ be a test function. For $X(t_{n})=X_{n}=x$, by
performing Taylor expansion of $\phi$ with integral remainder, it can be shown
that

$\displaystyle|\mathbb{E}(\phi(X(t_{n+1}))-\phi(X_{n+1}))|\leq\biggl{|}\mathbb{E}\left[\phi^{\prime}(x)\left(\int_{t_{n}}^{t_{n+1}}(f_{t/\epsilon}(X(t))-f_{t/\epsilon}(X(t_{n})))dt\right)\right]\biggr{|}$
$\displaystyle+\biggl{|}\mathbb{E}\left[\phi^{\prime\prime}(x)\left(\int_{t_{n}}^{t_{n+1}}(\sigma(X(t))-\sigma(X(t_{n})))dW(t)\right)\left(\int_{t_{n}}^{t_{n+1}}(\sigma(X(t))+\sigma(X(t_{n})))dW(t)\right)\right]\biggr{|}$
$\displaystyle+\frac{1}{6}\biggl{|}\mathbb{E}\left[\int_{0}^{1}(1-\tau)^{3}\phi^{(4)}(x+\tau(X(t_{n+1})-x))d\tau\,(X(t_{n+1})-x)^{4}\right]\biggr{|}$
$\displaystyle+\frac{1}{6}\biggl{|}\mathbb{E}\left[\int_{0}^{1}(1-\tau)^{3}\phi^{(4)}(x+\tau(X_{n+1}-x))d\tau\,(X_{n+1}-x)^{4}\right]\biggr{|}+Ch^{2}$
$\displaystyle\coloneqq\text{I+II+III+IV}+Ch^{2},$

where the term $Ch^{2}$ comes from the remaining expectations of the second
and third derivatives of $\phi$ applied to the integrals, which are zero for
odd number of stochastic integrals, and bounded by $Ch^{2}$ otherwise. The
constant C is independent of $\epsilon$ thanks to the boundedness of the
moments of the exact and numerical solutions, and to the regularity
assumptions made in Assumption 1.2. By Lemma 2.2, Lemma 2.3 and the
boundedness of the moments, the terms I, III and IV are $\mathcal{O}(h^{2})$.
As for the integrand of the second term, we have using Itô formula

$\displaystyle\sigma(X(t))-\sigma(X(t_{n}))$
$\displaystyle=\int_{t_{n}}^{t}\biggl{(}\sigma^{\prime}(X(s))f_{s/\epsilon}(X(s))\bigg{.}$
$\displaystyle\bigg{.}+\frac{1}{2}\sigma^{\prime\prime}(X(s))(\sigma(X(s)),\sigma(X(s)))\biggr{)}ds$
$\displaystyle+\int_{t_{n}}^{t}\sigma^{\prime}(X(s))\sigma(X(s))dW(s).$

Hence, using the Lipschitz-continuity of $\sigma$, the polynomial growth of
$\phi^{\prime\prime}$, Lemma 2.3, and the following consequence of Itô
isometry

$\forall g,h\in
L^{2}_{ad}([a,b]\times\Omega),~{}\mathbb{E}\left[\left(\int_{a}^{b}g(t)dW(t)\right)\left(\int_{a}^{b}h(t)dW(t)\right)\right]=\mathbb{E}\left[\int_{a}^{b}g(t)h(t)dt\right]~{}~{}$

we have

II $\displaystyle\leq
C(x)\biggl{|}\mathbb{E}\left[\left(\int_{t_{n}}^{t_{n+1}}(\sigma(X(t))-\sigma(X(t_{n})))dW(t)\right)\left(\int_{t_{n}}^{t_{n+1}}(\sigma(X(t))-\sigma(X(t_{n})))dW(t)\right.\right.$
$\displaystyle\left.\left.+2\int_{t_{n}}^{t_{n+1}}\sigma(X(t_{n}))dW(t)\right)\right]\biggr{|}$
$\displaystyle=C(x)\biggl{|}\mathbb{E}\left[\left(\int_{t_{n}}^{t_{n+1}}(\sigma(X(t))-\sigma(X(t_{n})))^{2}dt\right)\right]\biggr{|}$
$\displaystyle+2C(x)\biggl{|}\mathbb{E}\left[\left(\int_{t_{n}}^{t_{n+1}}\sigma(X(t_{n}))(\sigma(X(t))-\sigma(X(t_{n})))dt\right)\right]\biggr{|}$
$\displaystyle\leq
C(x)L\int_{t_{n}}^{t_{n+1}}C(t-t_{n})dt+2\left|\mathbb{E}\left[\sigma(X(t_{n}))\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{t}d\sigma(X(s))dt\right]\right|$
$\displaystyle=2\left|\mathbb{E}\left[\sigma(X(t_{n}))\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{t}\biggl{(}\sigma^{\prime}(X_{s})f_{s/\epsilon}(X_{s})+\sigma^{\prime\prime}(X_{s})(\sigma(X_{s}),\sigma(X_{s}))\biggr{)}dsdt\right]\right|$
$\displaystyle+2\left|\mathbb{E}\left[\sigma(X(t_{n}))\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{t}\sigma^{\prime}(X_{s})\sigma(X_{s})dW(s)dt\right]\right|+C_{1}h^{2}.$

The first expectation after the last equal sign is clearly
$\mathcal{O}(h^{2})$. The second one is equal to zero since $X(t_{n})$ is
independent of the increments $W(r)-W(s)$ for $t_{n}\leq s\leq r$. The
boundedness of the moments and the local weak order $2$ imply the global weak
convergence of order $1$ by a theorem from [17] (see also [18, Chapter 2.2]).

### Step 3

We first derive an upper-bound of
$\mathbb{E}(|X(t_{n+1})-X_{n+1}|^{2};X(t_{n})=x)$ as follows

$\displaystyle\mathbb{E}\biggl{(}\biggl{|}\int_{t_{n}}^{t_{n+1}}(f_{t/\epsilon}(X(t))-f_{t/\epsilon}(x))dt+\int_{t_{n}}^{t_{n+1}}(\sigma(X(t))-\sigma(x))dW(t)\biggr{|}^{2}\biggr{)}$
$\displaystyle\leq
2\mathbb{E}\biggl{(}\biggl{|}\int_{t_{n}}^{t_{n+1}}(f_{t/\epsilon}(X(t))-f_{t/\epsilon}(x))dt\biggr{|}^{2}\biggr{)}+2\mathbb{E}\biggl{(}\biggl{|}\int_{t_{n}}^{t_{n+1}}(\sigma(X(t))-\sigma(x))dW(t)\biggr{|}^{2}\biggr{)}$
$\displaystyle\leq
2h\mathbb{E}\biggl{(}\int_{t_{n}}^{t_{n+1}}\biggl{|}f_{t/\epsilon}(X(t))-f_{t/\epsilon}(X(t_{n}))\biggr{|}^{2}dt\biggr{)}+2\mathbb{E}\biggl{(}\int_{t_{n}}^{t_{n+1}}\biggl{|}\sigma(X(t))-\sigma(X(t_{n}))\biggr{|}^{2}dt\biggr{)}$
$\displaystyle\leq
2L^{2}(h+1)\int_{t_{n}}^{t_{n+1}}\mathbb{E}(|X(t)-X(t_{n})|^{2})dt\leq
2L^{2}(h+1)\int_{t_{n}}^{t_{n+1}}C(t-t_{n})dt\leq Ch^{2},$

where we have used Lemma (2.3). Finally, we conclude that
$\mathbb{E}(|X(t_{n+1})-X_{n+1}|^{2})^{\frac{1}{2}}\leq Ch$. A well-known
theorem by Milstein [16] then allows to establish inequality (2.3). ∎

Remark that for the integral scheme we can relax the assumption
$(f_{\theta}(.))_{i},\,i=1,\dots,d\in C_{poly}^{4}\cap C_{lin}^{2}$ to
$(f_{\theta}(.))_{i},\,i=1,\dots,d\in C_{poly}^{2}$.

##  3\. Micro-Macro method

The integral scheme, despite its uniform accuracy, does not generalize to
higher order methods. For this reason, and inspired by [8, 9], we will
introduce and analyze a micro-macro separation of scales that leads the
standard Euler-Maruyama method to retain its usual weak and strong orders of
convergence. We believe that this will be an excellent starting point and a
good insight to try to develop higher order micro-macro methods in future
works.

As in general averaging theory, the purpose is to to find a periodic, near-
identity and smooth change of variable $\Phi_{\theta}$, together with a flow
$\Psi_{t}$, the flow map of an autonomous non-stiff differential equation on
$\mathbb{R}^{d}$, such that the solution of the original equation (1.1) takes
the composed form [19, 15, 8, 9]

$X(t)=\Phi_{t/\epsilon}\circ\Psi_{t}(X_{0}).$ (3.1)

In the current work, we restrict ourselves to first order averaging, and we
show that the first order change of variable $\Phi_{\theta}$ defined, for all
$0<\epsilon\leq 1$, by the formula

$\displaystyle\Phi_{\theta}(x)=x+\epsilon
F_{\theta}(x)=x+\epsilon\int_{0}^{\theta}\left(f_{\tau}(x)-\left<f\right>(x)\right)d\tau,$
(3.2)

and the first order flow $\Psi_{t}$ satisfying the autonomous equation

$d\Psi_{t}(x)=\left<f\right>(\Psi_{t}(x))dt+\sigma(\Psi_{t}(x))dW(t),\qquad\Psi_{0}(x)=x,$
(3.3)

both borrowed from deterministic averaging, work well to construct a _uniform
accurate micro-macro scheme of weak order 1 and strong order 1/2_ for (1.1).
For the derivation of (3.2) and (3.3) from (3.1), we refer to [9].

Let us denote by $\underline{X}$ the solution of (3.3), separating slow and
fast scales as follows

$\displaystyle X=\Phi_{t/\epsilon}(\underline{X})+Y,$ (3.4)

leads, using Itô formula, to the micro-macro system of the form

$\displaystyle d\underline{X}$
$\displaystyle=\left<f\right>(\underline{X})dt+\sigma(\underline{X})dW,\quad\underline{X}(0)=X_{0},$
(3.5) $\displaystyle dY$
$\displaystyle=\left(f_{t/\epsilon}\left(\Phi_{t/\epsilon}(\underline{X})+Y\right)-f_{t/\epsilon}(\underline{X})-\epsilon
F_{t/\epsilon}^{\prime}(\underline{X})\left<f\right>(\underline{X})\right.$
(3.6)
$\displaystyle\left.-\frac{\epsilon}{2}F_{t/\epsilon}^{\prime\prime}(\underline{X})\left(\sigma(\underline{X}),\sigma(\underline{X})\right)\right)dt$
$\displaystyle+\left(\sigma(\Phi_{t/\epsilon}(\underline{X})+Y)-\sigma(\underline{X})-\epsilon
F_{t/\epsilon}^{\prime}(\underline{X})\sigma(\underline{X})\right)dW,\quad
Y(0)=0.$

Here $Y$ represents the averaging error that is characterized by the equation
(3.6) and added after applying the change of variable in order to recover the
solution $X(t)$ of (1.1). We will show later that the expected value of $Y(t)$
is of size $\epsilon$ (see Lemma 3.3), which is analogous to the averaging
error for deterministic oscillatory equations ($\sigma\equiv 0$).

######  Remark 3.1.

The existence and the uniqueness of the solution $Y(t)$ of (3.6) is an
immediate consequence of the existence and the uniqueness of the solutions
$X(t)$ and $\underline{X}(t)$ of (1.1) and (3.5), and the decomposition (3.4).

We use the same uniform discretization as for the integral scheme. Our aim is
now to prove a uniform (in $\epsilon$) convergence result for the following
micro-macro scheme, which is nothing but the Euler-Maruyama method applied to
(3.5, 3.6)

$\displaystyle\underline{X}_{n+1}$
$\displaystyle=\underline{X}_{n}+h\left<f\right>(\underline{X}_{n})+\sqrt{h}\sigma(\underline{X}_{n})\xi_{n},\quad\underline{X}_{0}=X_{0},$
(3.7) $\displaystyle Y_{n+1}$
$\displaystyle=Y_{n}+h\left(f_{t/\epsilon}\left(\Phi_{t_{n}/\epsilon}(\underline{X}_{n})+Y\right)-f_{t_{n}/\epsilon}(\underline{X}_{n})\right)$
(3.8) $\displaystyle-\epsilon
h\left(F_{t_{n}/\epsilon}^{\prime}(\underline{X}_{n})\left<f\right>(\underline{X}_{n})+\frac{1}{2}F_{t_{n}/\epsilon}^{\prime\prime}(\underline{X}_{n})\left(\sigma(\underline{X}_{n}),\sigma(\underline{X}_{n})\right)\right)$
$\displaystyle+\sqrt{h}\left(\sigma(\Phi_{t_{n}/\epsilon}(\underline{X}_{n})+Y_{n})-\sigma(\underline{X}_{n})-\epsilon
F_{t_{n}/\epsilon}^{\prime}(\underline{X}_{n})\sigma(\underline{X}_{n})\right)\xi_{n},\quad
Y_{0}=0,$

where the increment $\xi_{n}$ is a random quantity sampled from a normalized
Gaussian centered at zero and with variance $1$.

###  3.1. Main result

######  Theorem 3.2.

Consider the Euler-Maruyama scheme (3.7, 3.8) for solving the micro-macro
system (3.5, 3.6) and let

$X_{n}=\Phi_{t_{n}/\epsilon}(\underline{X}_{n})+Y_{n}$

for $n=0,\dots,N$. Let $X(t)$ be the solution of (1.1). Then, for all $\phi\in
C_{poly}^{4}$, there exists $C>0$ independent of $h$, $n$, and $\epsilon$,
such that, one has

$\displaystyle\forall n\in\\{0,\ldots,N\\},\quad$
$\displaystyle|\mathbb{E}(\phi(X(t_{n})))-\mathbb{E}(\phi(X_{n}))|\leq C\,h,$
(3.9) $\displaystyle\forall n\in\\{0,\ldots,N\\},\quad$
$\displaystyle\mathbb{E}(|X(t_{n})-X_{n}|^{2})^{\frac{1}{2}}\leq
Ch^{\frac{1}{2}}.$ (3.10)

Again, the proof of the theorem follows the usual steps from [20]. Note that
the main novelty of our result lies in the fact that estimate (3.9) is uniform
w.r.t. $\epsilon\in(0,1]$. The following lemma will be needed in the proof of
Theorem 3.2.

######  Lemma 3.3.

There exists $c>0$ such that for all $\epsilon\in(0,1]$ and all $t\in[0,T]$,
one has

$\mathbb{E}(|Y(t)|)\leq c\epsilon.$

###### Proof.

First, note that we have for all $t\in[0,T]$ and all positive integer $m$ (see
[13])

$\mathbb{E}(|\underline{X}(t)|^{2m})\leq K(1+|X_{0}|^{2m}).$

Moreover, we have

$\displaystyle\mathbb{E}(|Y(t)|^{2})\leq
2t\mathbb{E}\left(\int_{0}^{t}\biggl{|}f_{s/\epsilon}(\Phi_{s/\epsilon}(\underline{X}(s))+Y(s))-f_{s/\epsilon}(\underline{X}(s))\right.$
$\displaystyle\left.-\epsilon\biggl{(}F^{\prime}_{s/\epsilon}(\underline{X}(s))\left<f\right>(\underline{X}(s))+\frac{1}{2}F^{\prime\prime}_{s/\epsilon}(\underline{X}(s))(\sigma(\underline{X}(s)),\sigma(\underline{X}(s)))\biggr{)}\biggr{|}^{2}ds\right)$
$\displaystyle+2\mathbb{E}\left(\int_{0}^{t}\biggl{|}\sigma(\Phi_{s/\epsilon}(\underline{X}(s))+Y(s))-\sigma(\underline{X}(s))-\epsilon
F^{\prime}_{s/\epsilon}(\underline{X}(s))\sigma(\underline{X}(s))\biggr{|}^{2}ds\right)$
$\displaystyle\leq
4t\int_{0}^{t}\mathbb{E}\left(\biggl{|}f_{s/\epsilon}(\Phi_{s/\epsilon}(\underline{X}(s))+Y(s))-f_{s/\epsilon}(\underline{X}(s))\biggr{|}^{2}\right)ds$
$\displaystyle+4t\epsilon^{2}\int_{0}^{t}\mathbb{E}\left(\biggl{|}F^{\prime}_{s/\epsilon}(\underline{X}(s))\left<f\right>(\underline{X}(s))+\frac{1}{2}F^{\prime\prime}_{s/\epsilon}(\underline{X}(s))(\sigma(\underline{X}(s)),\sigma(\underline{X}(s)))\biggr{|}^{2}\right)ds$
$\displaystyle+4\int_{0}^{t}\mathbb{E}\left(\biggl{|}\sigma(\Phi_{s/\epsilon}(\underline{X}(s))+Y(s))-\sigma(\underline{X}(s))\biggr{|}^{2}\right)ds+4\epsilon^{2}\int_{0}^{t}\mathbb{E}\left(\biggl{|}F^{\prime}_{s/\epsilon}(\underline{X}(s))\sigma(\underline{X}(s))\biggr{|}^{2}\right)ds$
$\displaystyle\coloneqq 4tA+4t\epsilon^{2}B+4C+4\epsilon^{2}D.$

In the first inequality, we have used the Cauchy-Schwarz inequality for the
first term and Itô isometry for the second one. We have also used the
triangular inequality as well as Young’s inequality. Given that
$\mathbb{E}(|\underline{X}(s)|^{2m})<+\infty$ and that the components of
$F_{\theta}$ are in $C^{2}_{poly}$ by Assumption 1.2, we can conclude that the
terms $B$ and $D$ are uniformly bounded. Now, by definition of the change of
variables (3.2) and the uniform Lipschitz continuity of $f_{\theta}$ and
$\sigma$, we get

$\displaystyle 4tA+4C$ $\displaystyle\leq
4tL^{2}\int_{0}^{t}\mathbb{E}(|\epsilon
F^{\prime}_{s/\epsilon}(\underline{X}(s))+Y(s)|^{2})ds+4L^{2}\int_{0}^{t}\mathbb{E}(|\epsilon
F^{\prime}_{s/\epsilon}(\underline{X}(s))+Y(s)|^{2})ds$ $\displaystyle\leq
8L^{2}(T+1)\epsilon^{2}\int_{0}^{t}\mathbb{E}(|F^{\prime}_{s/\epsilon}(\underline{X}(s))|^{2})ds+8L^{2}(T+1)\int_{0}^{t}\mathbb{E}(|Y(s)|^{2})ds$
$\displaystyle\leq
C_{1}\epsilon^{2}+C_{2}\int_{0}^{t}\mathbb{E}(|Y(s)|^{2})ds.$

Therefore, $\mathbb{E}(|Y(t)|^{2})$ satisfies the following inequality

$\mathbb{E}(|Y(t)|^{2})\leq
C_{3}\epsilon^{2}+C_{4}\int_{0}^{t}\mathbb{E}(|Y(s)|^{2})ds.$

Gronwall’s lemma imply that $\mathbb{E}(|Y(t)|^{2})\leq c^{2}\epsilon^{2}$.
Finally, $\mathbb{E}(|Y(t)|)^{2}\leq\mathbb{E}(|Y(t)|^{2})=c^{2}\epsilon^{2}$,
thus $\mathbb{E}(|Y(t)|)\leq c\epsilon$ where $c$ is a positive constant
independent of $\epsilon$ and $t$.

∎

###### Proof of Theorem 3.2.

We begin by the proof of the weak convergence: Collecting the variables
$\underline{X}$ and $Y$ into $Z=(\underline{X},Y)$ and similarly
$z=(\underline{x},y)$, we may rewrite the micro-macro system (3.5, 3.6) as

$dZ=g_{t/\epsilon}(Z)dt+\Sigma_{t/\epsilon}(Z)dW$

where

$g_{\theta}(z)=\left(\begin{array}[]{c}\left<f\right>(\underline{x})\\\
f_{\theta}\left(\Phi_{\theta}(\underline{x})+y\right)-f_{\theta}(\underline{x})-\epsilon
F_{\theta}^{\prime}(\underline{x})\left<f\right>(\underline{x})-\frac{\epsilon}{2}F_{\theta}^{\prime\prime}(\underline{x})\left(\sigma(\underline{x}),\sigma(\underline{x})\right)\end{array}\right)$

and

$\Sigma_{\theta}(z)=\left(\begin{array}[]{c}\sigma(\underline{x})\\\
\sigma(\Phi_{\theta}(\underline{x})+y)-\sigma(\underline{x})-\epsilon
F_{\theta}^{\prime}(\underline{x})\sigma(\underline{x})\end{array}\right)$

and accordingly equations (3.7), (3.8) as

$Z_{n+1}=Z_{n}+hg_{t_{n}/\epsilon}(Z_{n})+\sqrt{h}\Sigma_{t_{n}/\epsilon}(Z_{n})\xi_{n}$

where the $\xi_{n}\sim\mathcal{N}(0,I_{m})$ are independent Gaussian random
variables.

### Step $1$

Since the vector fields $g_{t/\epsilon}$ and $\Sigma_{t/\epsilon}$ have
uniform linear growth in $z$ (this follows from the linear growth of $f$,
$\sigma$, $\Phi$, $F$, $F^{\prime}$, and $F^{\prime\prime}$), the proof of the
boundedness of the moments of the numerical solution can be obtained following
the same arguments as in Milstein’s lemma for the non-oscillating case
$(\epsilon=1)$. As a matter of fact, we have that

  * (i)

$\mathbb{E}(|Z_{0}|^{2m})<+\infty$ since $Z_{0}=(X_{0},0)\in{\mathbb{R}}^{2d}$
(deterministic);

  * (ii)

$|\mathbb{E}(Z_{n+1}-Z_{n}~{}|~{}Z_{n}=z)|=h|g_{t_{n}/\epsilon}(z)|\leq
K(1+|z|)h$;

  * (iii)

$|Z_{n+1}-Z_{n}|\leq
h|g_{t_{n}/\epsilon}(Z_{n})|+|\Sigma_{t_{n}/\epsilon}(Z_{n})||\Delta
W_{n}|\leq M_{n}(1+|Z_{n}|)\sqrt{h}$, with $\Delta
W_{n}=W(t_{n+1})-W(t_{n})\sim\mathcal{N}(0,h)$ and $M_{n}=K|\frac{\Delta
W_{n}}{\sqrt{h}}|$ has clearly bounded moments uniformly in $n$.

We have used the fact that $g$ and $\Sigma$ grow linearly in $Z$ thanks to
Assumption 1.2. Under the above conditions, Lemma 2.2 from [18, p. 102]
implies the boundedness of the moments of arbitrary order of $Z_{n}$.

### Step $2$

We define, for a function $\psi$ that depends explicitly on $\theta$ and $z$,

$({\mathcal{L}}_{\theta}\psi)(z)=\partial_{\theta}\psi(\theta,z)+\partial_{z}\psi(\theta,z)g_{\theta}(z)+\frac{1}{2}\partial^{2}_{z}{\psi}(\theta,z)\left(\Sigma_{\theta}(z),\Sigma_{\theta}(z)\right).$
(3.11)

We have for any test function $\phi\in C_{poly}^{4}$ (independent of theta)

$\displaystyle\mathbb{E}(\phi(Z_{n+1})|Z_{n}=z)$
$\displaystyle\leq\phi(z)+h(\mathcal{L}_{t_{n}/\epsilon}\phi)(z)+Ch^{2}$
$\displaystyle+\frac{1}{6}\mathbb{E}\left(\left.\int_{0}^{1}(1-\tau)^{3}\phi^{(4)}(z+\tau(Z_{n+1}-z))d\tau\;(Z_{n+1}-z)^{4}\right|Z_{n}=z\right).$

where, as in the proof of Theorem 2.1, the constant $C$ is independent of
$\epsilon$ and where the term $Ch^{2}$ comes from the remaining expectations
of the second and third derivatives of $\phi$ applied repeatedly to
$hg+\sqrt{h}\Sigma\xi_{n}$, which are zero for odd moments of $\xi_{n}$, and
bounded by $Ch^{2}$ otherwise. For the part of the remainder coming from the
above first order Taylor expansion of the expectation of $\phi(Z_{n+1})$, it
is clearly bounded by $Ch^{2}$, where the constant $C$ is independent of
$\epsilon$. This stems from the polynomial growth of the test function $\phi$
and its derivatives up to order $4$, as well as from the bound
$|Z_{n+1}-Z_{n}|\leq M_{n}(1+|Z_{n}|)\sqrt{h}$.

Performing the Taylor expansion of the expectation of the test function $\phi$
applied to the exact solution leads to

$\mathbb{E}(\phi(Z(t_{n+1}))|Z_{n}=z)=\phi(z)+h(\mathcal{L}_{t_{n}/\epsilon}\phi)(z)+R,$

where, according to [17, P. 26] 222There are three other terms with vanishing
expectations.

$R=\mathbb{E}\left(\int_{t_{n}}^{t_{n+1}}\left(\int_{t_{n}}^{s}({\mathcal{L}}_{\tau/\epsilon}^{2}\phi)(Z(\tau))d\tau\right)ds|Z(t_{n})=z\right).$

Let us notice that the derivative with respect to $\theta$ will not appear in
${\mathcal{L}}_{t_{n}/\epsilon}\phi$ because our test function $\phi$ does not
depend explicitly on time. However, it will appear in
${\mathcal{L}}_{t_{n}/\epsilon}^{2}\phi$ owing to the explicit dependence of
$g_{t/\epsilon}$ and $\Sigma_{t/\epsilon}$ on $t$, and only because of this
derivation with respect to $\theta$, we will have terms of order
$\mathcal{O}(\frac{1}{\epsilon})$ that need to be bounded.

Let $g^{1}_{\theta}$ and $g^{2}_{\theta}$ be the two components of
$g_{\theta}$ (both belong to ${\mathbb{R}}^{d}$), then
$\partial_{\theta}g_{\theta}(z)=(\partial_{\theta}g_{\theta}^{1}(z),\partial_{\theta}g_{\theta}^{2}(z))^{T}=(0,\partial_{\theta}g_{\theta}^{2}(z))^{T}$,
with

$\displaystyle\partial_{\theta}(g_{\theta}^{2}(z))=\partial_{\theta}(g_{\theta}^{2}(\underline{x},y))$
$\displaystyle=\frac{1}{\epsilon}\partial_{\theta}f_{t/\epsilon}(\Phi_{t/\epsilon}(\underline{x})+y)+\frac{1}{\epsilon}f^{\prime}_{t/\epsilon}(\Phi_{t/\epsilon}(\underline{x})+y)\partial_{\theta}\Phi_{t/\epsilon}(\underline{x})$
$\displaystyle-\frac{1}{\epsilon}\partial_{\theta}f_{t/\epsilon}(\underline{x})-\partial_{\theta}F^{\prime}_{t/\epsilon}(\underline{x})\left<f\right>(\underline{x})$
$\displaystyle-\frac{1}{2}\partial_{\theta}F^{\prime\prime}_{t/\epsilon}(\underline{x})(\sigma(\underline{x}),\sigma(\underline{x})).$

Using the same notations for $\Sigma_{\theta}(z)$, we have

$\displaystyle\partial_{\theta}(\Sigma_{t/\epsilon}^{2}(z))=\frac{1}{\epsilon}\sigma^{\prime}(\Phi_{t/\epsilon}(\underline{x})+y)\partial_{\theta}\Phi_{t/\epsilon}(\underline{x})-\partial_{\theta}F^{\prime}_{t/\epsilon}(\underline{x})\sigma(\underline{x}).$

It remains to bound the difference

$\frac{1}{\epsilon}\partial_{\theta}f_{\tau/\epsilon}(\Phi_{\tau/\epsilon}(\underline{X}(\tau))+Y(\tau))-\frac{1}{\epsilon}\partial_{\theta}f_{\tau/\epsilon}(\underline{X}(\tau)),$
(3.12)

since all other terms are bounded independently of $\epsilon$ thanks to the
smoothnes of $f_{\theta}$, $\sigma$ and their derivatives with respect to the
space variable, and using the fact that
$\partial_{\theta}\Phi_{\tau/\epsilon}(X)=\mathcal{O}(\epsilon)$. Now, taking
into account the uniform Lipschitz continuity of
$\partial_{\theta}f_{\theta}$, the polynomial growth of the test function
$\phi$, and the above lemma, we have

$\begin{split}\left|\mathbb{E}\left(\frac{1}{\epsilon}\partial_{\theta}f_{\tau/\epsilon}(\Phi_{\tau/\epsilon}(\underline{X}(\tau))+Y(\tau))-\frac{1}{\epsilon}\partial_{\theta}f_{\tau/\epsilon}(\underline{X}(\tau))\right)\right|&\leq\frac{L}{\epsilon}\mathbb{E}\left|\Phi_{\tau/\epsilon}(\underline{X}(\tau))+Y(\tau)-\underline{X}(\tau)\right|\\\
&\leq\frac{L}{\epsilon}\mathbb{E}\left|\epsilon
F_{\tau/\epsilon}(\underline{X}(\tau))\right|+\frac{L}{\epsilon}\mathbb{E}\left|Y(\tau)\right|\\\
&=L\mathbb{E}\left|F_{\tau/\epsilon}(\underline{X}(\tau))\right|+Lc\\\ &\leq
L(K(1+\mathbb{E}\left|\underline{X}(\tau)\right|)+c)<+\infty\end{split}$

independently of the value of $\epsilon$. Hence, we have

$\left|\mathbb{E}\left(\int_{t_{n}}^{t_{n+1}}\left(\int_{t_{n}}^{s}({\mathcal{L}}_{\tau/\epsilon}^{2}\phi)(Z(\tau))d\tau\right)ds|Z(t_{n})=z\right)\right|\leq
C^{\prime}h^{2},$

and the local order 2 is proved.

### Step $3$

The boundedness of the moments and the local weak order $2$ imply the global
weak convergence of order $1$ by a theorem from [17] (see also [18, chap2.2]).

We have proved, for any test function $\phi\in
C^{4}_{poly}({\mathbb{R}}^{2d},{\mathbb{R}})$, for all $n=0,1,\dots,N$, that

$\left|\mathbb{E}(\phi(Z(t_{n}))-\phi(Z_{n}))\right|\leq Ch,$

where $Z=(\underline{X},Y)$ and $C$ is independent of $n$ and $\epsilon$. We
need to prove that for any test function $\phi\in
C^{4}_{poly}({\mathbb{R}}^{d},{\mathbb{R}})$, for all $n=0,1,\dots,N$,

$\left|\mathbb{E}(\phi(X(t_{n}))-\phi(X_{n}))\right|=\left|\mathbb{E}(\phi(\Phi_{t_{n}/\epsilon}(\underline{X}(t_{n}))+Y(t_{n}))-\phi(\Phi_{t_{n}/\epsilon}(\underline{X}_{n})+Y_{n}))\right|\leq
C_{1}h.$

where $C_{1}$ is independent of $n$ and $\epsilon$. Let $\phi\in
C^{4}_{poly}({\mathbb{R}}^{d},{\mathbb{R}})$, for each fixed parameter
$\theta\in\mathbb{T}$, we consider the test function

$\psi_{\theta}(Z)=\phi(\Phi_{\theta}(\underline{X})+Y).$

Note that $\psi_{\theta}\in C^{4}_{poly}({\mathbb{R}}^{d},{\mathbb{R}})$,
since $\Phi_{\theta}\in C^{4}_{poly}({\mathbb{R}}^{d},{\mathbb{R}}^{d})$ (by
assumption on $f_{\theta}$). Hence, for each $n=0,1,\dots,N$, we have

$\begin{split}\left|\mathbb{E}(\psi_{t_{n}/\epsilon}(Z(t_{n}))-\psi_{t_{n}/\epsilon}(Z_{n}))\right|&=\left|\mathbb{E}(\phi(\Phi_{t_{n}/\epsilon}(\underline{X}(t_{n}))+Y(t_{n}))-\phi(\Phi_{t_{n}/\epsilon}(\underline{X}_{n})+Y_{n}))\right|\\\
&=\left|\mathbb{E}(\phi(X(t_{n}))-\phi(X_{n}))\right|\leq C_{1}h,\end{split}$

where $C_{1}=C$.

##### Strong convergence

For the strong convergence with order 1/2 (for $Z$), in addition to the
bounded moments and the local weak order $2$, using [16], it is sufficient to
show first order strong convergence after one step, i.e,

$\mathbb{E}(|Z(t_{n+1})-Z_{n+1}|^{2};Z(t_{n})=z)^{\frac{1}{2}}\leq Ch,$ (3.13)

where the generic constant $C$ is independent of $n,h$, and $\epsilon$. By
Wagner-Platen expansion [17, chap1.2.2], we have

$\displaystyle Z(t_{n+1})-Z_{n+1}$
$\displaystyle=\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{s}\Lambda_{\tau/\epsilon}\Sigma_{\tau/\epsilon}(Z(\tau))dW(\tau)dW(s)+\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{s}{\mathcal{L}}_{\tau/\epsilon}\Sigma_{\tau/\epsilon}(Z(\tau))d\tau
dW(s)$
$\displaystyle+\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{s}\Lambda_{\tau/\epsilon}g_{\tau/\epsilon}(Z(\tau))dW(\tau)d(s)+\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{s}{\mathcal{L}}_{\tau/\epsilon}g_{\tau/\epsilon}(Z(\tau))d\tau
d(s)$

where
$\Lambda\varphi_{\theta}(z)=\varphi^{\prime}_{\theta}(z)\Sigma_{\theta}(z)$,
${\mathcal{L}}_{\theta}$ is defined in (3.11) and

$(\Lambda_{\theta}\psi)(z)=\psi^{\prime}(z)\Sigma_{\theta}(z).$

The only terms we need to bound are the derivatives with respect to $\theta$
of $g_{\theta}(Z)$ and $\Sigma_{\theta}(Z)$ at $\theta=\tau/\epsilon$ arising
from the application of the differential operator
${\mathcal{L}}_{\tau/\epsilon}$ to $g_{\tau/\epsilon}(Z)$ and
$\Sigma_{\tau/\epsilon}(Z)$ respectively. It can be checked that the other
terms in the integrals are again uniformly bounded by the regularity
assumptions, Lemma 3.3, and the boundedness of the moments of $Z(t)$. Once
this is done, we can conclude that (3.13) is satisfied. The uniform
boundedness of $\partial_{\theta}g_{\tau/\epsilon}(Z)$ was already proved in
the weak convergence case. Now, we have
$\partial_{\theta}\Sigma_{\tau/\epsilon}(Z)=(0,\partial_{\theta}\Sigma_{\theta}^{2}(Z))$,
and

$\partial_{\theta}\Sigma^{2}_{\tau/\epsilon}(Z)=\frac{1}{\epsilon}\sigma^{\prime}(\Phi_{\tau/\epsilon}(\underline{X})+Y)\partial_{\theta}\Phi_{\tau/\epsilon}(\underline{X})-\partial_{\theta}F_{\tau/\epsilon}^{\prime}(\underline{X})\sigma(\underline{X})$

The right hand side is uniformly bounded since
$\partial_{\theta}\Phi_{\theta}(x)=\mathcal{O}(\epsilon)$ (thanks to the
periodicity of $f_{\theta}$ with respect to $\theta$ and the definition of
$\left<f\right>$), $\sigma^{\prime}(x)$ and
$\partial_{\theta}F^{\prime}_{\theta}(x)$ have at most polynomial growth, and
$\sigma$ has linear growth. It follows that, after several applications of the
Itô isometry, the Cauchy-Schwarz and the Young inequalities,

$\displaystyle\mathbb{E}(|Z(t_{n+1})-Z_{n+1}|^{2})\leq Ch^{2},$

and thus $\mathbb{E}(|Z(t_{n+1})-Z_{n+1}|^{2})^{\frac{1}{2}}\leq Ch$.

Now we conclude the local strong order $1$ for $X$. First, note that

$\mathbb{E}(|Z(t_{n+1})-Z_{n+1}|^{2})=\mathbb{E}(|\underline{X}(t_{n+1})-\underline{X}_{n+1}|^{2})+\mathbb{E}(|Y(t_{n+1})-Y_{n+1}|^{2})\leq
Ch^{2},$

which implies that each component of $Z_{n+1}=(\underline{X}_{n+1},Y_{n+1})$
converges strongly with local order $1$ to the corresponding component of
$Z(t_{n+1})=(\underline{X}(t_{n+1}),Y(t_{n+1}))$. Next, for a given
$X(t_{n})=X_{n}$, we have

$\begin{split}&\mathbb{E}\left(\left|X(t_{n+1})-X_{n+1}\right|^{2}\right)=\mathbb{E}\left(\left|\Phi_{t_{n+1}/\epsilon}(\underline{X}(t_{n+1}))+Y(t_{n+1})-\Phi_{t_{n+1}/\epsilon}(\underline{X}_{n+1})-Y_{n+1}\right|^{2}\right)\\\
&\leq\mathbb{E}\left(\left|\Phi_{t_{n+1}/\epsilon}(\underline{X}(t_{n+1}))-\Phi_{t_{n+1}/\epsilon}(\underline{X}_{n+1})\right|^{2}\right)+\mathbb{E}\left(\left|Y(t_{n+1})-Y_{n+1}\right|^{2}\right)\\\
&+2\mathbb{E}\left(\left|\left(\Phi_{t_{n+1}/\epsilon}(\underline{X}(t_{n+1}))-\Phi_{t_{n+1}/\epsilon}(\underline{X}_{n+1})\right)\left(Y(t_{n+1})-Y_{n+1}\right)\right|\right)\\\
&\leq
L^{2}\mathbb{E}\left(\left|\underline{X}(t_{n+1})-\underline{X}_{n+1}\right|^{2}\right)+\mathbb{E}\left(\left|Y(t_{n+1})-Y_{n+1}\right|^{2}\right)\\\
&+2\mathbb{E}\left(\left|\Phi_{t_{n+1}/\epsilon}(\underline{X}(t_{n+1}))-\Phi_{t_{n+1}/\epsilon}(\underline{X}_{n+1})\right|^{2}\right)^{\frac{1}{2}}\mathbb{E}\left(\left|Y(t_{n+1})-Y_{n+1}\right|^{2}\right)^{\frac{1}{2}}\\\
&\leq Ch^{2}+Ch^{2}+2(Ch)(Ch)\leq Ch^{2}.\end{split}$

Therefore,

$\mathbb{E}\left(\left|X(t_{n+1})-X_{n+1}\right|^{2}\right)^{\frac{1}{2}}\leq
Ch$

and the local strong order 1 is proved. ∎

######  Remark 3.4.

When applying the Euler-Maruyama method directly to (1.1), a term that
contains $\frac{1}{\epsilon}\partial_{\theta}f_{t/\epsilon}(X(t_{n}))$ will
appear in the remainder. This term causes the local truncation error to be of
order $\mathcal{O}(\frac{1}{\epsilon})$ which causes instability in the high
frequency regime ($\epsilon\ll 1$) even when Assumption 1.2 is satisfied. In
contrast, applying Euler Maruyama method to the micro-macro system (3.5,3.6)
leads instead of $\frac{1}{\epsilon}\partial_{\theta}f_{t/\epsilon}(X(t_{n}))$
to the term (3.12) which can be bounded independently of $\epsilon$ using
Assumption 1.2.

##  4\. Numerical experiments

Throughout this section (except for the last experiment) we will focus on the
Hénon-Heiles model (see [7, 14]). We consider the Hamiltonian

$H(p,q)=\frac{p_{1}^{2}}{2\epsilon}+\frac{p_{2}^{2}}{2}+\frac{q_{1}^{2}}{2\epsilon}+\frac{q_{2}^{2}}{2}+q_{1}^{2}q_{2}-\frac{1}{3}q_{2}^{3}.$

Let

$\begin{cases}X_{1}(t)=\cos\left(\frac{t}{\epsilon}\right)q_{1}(t)-\sin\left(\frac{t}{\epsilon}\right)p_{1}(t),\\\
X_{2}(t)=q_{2}(t),\\\
X_{3}(t)=\sin\left(\frac{t}{\epsilon}\right)q_{1}(t)+\cos\left(\frac{t}{\epsilon}\right)p_{1}(t),\\\
X_{4}(t)=p_{2}(t).\end{cases}$

It can be checked that the variable $X(t)$ satisfies the following ODE [9]

$\frac{dX}{dt}(t)=f_{t/\epsilon}(X(t)),$

with,

$\begin{cases}f^{1}_{\theta}(X)=2\sin\theta(X_{1}\cos\theta+X_{3}\sin\theta)X_{2},\\\
f^{2}_{\theta}(X)=X_{4},\\\
f^{3}_{\theta}(X)=-2\cos\theta(X_{1}\cos\theta+X_{3}\sin\theta)X_{2},\\\
f^{4}_{\theta}(X)=-2(X_{1}\cos\theta+X_{3}\sin\theta)^{2}+X_{2}^{2}-X_{2}.\end{cases}$

Now, we consider the SDE

$dX=f_{t/\epsilon}(X)dt+\sigma(X)dW(t).$ (4.1)

In all our experiments $\epsilon=2^{-2i},~{}i=2,3,4,5$, final time $T=1$; $M$
denotes the number of computed samples and $\Delta t$ denotes the time step
size.

###  4.1. Weak convergence

In this section we use $X^{0}=(0.7,0.7,0.7,0.7)$, $\Delta
t=2^{-i},~{}i=1,\dots,5$, $M=10^{4}$.

####  4.1.1. Multiplicative noise

We consider the above SDE (4.1) with multiplicative noise where
$\sigma(X)=0.2(0,0,X_{1},X_{2})^{T}$. We use the test function $\phi(X)=X_{1}$
to measure the weak convergence. In Figure 1A we plot the weak error with
respect to the time step for different values of $\epsilon$ (left figure)
using the micro-macro method (3.7)-(3.8). We can see that the convergence
behavior looks almost the same, with weak order one, for all the different
values of $\epsilon$. The right picture of Figure 1A shows that for a fixed
time step, the weak error remains almost constant when varying $\epsilon$. The
above description applies also to the integral scheme (see Figure 1B).

(a) Micro-Macro scheme.

(b) Integral scheme.

Figure 1. Weak convergence with multiplicative noise for the micro-macro
method (3.7)-(3.8) and the integral scheme(2.1) .

####  4.1.2. Additive noise

The weak error of the micro-macro scheme (3.7)-(3.8) applied to the SDE (4.1)
with additive noise is shown in Figure 2. We set
$\sigma(X)=(0,0,0.2,0.2)^{T}$, and we perform the test with two different test
functions. We see again the uniform weak order one.

(a) $\phi(X)=X_{1}$.

(b) $\phi(X)=\sum_{i=1}^{4}X_{i}$.

Figure 2. Weak convergence with additive noise for the micro-macro method
(3.7)-(3.8).

###  4.2. Strong convergence

In this section we use $X^{0}=(0.12,0.12,0.12,0.12)$, $\Delta
t=2^{-i},~{}i=4,\dots,8$, $M=10^{2}$.

####  4.2.1. Multiplicative noise

We consider the above SDE (4.1) with multiplicative noise where
$\sigma(X)=0.5X$. Figure 3 shows the uniform strong order $\frac{1}{2}$ for
both methods.

(a) Micro-Macro.

(b) Integral scheme.

Figure 3. Strong convergence with multiplicative noise for methods (3.7)-(3.8)
and (2.1).

####  4.2.2. Additive noise

We consider the above SDE (4.1) with additive noise where
$\sigma(X)=(0,0,0.5,0.5)^{T}$. In addition to the uniform convergence, Figure
4 shows strong order one for the micro-macro method (3.7)-(3.8) since when the
noise is additive, Euler-Maruyama method coincides with Milstein method of
strong order one. This applies to uniformly accurate methods too.

Figure 4. Strong convergence with additive noise for the micro-macro method
(3.7)-(3.8).

###  4.3. Inefficiency of Euler-Maruyama method for particular time steps

Although Euler-Maruyama method seems to work quite well, it still fails for
some particular choices of time steps, while the Micro-Macro method
(3.7)-(3.8) does not. We recall that the importance of uniformly accurate
methods appears more when using higher order schemes. See Figure 5. We
consider the logistic SDE

$dX=\left(X(1-X)+\sin\frac{t}{\epsilon}\right)dt+0.2XdW(t),\qquad X(0)=2.$

We set $\epsilon=0.1$ and we plot in Figure 5A the reference solution (in
blue) calculated with very small time step using the integral scheme and the
solution obtained using EM with time step $\Delta t=0.99(2\pi\epsilon)$. In
Figure 5B, we plot the reference solution (in blue) calculated with very small
time step using the integral scheme and the solution obtained using the
uniformly accurate method (3.7)-(3.8) with time step $h=0.99(2\pi\epsilon)$.

(a) Euler Maruyama

(b) Uniformly accurate method (3.7, 3.8)

Figure 5. Failure of EM for particular time steps.

##  5\. Conclusion

In this work, we have introduced two uniformly accurate methods for solving
numerically stochastic differential equations with oscillatory drift. The
first one is the so-called integral scheme (2.1) and can be derived quite
straightforwardly, whereas the second one is obtained through a more elaborate
transformation, namely a micro-macro decomposition (3.5, 3.6). Both schemes
exhibit weak-order $1$ and strong-order $1/2$, as proved in the corresponding
sections and confirmed numerically in Section 4. Given their comparable
performance, the first scheme is arguably better for its simplicity. However,
it is our belief that the micro-macro scheme exposed here could be generalized
to higher order methods which would be the stochastic counterpart of existing
deterministic uniformly accurate methods [6, 9].

##### Acknowledgement.

The authors would like to thank Gilles Vilmart for useful discussions and
comments.

## References

  * [1] A. Abdulle, I. Almuslimani, and G. Vilmart.  Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations.  SIAM/ASA J. Uncertain. Quantif., 6(2):937–964, 2018.
  * [2] A. Abdulle and T. Li.  S-ROCK methods for stiff Ito SDEs.  Commun. Math. Sci., 6(4):845–868, 2008.
  * [3] A. Abdulle and G. Vilmart.  PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise.  J. Comput. Phys., 242:869–888, 2013.
  * [4] A. Abdulle, G. Vilmart, and K. C. Zygalakis.  Weak second order explicit stabilized methods for stiff stochastic differential equations.  SIAM J. Sci. Comput., 35(4):A1792–A1814, 2013.
  * [5] W. Bao and X. Dong.  Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime.  Numerische Mathematik, 120:189–229, 2012.
  * [6] S. Baumstark, E. Faou, and K. Schratz.  Uniformly accurate exponential-type integrators for Klein-Gordon equations with asymptotic convergence to the classical NLS splitting.  Math. Comp., 87:1227–1254, 2018.
  * [7] P. M. Burrage and K. Burrage.  Structure-preserving Runge-Kutta methods for stochastic Hamiltonian equations with additive noise.  Numer. Algorithms, 65(3):519–532, 2014.
  * [8] F. Castella, P. Chartier, F. Méhats, and A. Murua.  Stroboscopic averaging for the nonlinear schrödinger equation.  Found. Comput. Math., 15:519–559, 2015.
  * [9] P. Chartier, M. Lemou, F. Méhats, and G. Vilmart.  A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations.  Found. Comput. Math., 20(1):1–33, 2020.
  * [10] P. Chartier, J. Makazaga, A. Murua, and G. Vilmart.  Multi-revolution composition methods for highly-oscillatory differential equations.  Numerische Mathematik, 128:167–192, 2014.
  * [11] N. Crouseilles, M. Lemou, and F. Méhats.  Asymptotic preserving schemes for highly oscillatory kinetic equations.  J. Comp. Phys., 248:287–308, 2013.
  * [12] E. Faou and K. Schratz.  Asymptotic preserving schemes for the klein-gordon equation in the non-relativistic limit regime.  Numerische Mathematik, 126:441–469, 2014.
  * [13] I. I. Gīhman and A. V. Skorohod.  Stochastic differential equations.  Springer-Verlag, New York-Heidelberg, 1972.  Translated from the Russian by Kenneth Wickwire, Ergebnisse der Mathematik und ihrer Grenzgebiete, Band 72.
  * [14] M. Han, Q. Ma, and X. Ding.  High-order stochastic symplectic partitioned Runge-Kutta methods for stochastic Hamiltonian systems with additive noise.  Appl. Math. Comput., 346:575–593, 2019.
  * [15] P. Lochak and C. Meunier.  Multiphase Averaging for Classical Systems, volume 72 of Applied Mathematical Sciences.  Springer New York, NY, first edition, 1988.
  * [16] G. N. Mil′shteĭn.  A theorem on the order of convergence of mean-square approximations of solutions of systems of stochastic differential equations.  Teor. Veroyatnost. i Primenen., 32(4):809–811, 1987.
  * [17] G. Milstein.  Numerical integration of stochastic differential equations, volume 313 of Mathematics and its Applications.  Kluwer Academic Publishers Group, Dordrecht, 1995.  Translated and revised from the 1988 Russian original.
  * [18] G. N. Milstein and M. V. Tretyakov.  Stochastic numerics for mathematical physics.  Scientific Computation. Springer-Verlag, Berlin, 2004.
  * [19] J. A. Sanders, F. Verhulst, and J. Murdock.  Averaging Methods in Nonlinear Dynamical Systems, volume 59 of Applied Mathematical Sciences.  Springer New York, NY, second edition, 2007.
  * [20] D. Talay.  Discrétisation d’une équation différentielle stochastique et calcul approché d’espérances de fonctionnelles de la solution.  ESAIM: Mathematical Modelling and Numerical Analysis - Modélisation Mathématique et Analyse Numérique, 20(1):141–179, 1986.
  * [21] G. Vilmart.  Weak second order multirevolution composition methods for highly oscillatory stochastic differential equations with additive or multiplicative noise.  SIAM J. Sci. Comput., 36:1770––1796, 2014.