# Actuarial strategy for pricing Asian options under a mixed fractional
Brownian motion with jumps

Foad Shokrollahi  Department of Mathematics and Statistics, University of
Vaasa, P.O. Box 700, FIN-65101 Vaasa, FINLAND<EMAIL_ADDRESS>,
Davood Ahmadian  Faculty of Mathematical Sciences, University of Tabriz
29Blvd, Tabriz, Iran<EMAIL_ADDRESS>and  Luca Vincenzo Ballestra
Department of Statistical Sciences, Alma Mater Studiorum University of
Bologna, Via delle Belle Arti 41, 40126 Bologna, Italy
<EMAIL_ADDRESS>
(Date: August 28, 2024)

###### Abstract.

The mixed fractional Brownian motion ($mfBm$) has become quite popular in
finance, since it allows one to model long-range dependence and self-
similarity while remaining, for certain values of the Hurst parameter,
arbitrage-free. In the present paper, we propose approximate closed-form
solutions for pricing arithmetic Asian options on an underlying described by
the $mfBm$. Specifically, we consider both arithmetic Asian options and
arithmetic Asian power options, and we obtain analytical formulas for pricing
them based on a convenient approximation of the strike price. Both the
standard $mfBm$ and the $mfBm$ with Poisson log-normally distributed jumps are
taken into account.

###### Key words and phrases:

Mixed fractional Brownian motion; Asian options; Power options; Jump-
diffusion; Option pricing

###### 2010 Mathematics Subject Classification:

91G20; 91G80; 60G22

##  1\. Introduction

As documented by several empirical studies, returns of financial assets are
often affected by long-range dependence and self-similarity, and thus several
scholars have proposed to model them by using the so-called fractional
Brownian motion ($fBm$) with Hurst exponent $H\in(1/2,1)$, see, e.g., [9],
[12], [17]. Nevertheless, as shown in [3], and [7], the $fBm$ is not
arbutrage-free, and thus, some authors [6], [13], [14], [28] have proposed to
use the so-called mixed fractional Brownian motion ($mfBm$), i.e., a linear
combination of the Wiener process and the $fBm$ itself. When the Hurst
exponent of the $fBm$ is greater than $1/2$, the $mfBm$ results in a Gaussian
long-memory process suitable for modeling the returns of financial assets.
Furthermore, as shown in [6], if $H$ is greater than $3/4$ the $mfBm$ is
arbitrage-free. Therefore, the $mfBm$ has been largely employed for pricing
various kinds of derivatives, including stock options [2], [5], [25], equity
warrants [26], currency options [22, 20, 19], and credit derivatives [10].

The $mfBm$ has also been used for pricing Asian options, i.e. options whose
payoff depends on the average, either geometric or arithmetic, of the
underlying asset over the whole option’s lifespan [18]. In particular, [16]
obtained an exact closed-form solution for geometric Asian options, [27]
derived an exact analytical formula for geometric Asian power options, [25]
proposed an exact closed-form solution for geometric Asian rainbow options on
two underlying assets, which was generalized to the case on an arbitrarily
large number of underlying assets in [1]. Furthermore, [15] considered a
$mfBm$ with jumps and obtained an exact analytical formula for valuing
geometric Asian power options.

All the exact closed-form solutions proposed in the aforementioned papers are
valid for Asian options of geometric type. By contrast, for arithemtic Asian
options, exact closed-form solutions cannot be found, and thus some numerical
approximation is required. To this aim, some researchers have proposed to
apply the Monte Carlo simulation method, which could also be enhanced by using
the price of a geometric Asian option as a control variate (see [15]).
Nevertheless, this approach, albeit flexible and relatively simple to
implement, has the disadvantage of being very slow to converge.

Therefore, in this paper we propose the use of approximate closed-form
solutions. Specifically, we consider both arithmetic Asian options and
arithmetic Asian power options, and we value them by applying an analytical
formula that we obtain by a convenient approximation of the strike price. Both
the standard $mfBm$ and the $mfBm$ with Poisson log-normally distributed jumps
are taken into account.

The remainder of the paper is organized as follows: in Section 2 we describe
the $mfBm$ with jumps and present exact analytical formulas for pricing Asian
options of geometric type. In Section 3 we derive an approximate closed-form
solution for pricing arithmetic Asian options, which in Section 4 we extend to
arithmetic Asian power options. Finally, Section 5 concludes.

##  2\. Preliminaries

We assume that the following assumptions hold:

  1. (i)

the dynamics of the stock price is governed by the following equation

(2.1) $\displaystyle S_{t}=S_{0}\exp\\{f(t)+\sigma B_{t}+\varepsilon
B_{t}^{H}-\frac{1}{2}\sigma^{2}t-\frac{1}{2}\varepsilon^{2}t^{2H}+\sum_{i=1}^{N_{t}}J_{i}\\},$

where $f(t)$ is a non-random function of $t$, and, $\sigma$, $\varepsilon$ are
constants; $B_{t}$ and $B_{t}^{H}$ are a standard Brownian motion and a $fBm$
respectively; $N_{t}$ is a poisson process with rate $\lambda$; $e^{J_{i}}-1$
is jump size which is a sequence of independent identically distributed and
$J_{i}\sim N(\mu_{J},\sigma_{J}^{2})$. Moreover, all three sources of
randomness, the $B_{t}$ and $B_{t}^{H}$ and $N_{t}$ are supposed to be
independent;

  2. (ii)

the risk free interest rate $r$ and dividend rate $q$ are known and constant
through time;

  3. (iii)

the option can be exercised only at the maturity time;

  4. (iv)

there are no transaction costs in buying or selling the stocks or option.

Based on option valuation strategy and risk neutral measure (cf.Refs.[8, 21]),
the dynamics of the stock price process $S(t)$ can be written as:

(2.2) $\displaystyle S_{t}=S_{0}\exp[(r-q)t+\sigma B_{t}+\varepsilon
B_{t}^{H}-\frac{1}{2}\sigma^{2}t-\frac{1}{2}\varepsilon^{2}t^{2H}+\sum_{i=1}^{N_{t}}J_{i}].$

######  Lemma 2.1.

(2.3) $\displaystyle E[S_{T}]=S_{0}e^{(r-q)T+\lambda(\rho-1)T}.$

where $\rho=E[e^{J_{1}}]$.

###### Proof.

(2.4) $\displaystyle E[S_{T}]$ $\displaystyle=$ $\displaystyle
S_{0}E[e^{(r-q)t+\sigma B_{t}+\varepsilon
B_{t}^{H}-\frac{1}{2}\sigma^{2}t-\frac{1}{2}\varepsilon^{2}t^{2H}+\sum_{i=1}^{N_{T}}J_{i}}$
$\displaystyle=$ $\displaystyle S_{0}E[e^{(r-q)t+\sigma B_{t}+\varepsilon
B_{t}^{H}-\frac{1}{2}\sigma^{2}t-\frac{1}{2}\varepsilon^{2}t^{2H}}]E[\prod_{i=1}^{N_{T}}e^{J_{i}}]$
$\displaystyle=$ $\displaystyle
S_{0}e^{(r-q)T}E[\prod_{i=0}^{N_{T}}e^{J_{i}}],$

and

(2.5) $\displaystyle E\left[\prod_{i=0}^{N_{T}}e^{J_{i}}\right]$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda T}}{n!}E[\prod_{i=1}^{n}e^{J_{i}}]$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\lambda
T}}{n!}\rho^{n}$ $\displaystyle=$ $\displaystyle
e^{\lambda(\rho-1)T}\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\rho\lambda
T}}{n!}$ $\displaystyle=$ $\displaystyle e^{\lambda(\rho-1)T}.$

Then,

(2.6) $\displaystyle E[S_{T}]=S_{0}e^{(r-q)T+\lambda(\rho-1)T}.$

∎

This section deals with the new pricing model for the currency options using
the actuarial approach, when the spot exchange rate follows the $MFBM$ with
jumps process. This model can be applied to different financial markets such
as: in the arbitrage-free, equilibrium and complete markets and also in the
arbitrage,non-equilibrium and incomplete markets.

######  Definition 2.1.

([4]) The expectation return rate $\theta(t)$ of $S_{t}$ on $t\in[0,T]$ is
defined to $\int_{0}^{T}\theta(s)ds$ as follows:

(2.7)
$\displaystyle\frac{E\Big{(}S(T)\Big{)}}{S_{0}}=\exp\Big{(}\int_{0}^{T}\theta(s)ds\Big{)},$

where $\theta(t)$ is the continuously compounding ram of return of $S_{t}$ at
time $t$.

######  Definition 2.2.

Let $A(T)$ be the average price of the underlying asset over the predetermined
interval. Then, the value of geometric Asian options with exercise price $K$
and maturity time $T$ by actuarial approach are as follows:

(2.8) $\displaystyle
C(S_{T},K)=E\Big{[}\Big{(}\exp\Big{(}-\int_{0}^{T}\theta(t)dt\Big{)}A(T)-Ke^{-rT}\Big{)}I_{B}\Big{]},$

for Asian call option, where
$B=\exp\Big{(}-\int_{0}^{T}\theta(t)dt\Big{)}A(T)>Ke^{-rT}$.

(2.9) $\displaystyle
P((S_{T},K)=E\Big{[}\Big{(}Ke^{-rT}-\exp\Big{(}-\int_{0}^{T}\theta(t)dt\Big{)}A(T)\Big{)}I_{B^{\prime}}\Big{]},$

for Asian put option,
$B^{\prime}=\exp\Big{(}-\int_{0}^{T}\theta(t)dt\Big{)}A(T)<Ke^{-rT}$.

######  Theorem 2.1.

Let the stock price $S_{T}$ satisfy Equation (2.2). Then, the value of
geometric Asian options with exercise price $K$ is given by

(2.10) $\displaystyle C(S_{T},K)$ $\displaystyle=$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\lambda
T}}{n!}\left[e^{-(r-q)T-\lambda(\rho-1)T+\hat{\mu}+\frac{1}{2}\hat{\sigma}^{2}}\Phi(d_{1})-Ke^{-rT}\Phi(d_{2})\right],$

for the call options, and

(2.11) $\displaystyle P(S_{T},K)$ $\displaystyle=$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\lambda
T}}{n!}\left[Ke^{-rT}\Phi(-d_{2})-\mathrm{e}^{-(r-q)T-\lambda(\rho-1)T+\hat{\mu}+\frac{1}{2}\hat{\sigma}^{2}}\Phi(-d_{1})\right],$

for the put options, where

$\displaystyle d_{2}=\frac{\hat{\mu}-\ln
U}{\hat{\sigma}},\,d_{1}=d_{2}+\hat{\sigma},\,U=Ke^{-qT-(1-\rho)\lambda T},$
$\displaystyle\rho=E(e^{J})=e^{\mu_{J}+\frac{1}{2}\sigma_{J}^{2}},\,\hat{\mu}=\ln
S_{0}+\frac{1}{2}(r_{n}-\frac{1}{2}\sigma_{n}^{2})T,\,\hat{\sigma}^{2}=\sigma_{n}^{2}\frac{T}{3},$
(2.12) $\displaystyle
r_{n}=r-q+\frac{n}{T}(\mu_{J}+\frac{1}{2}\sigma_{J}^{2}),\,\sigma_{n}^{2}=\sigma^{2}+\frac{1}{T}(\varepsilon^{2}T^{2H}+n\sigma_{J}^{2}).$

and $\Phi(.)$ is cumulative normal density function.

###### Proof.

We know that $B_{T}\sim N(0,T)$ and $B_{T}^{H}\sim N(0,T^{2H})$. If $Z_{1}\sim
N(0,1)$ and $Z_{2}\sim N(0,1)$, then $B_{T}=\sqrt{T}Z_{1}$ and
$B_{T}^{H}=\sqrt{T^{2H}}Z_{2}$. Given that $J_{i}\sim
N(\mu_{J},\sigma_{J}^{2})$ and $e^{J_{i}}-1$ is a sequence of independent
identically distributed, we have $\sum_{i=1}^{n}J_{i}\sim
N\left(n\mu_{J},n\sigma_{J}^{2}\right)$; and for $Z_{3}\sim N(0,1)$ we have
$\sum_{i=1}^{n}J_{i}=n\mu_{J}+\sqrt{n}\sigma_{J}Z_{3}$. Then for
$Z_{n}\sim(0,1)$, it follows

(2.13)
$\displaystyle\sigma_{n}\sqrt{T}Z_{n}=\sigma\sqrt{T}Z_{1}+\varepsilon\sqrt{T^{2H}}Z_{2}+\sqrt{n}\sigma_{J}Z_{3},$

where
$\sigma_{n}^{2}=\sigma^{2}+\frac{1}{T}\left(\varepsilon^{2}T^{2H}+n\sigma_{J}^{2}\right)$.
Consequently, assuming $N_{T}=n$, from Eqs. (2.2) and (2.13) we obtain

(2.14) $\displaystyle S_{T}^{n}$ $\displaystyle=$ $\displaystyle
S_{0}\exp\left[(r-q)T-\frac{1}{2}\sigma^{2}T-\frac{1}{2}\varepsilon^{2}T^{2H}+n\mu_{J}+\sigma_{n}\sqrt{T}Z_{n}\right]$
$\displaystyle=$ $\displaystyle
S_{0}\exp\left[(r_{n}-\frac{1}{2}\sigma_{n}^{2})T+\sigma_{n}\sqrt{T}Z_{n}\right],$

here

(2.15) $\displaystyle
r_{n}=r-q+\frac{n}{T}(\mu_{J}+\frac{1}{2}\sigma_{J}^{2}).$

Let,

(2.16) $\displaystyle L(T)=\frac{1}{T}\int_{0}^{T}\ln S_{t}dt,$

and

(2.17) $\displaystyle G(T)=\exp\left\\{L(T)\right\\}.$

We know that the random variable $L(T)$ has Gaussian distribution under the
risk-neutral probability measure. Let $\hat{\mu}$ and $\hat{\sigma}^{2}$
denote the mean and the variance of the random variable $L(T)$, respectively.
If $N_{T}=n$, then

(2.18) $\displaystyle\hat{\mu}$ $\displaystyle=$ $\displaystyle
E[L(T)]=\frac{1}{T}\int_{0}^{T}E\left[\ln S_{t}^{n}\right]dt=\ln
S_{0}+\frac{1}{T}\int_{0}^{T}(r_{n}-\frac{1}{2}\sigma_{n}^{2})tdt$
$\displaystyle=$ $\displaystyle\ln
S_{0}+\frac{1}{2}(r_{n}-\frac{1}{2}\sigma_{n}^{2})T,$

and

(2.19) $\displaystyle\hat{\sigma}^{2}$ $\displaystyle=$ $\displaystyle
Var[L(T)]=E[L(T)-\hat{\mu}]^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{T^{2}}\sigma_{n}^{2}\int_{0}^{T}\int_{0}^{T}E\left(\sqrt{t}Z_{n}(t)\sqrt{s}Z_{n}(s)\right)dtds$
$\displaystyle=$
$\displaystyle\frac{1}{T^{2}}\sigma_{n}^{2}\int_{0}^{T}\int_{0}^{T}[|t|+|s|-|t-s|]dtds$
$\displaystyle=$ $\displaystyle\sigma_{n}^{2}\frac{T}{3}$

The random variable $G(T)$ is log-normally distributed and the random variable
$\ln G(T)$ has the Gaussian distribution with the mean $\hat{\mu}$ and the
variance $\hat{\sigma}^{2}$ as obtained above. Then, from Eq. (2.8), the price
of geometric Asian call option is:

(2.20) $\displaystyle C(S_{T},K)$ $\displaystyle=$ $\displaystyle
E\left[\left(\exp\left(-\int_{0}^{T}\theta(t)dt\right)G(T)-Ke^{-rT}\right)I_{B}\right]$
$\displaystyle=$ $\displaystyle
E\left[\left(e^{-(r-q)T-\lambda(\rho-1)T}G(T)-Ke^{-rT}\right)I_{B}\right]$
$\displaystyle=$ $\displaystyle
e^{-(r-q)T-\lambda(\rho-1)T}E\left[G(T)I_{B}\right]-Ke^{-rT}E\left[I_{B}\right]$
$\displaystyle=$ $\displaystyle
e^{-(r-q)T-\lambda(\rho-1)T}I_{1}-Ke^{-rT}I_{2},$

where $I_{1}=E\left[G(T)I_{B}\right]$ and $I_{2}=E\left[I_{B}\right]$. Let

(2.21) $\displaystyle B$ $\displaystyle=$
$\displaystyle\left\\{x:e^{x}>U=Ke^{-qT-(1-\rho)\lambda T}\right\\}$
$\displaystyle=$ $\displaystyle\left\\{x:G(T)>U\right\\}$ $\displaystyle=$
$\displaystyle\left\\{y:e^{\hat{\mu}+\hat{\sigma}y}>U\right\\}=\left\\{y:\hat{\mu}+\hat{\sigma}y>\ln
U\right\\}$ $\displaystyle=$ $\displaystyle\left\\{y:y>-d_{2}=\frac{\ln
U-\hat{\mu}}{\hat{\sigma}}\right\\}.$

Observe that

(2.22) $\displaystyle I_{2}$ $\displaystyle=$ $\displaystyle
E\left[I_{B}\right]=\sum_{n=0}^{\infty}E\left[I_{G(T)>-d_{2}}\right]P\left[N_{T}=n\right]$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda T}}{n!}E\left[I_{G(T)>-d_{2}}\right]$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\lambda
T}}{n!}\int_{-d_{2}}^{\infty}\phi(y)dy$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\lambda
T}}{n!}\Phi(d_{2}),$

and

(2.23) $\displaystyle I_{1}$ $\displaystyle=$ $\displaystyle
E\left[G(T)I_{B}\right]=\sum_{n=0}^{\infty}E\left[G(T)I_{G(T)>-d_{2}}\right]P\left[N_{T}=n\right]$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda T}}{n!}E\left[e^{L(T)}I_{G(T)>-d_{2}}\right]$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda
T}}{n!}\int_{B}e^{x}\frac{1}{\sqrt{2\pi}\hat{\sigma}}e^{-\frac{(x-\hat{\mu})^{2}}{2\hat{\sigma}^{2}}}dx$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda
T}}{n!}\int_{B}e^{\hat{\mu}+\hat{\sigma}y}\frac{1}{\sqrt{2\pi}\hat{\sigma}}e^{-\frac{(x-\hat{\mu})^{2}}{2\hat{\sigma}^{2}}}dy$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda
T}}{n!}e^{\hat{\mu}+\frac{1}{2}\hat{\sigma}^{2}}\int_{-d_{2}}^{\infty}\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}(y-\hat{\sigma})^{2}}dy$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda
T}}{n!}e^{\hat{\mu}+\frac{1}{2}\hat{\sigma}^{2}}\int_{-d_{2}-\hat{\sigma}}^{\infty}\phi(y)dy$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda
T}}{n!}e^{\hat{\mu}+\frac{1}{2}\hat{\sigma}^{2}}\Phi(d_{1}=d_{2}+\hat{\sigma}).$

From Eqs. (2.23), (2.22) and (2.20), we get desired result. ∎

Figure 1. The variation of the geometric Asian option with respect to the
Husrt parameter $H$.

Figure 2. The variation of the geometric Asian option with respect to the
$\sigma_{J}$.

Figure 3. The variation of the geometric Asian option with respect to the
$\mu_{J}$.

In Figs 1-3 we show the price of the geometric Asian option computed using
formula (2.10) as a function of the Hurst and the jump parameters.

Figure 4. The variation of the geometric Asian option with respect to the
strike price $K$ and time to maturity (year) $T$.

In Fig 4 we show the price of the geometric Asian option computed using
formula (2.10) as a function of the strike price and of the maturity.

##  3\. Arithmetic Asian option

######  Lemma 3.1.

([24])

An approximate value of the arithmetic Asian option is given by

(3.1) $\displaystyle\tilde{C}(S_{T},K)$ $\displaystyle=$ $\displaystyle
e^{-rT}E\left[(G(T)-K^{\prime})^{+}\right],$

where $G(T)$ ia the geometric Asian option and $K^{\prime}$ is the adjusted
strike price as follows

(3.2) $\displaystyle K^{\prime}=K-\left(E[A]-E[G]\right).$

Now, based on the actuarial approach an approximate of the arithmetic Asian
options are proposed [23, 11].

######  Theorem 3.1.

An approximate value of the arithmetic Asian option is given by

(3.3) $\displaystyle\tilde{C}(S_{T},K)$ $\displaystyle=$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\lambda
T}}{n!}\left[e^{-(r-q)T-\lambda(\rho-1)T+\hat{\mu}+\frac{1}{2}\hat{\sigma}^{2}}\Phi(d_{1})-K^{\prime}e^{-rT}\Phi(d_{2})\right],$

for the call options, and

(3.4) $\displaystyle\tilde{P}(S_{T},K)$ $\displaystyle=$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\lambda
T}}{n!}\left[K^{\prime}e^{-rT}\Phi(-d_{2})-\mathrm{e}^{-(r-q)T-\lambda(\rho-1)T+\hat{\mu}+\frac{1}{2}\hat{\sigma}^{2}}\Phi(-d_{1})\right],$

for the put options, where

$\displaystyle d_{2}=\frac{\hat{\mu}-\ln
U}{\hat{\sigma}},\,d_{1}=d_{2}+\hat{\sigma},$
$\displaystyle\rho=E(e^{J})=e^{\mu_{J}+\frac{1}{2}\sigma_{J}^{2}},\,U=K^{\prime}e^{-qT-(1-\rho)\lambda
T},$ $\displaystyle\hat{\mu}=\ln
S_{0}+\frac{1}{2}(r_{n}-\frac{1}{2}\sigma_{n}^{2})T,\,\bar{\mu}=\frac{S_{0}}{r_{n}T}\left(e^{r_{n}T}-1\right)$
$\displaystyle\hat{\sigma}^{2}=\sigma_{n}^{2}\frac{T}{3},\,\sigma_{n}^{2}=\sigma^{2}+\frac{1}{T}(\varepsilon^{2}T^{2H}+n\sigma_{J}^{2})$
(3.5) $\displaystyle
K^{\prime}=K+E[G]-\bar{\mu},\,r_{n}=r-q+\frac{n}{T}(\mu_{J}+\frac{1}{2}\sigma_{J}^{2}),\leavevmode\nobreak\
E[G]=e^{\hat{\mu}+\frac{1}{2}\hat{\sigma}^{2}}.$

and $\Phi(.)$ is cumulative normal density function.

###### Proof.

Suppose

(3.6) $\displaystyle A(T)=\frac{1}{T}\int_{0}^{T}S_{t}dt$

be an arithmetic average with mean $\bar{\mu}=E\left[A(T)\right]$, then from
Lemma 2.1 we can get

(3.7) $\displaystyle\bar{\mu}$ $\displaystyle=$ $\displaystyle
E\left[A(T)\right]=\frac{1}{T}\int_{0}^{T}E\left[S_{T}\right]dt$
$\displaystyle=$ $\displaystyle\frac{S_{0}}{T}\int_{0}^{T}e^{{r_{n}t}}dt$
$\displaystyle=$ $\displaystyle\frac{S_{0}}{r_{n}T}\left(e^{r_{n}T}-1\right).$

From Theorem 3.1, to provide an approximation of the arithmetic Asian options,
we need to calculate the new exercise price $K^{\prime}$,

(3.8) $\displaystyle K^{\prime}$ $\displaystyle=$ $\displaystyle
K-\left(E[A]-E[G]\right)$ $\displaystyle=$ $\displaystyle K+E[G]-\bar{\mu},$

where $E[G]=e^{\hat{\mu}+\frac{1}{2}\hat{\sigma}^{2}}$ is calculated in
Theorem 2.1. Now, by replacing $K$ with $K^{\prime}$ in the proof of Theorem
2.1, the proof is completed. ∎

Suppose, $G$ and $A$ be the geometric and arithmetic average, respectively.
Since the geometric average is less than the arithmetic average, $G<A$, the
geometric option price is a lower bound for the arithmetic option price

(3.9) $\displaystyle C_{G}=e^{-rT}E\left[(G-K)^{+}\right]\leq
e^{-rT}E\left[(A-K)^{+}\right]=C_{A}.$

The inequality

(3.10) $\displaystyle\max(A-K,0)$ $\displaystyle=$
$\displaystyle\max(G-K,G-A)+A-G$ $\displaystyle\leq$
$\displaystyle\max(G-K,0)+A-G,$

leads to the upper bound

$\displaystyle e^{-rT}E\left[\max(G-K,0)+A-G\right]$
$\displaystyle=e^{-rT}\left\\{E\left[\max(G-K)^{+}\right]+E[A]-E[G]\right\\}$
(3.11) $\displaystyle=C_{G}+e^{-rT}\left(E[A]-E[G]\right).$

Therefore, we have the following bounds on the arithmetic option price

(3.12) $\displaystyle C_{G}\leq C_{A}\leq
C_{G}+e^{-rT}\left(E[A]-E[G]\right).$

To compute the upper bound we need to compute the mean of arithmetic average,
which is obtained in Theorem 3.1.

######  Theorem 3.2.

The pricing error of arithmetic approximate vale $\tilde{C}$ is given by

(3.13) $\displaystyle|\tilde{C}-C_{A}|$ $\displaystyle=$ $\displaystyle
e^{-rT}\left(E[A]-E[G]\right).$

###### Proof.

Sine

(3.14) $\displaystyle\max(G-K,0)$ $\displaystyle\leq$
$\displaystyle\max\left(G-K+E[A]-E[G],0\right)$ $\displaystyle=$
$\displaystyle\max\left(G-K+E[G]-E[A]\right)+E[A]-E[G]$ $\displaystyle\leq$
$\displaystyle\max\left(G-K,0\right)+E[A]-E[G],$

it follows that

(3.15) $\displaystyle C_{G}$ $\displaystyle=$ $\displaystyle
e^{-rT}E\left[(G-K)^{+}\right]$ $\displaystyle\leq$ $\displaystyle
e^{-rT}E\left[(G-K+E[A]-E[G])^{+}\right]$ $\displaystyle=$
$\displaystyle\tilde{C}$ $\displaystyle\leq$ $\displaystyle
e^{-rT}E\left[(G-K)^{+}+E[A]-E[G]\right]$ $\displaystyle=$ $\displaystyle
C_{G}+e^{-rT}\left(E[A]-E[G]\right).$

The analogous bounds exert for the exact value $C_{A}$, then the absolute
difference between $\tilde{C}$ and $C_{A}$ is less than the difference between
the upper and lower bounds

(3.16) $\displaystyle|\tilde{C}-C_{A}|$ $\displaystyle\leq$ $\displaystyle
C_{G}+e^{-rT}\left(E[A]-E[G]\right)-C_{G}$ $\displaystyle=$ $\displaystyle
e^{-rT}\left(E[A]-E[G]\right).$

∎

##  4\. Asian power options

This section is deals with evaluate Asian power options by using actuarial
approach in a mixed fractional Brownian motion with jumps environment. From
[some references], we know that the payoff function of Asian options is given
by $(G^{m}(T)-K)^{+}$ and $(A^{m}(T)-K)^{+}$ for the geometric and arithmetic
call options with some integer $m\geq 1$, respectively.

######  Theorem 4.1.

Let the stock price $S(T)$ satisfy Equation (2.2). Then, the value of
geometric Asian power options with exercise price $K$ is given by

(4.1) $\displaystyle C(S_{T},K)$ $\displaystyle=$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\lambda
T}}{n!}\left[e^{-(r-q)T-\lambda(\rho-1)T+m\hat{\mu}+\frac{1}{2}m^{2}\hat{\sigma}^{2}}\Phi(d_{1})-Ke^{-rT}\Phi(d_{2})\right],$

for the call options, and

(4.2) $\displaystyle P(S_{T},K)$ $\displaystyle=$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\lambda
T}}{n!}\left[Ke^{-rT}\Phi(-d_{2})-\mathrm{e}^{-(r-q)T-\lambda(\rho-1)T+m\hat{\mu}+\frac{1}{2}m^{2}\hat{\sigma}^{2}}\Phi(-d_{1})\right],$

for the put options, where $d_{2}=\frac{\hat{\mu}-\ln
U^{\frac{1}{n}}}{\hat{\sigma}}$ and other parameters are defined in Theorem
2.1.

###### Proof.

For the geometric Asian power options, the payoff function is
$(G^{m}(T)-K)^{+}=(e^{mL(T)}-K)^{+}$. Using the strategy given in proof of
Theorem 2.1, then we have

(4.3) $\displaystyle C(S_{T},K)$ $\displaystyle=$ $\displaystyle
E\left[\left(\exp\left(-\int_{0}^{T}\theta(t)dt\right)G^{m}(T)-Ke^{-rT}\right)I_{B}\right]$
$\displaystyle=$ $\displaystyle
E\left[\left(e^{-(r-q)T-\lambda(\rho-1)T}G^{m}(T)-Ke^{-rT}\right)I_{B}\right]$
$\displaystyle=$ $\displaystyle
e^{-(r-q)T-\lambda(\rho-1)T}E\left[G^{m}(T)I_{B}\right]-Ke^{-rT}E\left[I_{B}\right]$
$\displaystyle=$ $\displaystyle
e^{-(r-q)T-\lambda(\rho-1)T}I_{1}-Ke^{-rT}I_{2},$

where $I_{1}=E\left[A(T)I_{B}\right]$ and $I_{2}=E\left[I_{B}\right]$. Let

(4.4) $\displaystyle B$ $\displaystyle=$
$\displaystyle\left\\{x:e^{mx}>U=Ke^{-qT-(1-\rho)\lambda T}\right\\}$
$\displaystyle=$ $\displaystyle\left\\{x:G^{m}(T)>U\right\\}$ $\displaystyle=$
$\displaystyle\left\\{y:e^{m(\hat{\mu}+\hat{\sigma}y)}>U\right\\}=\left\\{y:m(\hat{\mu}+\hat{\sigma}y)>\ln
U\right\\}$ $\displaystyle=$ $\displaystyle\left\\{y:y>-d_{2}=\frac{\ln
U^{\frac{1}{m}}-\hat{\mu}}{\hat{\sigma}}\right\\}.$

Observe that

(4.5) $\displaystyle I_{2}$ $\displaystyle=$ $\displaystyle
E\left[I_{B}\right]=\sum_{n=0}^{\infty}E\left[I_{G^{m}(T)>-d_{2}}\right]P\left[N_{T}=n\right]$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda T}}{n!}E\left[I_{G^{m}(T)>-d_{2}}\right]$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\lambda
T}}{n!}\int_{-d_{2}}^{\infty}\phi(y)dy$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\lambda
T}}{n!}\Phi(d_{2}),$

and

(4.6) $\displaystyle I_{1}$ $\displaystyle=$ $\displaystyle
E\left[G^{m}(T)I_{B}\right]=\sum_{n=0}^{\infty}E\left[G^{m}(T)I_{G^{m}(T)>-d_{2}}\right]P\left[N_{T}=n\right]$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda T}}{n!}E\left[e^{mL(T)}I_{G^{m}(T)>-d_{2}}\right]$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda
T}}{n!}\int_{B}e^{mx}\frac{1}{\sqrt{2\pi}\hat{\sigma}}e^{-\frac{(x-\hat{\mu})^{2}}{2\hat{\sigma}^{2}}}dx$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda
T}}{n!}\int_{B}e^{m(\hat{\mu}+\hat{\sigma}y)}\frac{1}{\sqrt{2\pi}\hat{\sigma}}e^{-\frac{(x-\hat{\mu})^{2}}{2\hat{\sigma}^{2}}}dy$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda
T}}{n!}e^{m\hat{\mu}+\frac{1}{2}m^{2}\hat{\sigma}^{2}}\int_{-d_{2}}^{\infty}\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}(y-m\hat{\sigma})^{2}}dy$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda
T}}{n!}e^{m\hat{\mu}+\frac{1}{2}m^{2}\hat{\sigma}^{2}}\int_{-d_{2}-\hat{\sigma}}^{\infty}\phi(y)dy$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda
T)^{n}e^{-\lambda
T}}{n!}e^{m\hat{\mu}+\frac{1}{2}m^{2}\hat{\sigma}^{2}}\Phi(d_{1}=d_{2}+\hat{\sigma}).$

From Eqs. (4.3), (4.5) and (4.6), we get desired result.

∎

In Figs 5-7 we show the price of the power Asian option computed using formula
4.1 as a function of the Hurst and the jump. parameters.

Figure 5. The variation of the power Asian option with respect to the Husrt
parameter $H$.

Figure 6. The variation of the power Asian option with respect to the
$\sigma_{J}$

.

Figure 7. The variation of the power Asian option with respect to the
$\mu_{J}$

.

Figure 8. The variation of the power Asian option with respect to the both
strike price $K$ and time to maturity (year) $T$.

In Fig 8 we show the price of the power Asian option computed using formula
4.1 as a function of the strike price and of the maturity.

The payoff function of the power Asian power call option is $(A^{m}(T)-K)$,
then its approximate value is given by the following theorem.

######  Theorem 4.2.

An approximate value of the arithmetic Asian power option is given by

(4.7) $\displaystyle\tilde{C}(S_{T},K)$ $\displaystyle=$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\lambda
T}}{n!}\left[e^{-(r-q)T-\lambda(\rho-1)T+\hat{\mu}+\frac{1}{2}\hat{\sigma}^{2}}\Phi(d_{1})-K^{\prime}e^{-rT}\Phi(d_{2})\right],$

for the call options, and

(4.8) $\displaystyle\tilde{P}(S_{T},K)$ $\displaystyle=$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(\lambda T)^{n}e^{-\lambda
T}}{n!}\left[K^{\prime}e^{-rT}\Phi(-d_{2})-\mathrm{e}^{-(r-q)T-\lambda(\rho-1)T+\hat{\mu}+\frac{1}{2}\hat{\sigma}^{2}}\Phi(-d_{1})\right],$

for the put options, where $d_{2}=\frac{\hat{\mu}-\ln
U^{\frac{1}{m}}}{\hat{\sigma}}$ and other parameters are defined in Theorem
3.1.

###### Proof.

Suppose $\bar{\mu}=E\left[A^{m}(T)\right]$, then

(4.9) $\displaystyle\bar{\mu}$ $\displaystyle=$ $\displaystyle
E\left[A^{m}(T)\right]=\frac{1}{T}\int_{0}^{T}E\left[(S_{t}^{n})^{m}\right]dt$
$\displaystyle=$
$\displaystyle\frac{S_{0}}{T}\int_{0}^{T}e^{m((r_{n}-\frac{1}{2}\sigma_{n}^{2})t+\frac{1}{2}m\sigma_{n}^{2}t)}=\frac{S_{0}^{m}}{T}\int_{0}^{T}e^{mr_{n}t}dt$
$\displaystyle=$
$\displaystyle\frac{S_{0}^{m}}{mr_{n}T}\left(e^{mr_{n}T}-1\right).$

From Theorem [24], to provide an approximation of the arithmetic Asian power
options, we need to calculate the new exercise price $K^{\prime}$,

(4.10) $\displaystyle K^{\prime}$ $\displaystyle=$ $\displaystyle
K-\left(E[A^{m}]-E[G^{m}]\right)$ $\displaystyle=$ $\displaystyle
K+E[G^{m}]-\bar{\mu},$

where $E[G^{m}]=e^{m\hat{\mu}+\frac{1}{2}m^{2}\hat{\sigma}^{2}}$ is calculated
in Theorem 4.1. Now, by replacing $K$ with $K^{\prime}$ in the proof of
Theorem 4.1, the proof is completed. ∎

##  5\. Conclusions

We derived approximate closed-form solutions for pricing arithmetic Asian
options on an underlying described by the mixed fractional Brownian motion
($mfBm$). In the last decade, the $mfBm$ has gained increasing popularity in
finance, since it allows one to incorporate long-range dependence and self-
similarity and, for certain values of the Hurst parameter, it is arbitrage-
free. In the present paper, suitable approximate closed-form solutions are
obtained for both arithmetic Asian options and arithmetic Asian power options,
by applying a convenient approximation of their strike prices. Both the
standard $mfBm$ and the $mfBm$ with Poisson log-normally distributed jumps are
considered.

## References

  * [1] D. Ahmadian and L. V. Ballestra.  Pricing geometric Asian rainbow options under the mixed fractional Brownian motion.  Physica A: Statistical Mechanics and its Applications, 555:124458, 2020.
  * [2] L. V. Ballestra, G. Pacelli, and D. Radi.  A very efficient approach for pricing barrier options on an underlying described by the mixed fractional Brownian motion.  Chaos, Solitons & Fractals, 87:240–248, 2016.
  * [3] C. Bender and R. J. Elliott.  Arbitrage in a discrete version of the Wick-fractional Black-Scholes market.  Mathematics of Operations Research, 29:935–945, 2004.
  * [4] Mogens Bladt and Tina Hviid Rydberg.  An actuarial approach to option pricing under the physical measure and without market assumptions.  Insurance: Mathematics and Economics, 22(1):65–73, 1998\.
  * [5] J. H Chen, F. Y. Ren, and W. Y. Qiu.  Option pricing of a mixed fractional-fractional version of the Black-Scholes model.  Chaos, Solitons & Fractals, 21:1163–1174, 2004.
  * [6] P. Cheridito.  Mixed fractional Brownian motion.  Bernoulli, 7:913–934, 2001.
  * [7] P. Cheridito.  Arbitrage in fractional Brownian motion models.  Finance and Stochastics, 7:533–553, 2003.
  * [8] Patrick Cheridito.  Arbitrage in fractional Brownian motion models.  Finance and Stochastics, 7(4):533–553, 2003.
  * [9] Z. Ding, C. W. J. Granger, and R. F. Engle.  A long memory property of stock market returns and a new model.  Journal of Empirical Finance, 1:83–106, 1993.
  * [10] X. He and W. Chen.  The pricing of credit default swaps under a generalized mixed fractional Brownian motion.  Physica A: Statistical Mechanics and its Applications, 404:26–33, 2014.
  * [11] Angelien GZ Kemna and Antonius CF Vorst.  A pricing method for options based on average asset values.  Journal of Banking & Finance, 14(1):113–129, 1990.
  * [12] A. W. Lo.  Long-term memory in stock market prices.  Econometrica, 59:1279–1313, 1991.
  * [13] Y. S. Mishura.  Stochastic Calcular for Fractional Brownian Motion and Related Process.  Springer–Verlag, Berlin, 2008.
  * [14] Y. S. Mishura and E. Valkeila.  The absence of arbitrage in a mixed Brownian fractional Brownian model.  In Proceedings of the Steklov Institute of Mathematics, number 237, pages 224–233. Trudy Matematicheskogo Instituta Imeni V.A. Steklova, 2002.
  * [15] B. Peng and F. Peng.  Pricing Asian power options under jump-fraction process.  Journal of Economics Finance and Administrative Science, 17:2–9, 2012.
  * [16] B. L. S. Prakasa Rao.  Pricing geometric Asian power options under mixed fractional Brownian motion environment.  Physica A: Statistical Mechanics and its Applications, 446:92–99, 2012.
  * [17] S. Rostek and R. Schöbel.  A note on the use of fractional Brownian motion for financial modeling.  Economic Modelling, 30:30–35, 2013.
  * [18] Foad Shokrollahi.  The evaluation of geometric Asian power options under time changed mixed fractional Brownian motion.  Journal of Computational and Applied Mathematics, 344:716–724, 2018.
  * [19] Foad Shokrollahi and Adem Kılıçman.  Pricing currency option in a mixed fractional Brownian motion with jumps environment.  Mathematical Problems in Engineering, 2014, 2014.
  * [20] Foad Shokrollahi and Adem Kılıçman.  Actuarial approach in a mixed fractional Brownian motion with jumps environment for pricing currency option.  Advances in Difference Equations, 2015(1):1–8, 2015.
  * [21] Steven E Shreve.  Stochastic calculus for finance II: Continuous-time models, volume 11.  Springer Science & Business Media, 2004.
  * [22] L. Sun.  Pricing currency options in the mixed fractional Brownian motion.  Physica A: Statistical Mechanics and its Applications, 392:3441–3458, 2013.
  * [23] Stuart M Turnbull and Lee Macdonald Wakeman.  A quick algorithm for pricing European average options.  Journal of Financial and Quantitative Analysis, pages 377–389, 1991.
  * [24] Ton Vorst.  Prices and hedge ratios of average exchange rate options.  International Review of Financial Analysis, 1(3):179–193, 1992.
  * [25] X.-T. Wang, E.-H. Zhu, M.-M. Tang, and H.-G. Yan.  Scaling and long-range dependence in option pricing II: Pricing European option with transaction costs under the mixed Brownian–fractional Brownian model.  Physica A: Statistical Mechanics and its Applications, 3:445–451, 2010.
  * [26] W. L. Xiao, W. G. Zhang, X. Zhang, and X. Zhang.  Pricing model for equity warrants in a mixed fractional Brownian environment and its algorithm.  Physica A: Statistical Mechanics and its Applications, 391:6418–6431, 2012.
  * [27] W. G. Zhang, Z. Li, and Y. J. Liu.  Analytical pricing of geometric Asian power options on an underlying driven by a mixed fractional Brownian motion.  Physica A: Statistical Mechanics and its Applications, 490:402–418, 2018.
  * [28] M. Zili.  On the mixed fractional Brownian motion.  Journal of Applied Mathematics and Stochastic Analysis, ID 32435:1–9, 2006.


# Temperature coefficients of crystalline-quartz elastic constants over the
cryogenic range [4 K, 15 K]

Jérémy Bon<EMAIL_ADDRESS>FEMTO-ST Institute, Univ. Bourgogne
franche-Comté, CNRS, ENSMM, 26 rue de l’Épitaphe, 25030 Besançon, France
Serge Galliou<EMAIL_ADDRESS>FEMTO-ST Institute, Univ. Bourgogne
franche-Comté, CNRS, ENSMM, 26 rue de l’Épitaphe, 25030 Besançon, France
Roger Bourquin<EMAIL_ADDRESS>FEMTO-ST Institute, Univ. Bourgogne
franche-Comté, CNRS, ENSMM, 26 rue de l’Épitaphe, 25030 Besançon, France

###### Abstract

This paper brings out results of a measurement campaign aiming to determine
the temperature coefficients of synthetic quartz elastic constants at liquid
helium temperature. The method is based on the relationship between the
resonance frequencies of a quartz acoustic cavity and the elastic constants of
the material. The temperature coefficients of the elastic constants are
extracted from experimental frequency-temperature data collected from a set of
resonators of various cut angles, because of the anisotropy of quartz,
measured on the very useful cryogenic range [4 K-15 K]. The knowledge of these
temperature coefficients would allow to further design either quartz
temperature sensors or conversely frequency-temperature compensated quartz
cuts. With extremely low losses, lower than $10^{-9}$ for the best ones, key
applications of such devices are ultra-low loss mechanical systems used in
many research areas including frequency control and fundamental measurements.
The Eulerian formalism is used in this study to identify the temperature
coefficients.

##  I Introduction

Bulk acoustic wave (BAW) piezoelectric devices are widely used in research and
industry to make temperature sensors, acoustic cavities or filters for several
applications. Recently it has been demonstrated that quartz acoustic cavities
can exhibit very high quality factors, up to a few billions for the best ones,
at liquid helium temperature Galliou  _et al._ (2013). In these conditions
they become very attractive for various experiments in fundamental physics Lo
_et al._ (2016); Goryachev and Tobar (2014), as well as for hybrid quantum
systems Goryachev  _et al._ (2013); Aspelmeyer (2010), optomechanics
Aspelmeyer  _et al._ (2014), etc. This interest is amplified by the fact that
commercial cryorefrigerators are now easy to use. In the field of ultra stable
frequency standards, the advent of superconducting devices also offers new
interesting opportunities. The coupling of those cavities exhibiting extremely
low mechanical losses with SQUID (Superconducting QUantum Interference Device)
is a first step towards cryogenic ultra-stable frequency sources Goryachev
_et al._ (2014). Nevertheless, progress in this field suffer of the lack of
knowledge regarding the physical constants of the acoustic cavity material,
the crystalline quartz, at cryogenic temperatures. The very first topic of
interest is the behavior of elastic constants versus temperature. Although a
set of temperature coefficients of the elastic constants of quartz already
exists since 1975 Smagin (1995), it covers the very large temperature range
$4{\rm K}-400{\rm K}$, and it turns out that it is not reliable within the
very useful range $4{\rm K}-15{\rm K}$ commonly achieved by standard
cryocoolers. In addition, these coefficients have been determined from the
fast and slow (quasi-)shear propagation modes usually denoted B and C
respectively, and not from the (quasi-)longitudinal mode, denoted A, the one
which exhibits the highest quality factors at low temperatures.

Several methods can be used to determine the thermal coefficients. It can be
from phase velocity measurements by pulse echo overlap methods, resonant
ultrasound spectroscopy or combined resonance techniques, or by means of laser
ultrasonic techniques. However, application of these methods at low
temperature arises technical problems. Alternatively, one way consists in
probing the temperature sensitivity of a resonance frequency from a dedicated
acoustic device. Indeed, the acoustic resonance frequency of a mechanical
system is a function of its elastic coefficients which are themselves
temperature dependent. This method is particularly relevant for a quartz
crystal device: the temperature sensitivity of the elastic coefficients of the
quartz material can be determined by identifying the frequency-to-temperature
coefficients of a plano-convex quartz resonator, i.e. a trapped energy
resonator, whose frequency of overtone order $n$ derives from the basic
relationship $f_{n}=\frac{n}{2h}\sqrt{\frac{C_{e}}{\rho}}$ where $2h$ is a
thickness, $C_{e}$ an elastic coefficient as described below, and $\rho$ the
quartz density. This latter approach is used here to extract the temperature
coefficients of quartz stiffness within $4{\rm K}-15{\rm K}$ from various
quartz-cut with different orientations relative to the crystallographic
reference.

##  II Method

The Eulerian formalism can be used to determine the thermal coefficients of
the elastic coefficients $C_{ijkl}$ because the linear thermal expansion
coefficients of quartz are known at the operating temperatures Barron  _et
al._ (1982).

Strictly speaking, according to the Eulerian description, the problem should
be solved in the so-called final state Dulmet and Bourquin (2001a), Dulmet and
Bourquin (2001b) by adding the mechanical strains caused by the acoustic
vibration to the temperature strains induced by the temperature change from
room temperature to cryogenic temperature. However, as an available
approximation, the problem can be solved in the intermediate state by taking
into account the strains induced by the very small vibration-caused
displacement magnitude, typically in the order of a nanometer.

In the Euler formalism the strain-stress relationship is given by the Cauchy
tensor expressed, for the intermediate state, by
$T_{ij}=C_{ijkl}\frac{\partial u_{k}}{\partial x_{l}}$ where $C_{ijkl}$ is the
elastic coefficient tensor and $u_{k}$ a component of the displacement vector.
The eigenfrequencies of the acoustical modes of the mechanical structure
obviously depend on the elastic coefficients $C_{ijkl}$ that are temperature
dependent. To bring up the thermal coefficients, the elastic coefficients can
be written in the form of a Taylor series:

$\begin{array}[]{l}C_{ijkl}(T)=C_{ijkl}(T_{ref})\left[1+\sum_{k}T^{(k)}C_{ijkl}\right.\vspace{8pt}\\\
\hskip 48.0pt\left.(T-T_{ref})^{k}\right]\end{array}$ (1)

where $C_{ijkl}(T_{ref})$ are the elastic coefficients at the reference
temperature $T_{ref}$, and $T^{(k)}C_{ijkl}$ the order $k$ thermal
coefficients defined by:

$T^{(k)}C_{ijkl}=\frac{1}{k!}\frac{1}{C_{ijkl}(T_{ref})}\frac{d^{k}C_{ijkl}}{dT^{k}}$
(2)

The eigenfrequencies of the plano-convex resonators measured in this study
(see Fig. 1) can be expressed according to the Stevens-Tiersten model Stevens
and Tiersten (1986):

$\begin{array}[]{l}f_{nmp}^{2}=f_{n}^{2}\left[1+\frac{1}{n\pi}\sqrt{\frac{2h}{R_{0}}}\left(\sqrt{\frac{M_{n}}{C_{e}}}(2m+1)\right.\right.\vspace{8pt}\\\
\hskip
32.0pt\left.\left.+\sqrt{\frac{P_{n}}{C_{e}}}(2p+1)\right)\right]\end{array}$
(3)

where index $n$ is the overtone number and $m,p$ the number of modal lines in
the resonator plane, $M_{n}$ and $P_{n}$ the dispersion constants, $2h$ the
resonator thickness, $R_{0}$ the curvature radius of the convex face, and
$f_{n}$ the eigenfrequencies of the corresponding infinite flat plate. As
shown in Fig. 2, the latter can usually be computed accurately from the
eigenfrequencies of the plano-convex resonator by combining both anharmonic
modes $n02$ and $n20$, and the corresponding quasi-harmonic one $n00$ as
follows:

$f_{n}^{2}=\frac{(6f_{n00}^{2}-f_{n02}^{2}-f_{n20}^{2})}{4}$ (4)

Figure 1: Plano-convex resonator. In such a trapped-energy device the
vibration displacement $u_{2}$ along the Y-axis has a gaussian distribution
Stevens and Tiersten (1986).

Figure 2: Fractional frequency deviation of $f_{5}$ $f_{502}$ and $f_{520}$
versus the temperature offset with respect to $T_{ref}=8\rm{K}$, in the case
of the single-rotated Y-cut called AT-cut, as an example. Plots of the
$5^{th}$ order resonance frequency $f_{5}$ of the infinite flat plate
calculated from the plano-convex resonator frequencies $f_{5mp}$ are in very
good agreement with (4).

As a consequence, the temperature coefficients $T^{(k)}C_{ijkl}$ of the
elastic constants can just be determined from the temperature coefficients of
the eigenfrequencies of the corresponding infinite flat plate. These
eigenfrequencies are given by Stevens and Tiersten (1986) :

${f_{n}^{\xi}}^{2}(T)=\frac{n^{2}}{4h^{2}(T)}\frac{\overline{C_{e}^{\xi}}(T)}{\rho(T)}\left(1-\frac{8k_{\alpha}^{2}(T)}{n^{2}\pi^{2}}-r\right)$
(5)

where $\overline{C_{e}^{\xi}}$ is the effective elastic coefficient,
$k_{\alpha}^{2}$ the coupling electromechanical coefficient, $h$ the resonator
thickness (about 1 mm; the electrode thickness is neglected because it is
small enough, a few hundreds of nanometers, to neglect the wave propagation
inside the electrodes), $\rho$ the mass density of quartz, $n$ the overtone
range and $r$ the ratio between masses of the resonator and electrodes
(typically $10^{-3}$). $\xi$ denotes A, B or C, referring to the extensional
propagation mode, the high-speed thickness mode, and the low-speed thickness
mode respectively.

Actually, terms inside the brackets of (5) can be reduced to unity in the
$T^{(k)}C_{ijkl}$ calculation. Indeed, the third term $r$ is a constant
independent of temperature, and then is not involved in this calculation.
Regarding the second term inside the brackets of (5), values of the
electromechanical coefficients $k_{\alpha}^{2}$ are known only at room
temperature, but are close to $10^{-4}$. In addition, as any physical
constant, their thermal sensitivity are smaller at low temperature than at
room temperature. Thus, this second term can also be reasonably neglected
especially for overtone orders greater than 3.

The effective elastic coefficients $\overline{C_{e}^{\xi}}$ are eigenvalues in
the Christoffel theory, given by Auld (1973) in the case of a piezoelectric
material:

$\overline{\Gamma_{il}}=C_{ijkl}^{E}n_{j}n_{k}+\frac{(e_{pij}n_{p})(e_{qkl}n_{q})}{(\epsilon_{jk}^{S}n_{j}n_{k})}$
(6)

where $C_{ijkl}^{E}$ are the elastic coefficients measured at constant
electric field, $e_{pij}$ and $e_{qkl}$ the piezoelectric constants,
$\epsilon_{jk}^{S}$ the dielectric constant at constant strain and $n_{\beta}$
(with $\beta=p,q,j$ or $k$) the directions of propagation. The second term of
(6) is linked to the electromechanical coefficient, and it turns out that it
can also be neglected in comparison with the first one.

To sum up, it results from the simplifications described above that
piezoelectricy is not involved anymore. According to these hypotheses the
issue is then limited to its pure mechanical nature, and it turns out that:
$\overline{C_{e}^{\xi}}\approx C_{e}^{\xi}$ and
$\overline{\Gamma_{il}}\approx\Gamma_{il}$.

Coefficients $C_{e}^{\xi}$ are eigenvalues of the Christoffel tensorRoyer and
Dieulesaint (2000):

$\Gamma_{il}{}^{\circ}u_{l}=\rho V^{2}{}^{\circ}u_{i}$ (7)

with $\Gamma_{il}=C_{ijkl}n_{j}n_{k}$, and where $V$ and ${}^{\circ}u_{l}$ are
respectively the phase velocity and polarization of the acoustic wave. Terms
of the symmetrical tensor $\Gamma_{il}$ are defined in the crystalline class
32 of quartz as:

$\begin{array}[]{lll}\Gamma_{il}=C_{i11l}n_{1}^{2}+C_{i22l}n_{2}^{2}+C_{i33l}n_{3}^{2}+(C_{i12l}+\vspace{4pt}\\\
\hskip
25.0ptC_{i21l})n_{1}n_{2}+(C_{i13l}+C_{i31l})n_{1}n_{3}+(C_{i23l}+\vspace{4pt}\\\
\hskip 25.0ptC_{i32l})n_{2}n_{3}\end{array}$ (8)

where the index $i$ and $l$ varies from 1 to 3.

The cubic shape of the frequency-temperature characteristic observed
experimentally implies that the thermal coefficients of the elastic constants
should be provided up to the $3^{rd}$ order. These thermal sensitivities can
be obtained by differentiating equation (5) with respect to temperature:

$\begin{array}[]{lll}T^{(1)}C_{e}\approx
2T^{(1)}f+T^{(1)}\rho+2T^{(1)}h\vspace{6pt}\\\
T^{(2)}C_{e}\approx\frac{1}{2}[T^{(1)}C_{e}]^{2}-[T^{(1)}f]^{2}-\frac{1}{2}[T^{(1)}\rho]^{2}-\vspace{6pt}\\\
\hskip 42.0pt[T^{(1)}h]^{2}+2T^{(2)}f+T^{(2)}\rho+2T^{(2)}h\vspace{6pt}\\\
T^{(3)}C_{e}\approx
2T^{(3)}f-2T^{(1)}fT^{(2)}f+\frac{2}{3}T^{(1)}f+\vspace{6pt}\\\ \hskip
42.0pt\frac{1}{3}[T^{(1)}\rho]^{3}-\frac{1}{3}[T^{(1)}C_{e}]^{3}+\frac{2}{3}[T^{(1)}h]^{3}-\vspace{6pt}\\\
\hskip 42.0ptT^{(1)}\rho
T^{(2)}\rho-2T^{(1)}hT^{(2)}h+2T^{(3)}h+\vspace{6pt}\\\ \hskip
42.0ptT^{(1)}C_{e}T^{(2)}C_{e}+T^{(3)}\rho\end{array}$ (9)

where $T^{(k)}\gamma=\frac{1}{\gamma!}\frac{d^{k}\gamma}{dT^{k}}$, with
$\gamma=C_{e},h,\rho$, or $f$.

The temperature coefficients of the thickness $2h$ are simply the expansion
coefficients of the quartz material within $4.5K-15K$, with $T_{ref}=8K$.
Their numerical values have been computed by a third order fitting of data
provided by Barron  _et al._ (1982) (see Tab. 1), and according to the proper
resonator cut angles.

Table 1: Interpolated values of expansion coefficients from Barron  _et al._ (1982) for the quartz material over the cryogenic range [4.5K-15K]. $\alpha_{\perp}$ and $\alpha_{//}$ are the linear expansion coefficients along the X-axis and Z-axis respectively. | $\alpha^{(1)},10^{-8}$ | $\alpha^{(2)},10^{-8}$ | $\alpha^{(3)},10^{-8}$
---|---|---|---
$\alpha_{\perp}$ $(\alpha_{11})$ | $0.9651$ | $0.5160$ | $0.1416$
$\alpha_{//}$ $(\alpha_{33})$ | $-0.5134$ | $-0.1202$ | $0.0241$

The temperature coefficients of the mass density $\rho$ also depend on the
previous expansion coefficients, and can be expressed as follows Dulmet and
Bourquin (1983):

$\begin{array}[]{lll}T^{(1)}\rho=-(2\alpha_{11}^{(1)}+\alpha_{33}^{(1)})\vspace{6pt}\\\
T^{(2)}\rho=-(2\alpha_{11}^{(2)}+\alpha_{33}^{(2)}-3(\alpha_{11}^{(1)})^{2}-\vspace{6pt}\\\
\hskip
44.0pt(\alpha_{33}^{(1)})^{2}-2\alpha_{11}^{(1)}\alpha_{33}^{(1)})\vspace{6pt}\\\
T^{(3)}\rho=-(2\alpha_{11}^{(3)}+\alpha_{33}^{(3)}-2\alpha_{11}^{(1)}(3\alpha_{11}^{(2)}+\vspace{6pt}\\\
\hskip
44.0pt\alpha_{33}^{(2)})-2\alpha_{33}^{(1)}(\alpha_{11}^{(2)}+\alpha_{33}^{(2)})+\vspace{6pt}\\\
\hskip
44.0pt3(\alpha_{11}^{(1)})^{2}\alpha_{33}^{(1)}+2\alpha_{11}^{(1)}(\alpha_{33}^{(1)})^{2}+\vspace{6pt}\\\
\hskip 44.0pt4(\alpha_{11}^{(1)})^{3}+(\alpha_{33}^{(1)})^{3})\end{array}$
(10)

Regarding the frequency differentiation, it is obtained from a $3^{rd}$ order
fitting of a set of experimental frequency-temperature plots recorded from
quartz resonators of various cut angles described in Tab. 2. In the Eulerian
description, six frequency-temperature graphs are necessary to determine the
three temperature coefficients of each of the seven elastic coefficients of
the tensor $C_{ijkl}$ written in matrix notation as:

$\begin{pmatrix}\@add@centering\centering C_{11}&C_{12}&C_{13}&C_{14}&0&0\\\
C_{12}&C_{11}&C_{13}&-C_{14}&0&0\\\ C_{13}&C_{13}&C_{33}&0&0&0\\\
C_{14}&-C_{14}&0&C_{44}&0&0\\\ 0&0&0&0&C_{44}&C_{14}\\\
0&0&0&0&C_{14}&C_{66}\\\ \end{pmatrix}$ (11)

In (11), the elastic coefficients are given according to the specific rules of
index contraction:

$(11)\leftrightarrow 1$ |  $(22)\leftrightarrow 2$ | $(33)\leftrightarrow 3$
---|---|---
$(23)=(32)\leftrightarrow 4$ |  $(13)=(31)\leftrightarrow 5$ | $(12)=(21)\leftrightarrow 6$

Cuts and resonance modes used for the calculation of the temperature
coefficients are summarized in Tab. 3.

The temperature coefficients $C_{e}^{\xi}$ can be linked to the temperature
behavior of the elastic coefficients by:

$T^{(k)}C_{e}^{\xi}=\left[A_{\xi j}\right]^{-1}T^{(k)}C_{j}$ (12)

Where $\left[A_{\xi j}\right]=\frac{C_{e}^{\xi}}{C_{j}}\frac{\partial
C_{e}^{\xi}}{\partial C_{j}}$ and $C_{j}$ are the elastic coefficients
$C_{ij}$ arranged in a column vector where index $j$ defines the component
location in this vector. Values of the elastic coefficients $C_{ij}$ of
$\left[A_{ij}\right]$ are those provided by Bechmann (1958). Although these
values are given at room temperature, they can also be used at cryogenic
temperature with a relative accuracy better than $10^{-3}$, the fractional
frequency shifts observed from room temperature to cryogenic temperature
operations being typically lower than $10^{-3}$.

Table 2: Elastic coefficients and quartz-cuts used for their respective determination. the X, Y, and Z-cuts are cuts whose planes are normal to the X, Y, and Z cristalline axis respectively. The AT-cut is a single rotated Y-cut, i.e. a cut whose cut plane, originally normal to the Y-axis, is rotated around the X-axis. The SC-cut is a doubly-rotated cut, i.e. a AT-cut rotated around $Z^{\prime}$ the new Z-axis after the first rotation. The elastic coefficient $C_{66}$ is related to $C_{11}$ and $C_{12}$ as $C_{66}=(C_{11}-C_{12})/2$ Royer and Dieulesaint (2000). Elastic coefficient | Quartz cut | Comments
---|---|---
$C_{11}$ | X |
$C_{12}$ | X | Calculated from $C_{11}$ and $C_{66}$
$C_{13}$ | SC |
$C_{14}$ | AT |
$C_{33}$ | SC |
$C_{44}$ | Z |
$C_{66}$ | Y |
Table 3: Cuts of crystal plates and modes used to measure the corresponding resonance-frequency coefficients from a $3^{rd}$ order polynomial fit, with $T_{ref}=8$K. A, B and C modes are the (quasi-)longitudinal, fast (quasi-)shear and slow (quasi-)shear thickness modes respectively. Quartz cut |  Mode | $T^{(1)}f,10^{-7}$ | $T^{(2)}f,10^{-8}$ | $T^{(3)}f,10^{-9}$
---|---|---|---|---
AT |  $C_{300}$, $C_{500}$ | 0.769 | 0.748 | 0.105
SC |  $A_{300}$, $A_{500}$ | -1.178 | -2.189 | -2.474
SC |  $B_{300}$, $B_{302}$, $B_{320}$, $B_{500}$ | -0.879 | -1.672 | -1.832
X |  $A_{100}$, $A_{1300}$, $A_{1500}$ | -0.143 | -0.451 | -0.581
Y |  $C_{300}$, $C_{500}$ | 1.339 | 1.869 | 1.289
Z |  $C_{300}$, $C_{500}$ | -1.465 | -2.662 | -3.571

##  III Experiment

The temperature coefficients of the elastic tensor of the quartz material at
liquid helium temperature have been calculated from the set of experimental
frequency-temperature characteristics of various cut angles. Actually the
resulting frequency coefficients in Tab. 3, used to determine these
temperature coefficients, are averaged values from measurements made with
several overtone numbers of each tested cut.

The general overview of the measurement procedure can be summed up by the
following steps: $1/$ Resonators are realized in a set of plates of various
cut angles according to the conventional process: plate orientation with
respect to the quartz crystal axis, plate cutting, rounding, lapping and
polishing to obtain a 13 mm diameter plano-convex device, gold-electrode
coating for piezo-electrically exciting/probing the acoustic cavity. The
thickness of these resonators is such that the $3^{rd}$ overtone of the C-mode
is close to 10 MHz at room temperature. Cuts used to take into account the
quartz anisotropy comply with the IEEE-std 176-1978 IEE (1978): AT
($0^{\circ},-35.3^{\circ}$), SC ($22^{\circ},-34^{\circ}$), X
($30^{\circ},0^{\circ}$), Y ($0^{\circ},0^{\circ}$), Z
($0^{\circ},90^{\circ}$). $2/$ To be tested, the resonators are embedded
inside an oxygen-free cooper block screwed on the cold stage of a pulse-tube
cryorefrigerator (a Sumitomo cryocooler SRP082B). This block is cooled down in
the vicinity of 4 K and temperature controlled within $\pm 3$ mK around the
set temperature value when the steady state is reached. $3/$ For a given
temperature, the resonator admittance $Y$ is measured around the interesting
resonance frequencies, and over a frequency span properly chosen, by means of
a network analyzer. However a calibration procedure should be followed before
recording, to cancel the influence of the coaxial cable feeding the resonator
inside the vacuum chamber from the network analyzer: this calibration involves
three identical coaxial cables respectively ended by a standard 50 ohms
resistor, an open circuit and a short circuit. The network analyzer (HP4195A)
is locked to a hydrogen maser to get the frequency-stability performances of
the latter. The analyzer provides a $0.001$ dB and $0.01^{\circ}$ resolution
in magnitude and phase respectively, and its frequency resolution reaches $1$
mHz. In addition, amplitude-to-frequency nonlinear effects are avoided by
limiting the excitation power. $4/$ The post-processing consists in plotting
the admittance circle $Y=G+jB$ from the experimental data, magnitude and phase
of $Y$, to get the resonance frequency with the best accuracy. Fig. 3 shows an
example of such a G-B plot.

Figure 3: a) G-B plot from the network analyzer for the $1^{st}$ OT of a X-cut
at 4.5K. $\Delta f$ is with respect to $f_{ref}=4.580957\rm{Hz}$. b)
Admittance circle from the experimental data above.

##  IV Results

As an example, Fig. 4 shows experimental results of the fractional frequency
deviation versus temperature for a X-cut operating around $T_{ref}=8K$. The
geometrical properties of this device are $h=600\mu m$ and $R_{0}=250$ mm. It
can be observed that all the A-mode overtones exhibit similar behaviors.The
temperature coefficients of the frequency are extracted from such experimental
data (see Tab. 3). As a final result, Tab. 4 gives the thermal coefficients
$T^{(k)}C_{ij}$ of the elastic coefficients $C_{ij}$ of quartz, up to the
$3^{rd}$ order, determined according to the Eulerian description.

Figure 4: Resonance fractional frequency shift
$\frac{f(T)-f(T_{ref})}{f(T_{ref})}$ for a X-cut as a function of temperature
offset $\delta T=T-T_{ref}$, where $T_{ref}=8{\rm K}$.

Table 4: $1^{st}$, $2^{nd}$, and $3^{rd}$ order temperature coefficients of the elastic constants of crystalline quartz within 4 K - 15 K, with $T_{ref}=8{\rm K}$. $C_{ij}$ | $T^{(1)}C_{ijkl},10^{-7}$ | $T^{(2)}C_{ijkl}10^{-8}$ | $T^{(3)}C_{ijkl},10^{-9}$
---|---|---|---
$C_{11}$ | $-0.235$ | $1.28$ | $-3.25$
$C_{12}$111Calculated coefficient | $-6.39$ | $-10.56$ | $-4.22$
$C_{13}$ | $-7.18$ | $-7.42$ | $-0.877$
$C_{14}$ | $-1.35$ | $-1.31$ | $-7.93$
$C_{33}$ | $-4.62$ | $-5.65$ | $-7.62$
$C_{44}$ | $-3.17$ | $-4.37$ | $-11.6$
$C_{66}$ | $3.08$ | $5.92$ | $0.486$

##  V Conclusion

Various cut-angles of synthetic quartz have been measured at cryogenic
temperature in order to determine the temperature coefficients of its elastic
coefficients. The Eulerian formalism have been used to provide these
$T^{(k)}C_{ij}$ up to the $3^{rd}$ order. Uncertainties on the elastic
coefficient values are not specified because most of the uncertainties linked
to the many steps of measurements and calculation are difficult to bring out.
Actually, in addition to the cut angle uncertainty, one of the higher
uncertainties is the result of the Tiersten model which is a theoretical
approach modelling plano-convex resonators. However, it can be stated that the
uncertainty on the temperature sensitivity of the elastic coefficient $C_{13}$
is higher than that of other coefficients because of its very slight
involvement in calculations.

## Acknowledgements

Thanks to Philippe Abbé and Xavier Vacheret for their help for making
resonators and experiments. This work is supported by Conseil Régional de
Franche-Comté. J.B. is thankful to Ministère de l’Enseignement Supérieur et de
la Recherche, France, for his grant.

## References

  * Galliou  _et al._ (2013) S. Galliou, M. Goryachev, R. Bourquin, P. Abbé, J. P. Aubry, and M. E. Tobar, Sci. Rep. 3, 2132 (2013).
  * Lo  _et al._ (2016) A. Lo, P. Haslinger, E. Mizrachi, L. Anderegg, H. Müller, M. Hohensee, M. Goryachev, and M. Tobar, Phys. Rev. X 6, 011018 (2016).
  * Goryachev and Tobar (2014) M. Goryachev and M. E. Tobar, Phys. Rev. D 90, 102005 (2014).
  * Goryachev  _et al._ (2013) M. Goryachev, D. L. Creedon, S. Galliou, and M. E. Tobar, Phys. Rev. Lett. 111, 085502 (2013).
  * Aspelmeyer (2010) M. Aspelmeyer, Nature 464, 685 (2010).
  * Aspelmeyer  _et al._ (2014) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Rev. Mod. Phys. 86, 1391 (2014).
  * Goryachev  _et al._ (2014) M. Goryachev, E. Ivanov, F. V. Kann, S. Galliou, and M. E. Tobar, Appl. Phys. Lett. 105, 153505 (2014).
  * Smagin (1995) A. Smagin, in  _1995 IEEE International Frequency Control Symposium_ (1995).
  * Barron  _et al._ (1982) T. Barron, J. Collins, T. Smith, and G. White, J. of Phys. C 15 (1982).
  * Dulmet and Bourquin (2001a) B. Dulmet and R. Bourquin, J. Acoust. Soc. Am. 110 (2001a).
  * Dulmet and Bourquin (2001b) B. Dulmet and R. Bourquin, J. Acoust. Soc. Am. 110 (2001b).
  * Stevens and Tiersten (1986) D. Stevens and H. Tiersten, J. Acoust. Soc. Am. 79 (1986).
  * Auld (1973) B. Auld,  _Acoustic Fields and Waves in Solids. Volume I Volume Edition_ (Wiley, 1973).
  * Royer and Dieulesaint (2000) D. Royer and E. Dieulesaint,  _Elastic Waves in Solids I, Free and Guided Propagation_ (Springer, 2000).
  * Dulmet and Bourquin (1983) B. Dulmet and R. Bourquin, Revue de physique appliquée 18 (1983).
  * Bechmann (1958) R. Bechmann, Phys. Rev. 110 (1958).
  * IEE (1978) _IEEE Std 176-1978_ (IEEE Inc., September 29, 1978).