# The Gaussian free-field as a stream function: continuum version of the
scale-by-scale homogenization result

Peter Morfe<EMAIL_ADDRESS>,  Felix Otto<EMAIL_ADDRESS>and
Christian Wagner<EMAIL_ADDRESS>
###### Abstract.

This note is about a drift-diffusion process $X$ with a time-independent,
divergence-free drift $b$, where $b$ is a smooth Gaussian field that
decorrelates over large scales. In two space dimensions, this just fails to
fall into the standard theory of stochastic homogenization, and leads to a
borderline super-diffusive behavior. In a previous paper by Chatzigeorgiou,
Morfe, Otto, and Wang (2022), precise asymptotics of the annealed second
moments of $X$ were derived by characterizing the asymptotics of the effective
diffusivity $\lambda_{L}$ in terms of an artificially introduced large-scale
cut-off $L$. The latter was carried out by a scale-by-scale homogenization,
and implemented by monitoring the corrector $\phi_{L}$ for geometrically
increasing cut-off scales $L^{+}=ML$. In fact, proxies
$(\tilde{\phi}_{L},\tilde{\sigma}_{L})$ for the corrector and flux corrector
were introduced incrementally and the residuum $f_{L}$ estimated.

In this short supplementary note, we reproduce the arguments of the above
paper in the continuum setting of $M\downarrow 1$. This has the advantage that
the definition of the proxies $(\tilde{\phi}_{L},\tilde{\sigma}_{L})$ becomes
more transparent – it is given by a simple Itô SDE with $\ln L$ acting as a
time variable. It also has the advantage that the residuum $f_{L}$, which is a
martingale, can be efficiently and precisely estimated by Itô calculus. This
relies on the characterization of the quadratic variation of the (infinite-
dimensional) Gaussian driver.

##  1\. Definition of the proxies via Itô SDEs

For the set-up, bibliography, and the notation, we refer the reader to [2].

In this note we adopt a differential-geometric language and basis-free
notation: Since the generic coordinate direction $\xi$ acts as a linear form
on $\mathbb{R}^{2}$, i. e. a cotangent vector, we combine the two (standard)
scalar corrector fields $\\{\phi^{i}\\}_{i=1,2}$ to a tangent vector field
$\phi$, which amounts to111in this note, we use Einstein’s notation in the
strict form: We add over repeated indices provided one is a super and the
other a sub-script $\phi=\phi^{i}e_{i}$, where the vectors
$\\{e_{i}\\}_{i=1,2}$ form the Cartesian basis. Hence the contraction
$\xi.\phi$ yields the corrector in direction $\xi$. Thinking of $\nabla$ as
the differential, which applied to a scalar field yields a one-form, that is,
a cotangent vector field, we interpret $\nabla\phi$ as a field of
endomorphisms of cotangent space. Interpreting also $a$ as such a field, the
product $a\nabla\phi$ makes sense.

Recall that the stream functions $\psi_{L}$’s are coupled such that the
process $[1,\infty)\ni L\mapsto\psi_{L}$ has independent increments. In view
of $\mathbb{E}\psi_{L}^{2}=\ln L$, see [2, (13)], $[0,\infty)\ni\ln
L=:s\mapsto\psi_{L}$ acts like Brownian motion, with values in the space of
(smooth) stationary Gaussian scalar fields. In terms of this continuum
variable $s$, the discrete relation [2, (49)] for the proxy $\tilde{\lambda}$
to the effective diffusivity turns into

(1) $\displaystyle
d\tilde{\lambda}=\frac{\varepsilon^{2}}{2\tilde{\lambda}}ds\quad\mbox{and}\quad\tilde{\lambda}|_{s=0}=1$

in the limit $M\searrow 1$ (see [2, (49)]). In view of the Brownian character
of $s\mapsto\psi_{L}$, its distributional derivative/infinitesimal increment
$s\mapsto d\psi$ acts like (stationary, function valued, and Gaussian) white
noise (in $s$). In line with [2, (38)], we define the (stationary and vector-
field valued) white noises $d\phi$ and $d\sigma$ via the Helmholtz
decomposition

(2) $\displaystyle\tilde{\lambda}\nabla d\phi+d\psi J=J\nabla
d\sigma\quad\mbox{and}\quad\mathbb{E}d\phi=\mathbb{E}d\sigma=0.$

As discussed above, the white noises $\nabla d\phi$ and $\nabla d\sigma$ have
values in the space of fields of endomorphisms of cotangent space, so that
provided we interpret $J$ as the counter-clockwise rotation by $\frac{\pi}{2}$
on (co)tangent space, (2) makes sense.

Equipped with the vectorial driver222in the jargon of SDEs $d\phi$, we now
define the vectorial proxy corrector $\tilde{\phi}$ via an SDE in Itô form
that encapsulates the two-scale expansion:

(3) $\displaystyle
d\tilde{\phi}=\varepsilon(1+\tilde{\phi}^{i}\partial_{i})d\phi,\quad\tilde{\phi}|_{s=0}=0.$

On the basis of the vectorial driver $d\sigma$ and the original scalar driver
$d\psi$ we define the proxy flux corrector $\tilde{\sigma}$ via

(4) $\displaystyle
d\tilde{\sigma}=\varepsilon(d\sigma+\tilde{\sigma}^{i}\partial_{i}d\phi+d\psi\,\tilde{\phi})+J\tilde{\phi}\,\frac{\varepsilon^{2}}{2\tilde{\lambda}}ds,\quad\tilde{\sigma}|_{s=0}=0.$

Definitions (3) & (4) correspond to [2, (36) & (37)], the motivation for (4)
lies in (7) below. Note that in view of the drift term in (4),
$\tilde{\sigma}$ is not a martingale, as opposed to $\tilde{\phi}$.

As in [2], we measure the quality of the proxies $\tilde{\phi}$ and
$\tilde{\sigma}$ in terms of the residuum $f$ in

(5) $\displaystyle a({\rm id}+\nabla\tilde{\phi})=\tilde{\lambda}\,{\rm
id}+J\nabla\tilde{\sigma}+f,$

where we recall the definition of $a$ from [2, (11)]

(6) $\displaystyle a={\rm id}+\varepsilon\psi
J\quad\mbox{with}\quad\psi|_{s=0}=0.$

In this note, we give a self-contained proof for333what amounts to a first-
order error estimate in $\varepsilon$

######  Proposition 1.1 (see [2, (79)]).

For $\varepsilon^{2}\ll 1$ it holds

$\mathbb{E}|f|^{2}\lesssim\varepsilon^{2}\tilde{\lambda}.$

In preparation, we will establish in the next section that (5) is consistent
with the tensor-valued444Note that as a tensor product of a cotangent vector
and a tangent vector, the last term in (7) canonically is an endomorphism of
cotangent vector space. SDE

(7) $\displaystyle df=\varepsilon\big{(}f\nabla
d\phi+(\tilde{\phi}^{i}a-\tilde{\sigma}^{i}J)\nabla\partial_{i}d\phi-(J\nabla
d\psi)\otimes\tilde{\phi}\big{)},\quad f|_{s=0}=0,$

see [2, (50)]. Note that (7) shows that the martingale $f$ has vanishing
expectation, which was part of its definition in [2, (48)]. This implies the
exact consistency $\mathbb{E}a({\rm id}+\nabla\tilde{\phi})$
$=\tilde{\lambda}{\rm id}$, which in [2] only holds in the limit $M\searrow
1$, see [2, (51)].

The approach of a continuum decomposition of an underlying Gaussian field,
often the Gaussian free field $\psi$ like here, according to spatial scales
$L$ is also used in quantum field theory (where the underlying noise is
thermal instead of quenched/environmental), and is known as variance
decomposition. For instance, [1] use it to control (exponential) moments of
the field $\phi$ under consideration, with help of suitable martingales in
$\ln L$, analogous to $\tilde{\phi}$. The variance decomposition also induces
a Hamilton-Jacobi-type evolution equation for the effective Hamiltonian in
$\ln L$, known as Polchinski equation. In [3, Section 2.3], to cite another
recent work, proxies for the solution of the Polchinski equation are
constructed by truncation of an expansion (analogous to an expansion in
$\varepsilon$ here), and the evolution of the residuum is monitored, like $f$
in (7).

##  2\. Proof of (5) by identification of quadratic variations

We show how the SDEs (3), (4), and (7) imply the identity (5). To this end, we
monitor the tensor-field valued residuum

(8) $\displaystyle r:=a({\rm id}+\nabla\tilde{\phi})-\tilde{\lambda}\,{\rm
id}-J\nabla\tilde{\sigma}-f$

and shall establish the homogeneous SDE

(9) $\displaystyle dr=\varepsilon r\nabla d\phi,\quad r|_{s=0}=0.$

The second item in (9) follows directly from the second items in (1), (3),
(4), (7), and (6).

In order to establish the first item in (9), and in view of (8), we need to
identify the infinitesimal increment of the flux $a({\rm
id}+\nabla\tilde{\phi})$, which is a bilinear term in the drivers and thus
involves a quadratic variation term. We claim that

(10) $\displaystyle d\big{(}a({\rm
id}+\nabla\tilde{\phi})\big{)}=\varepsilon(-\tilde{\lambda}\nabla
d\phi+J\nabla d\sigma)({\rm id}+\nabla\tilde{\phi})+a\nabla
d\tilde{\phi}+\varepsilon^{2}({\rm
id}-J\nabla\tilde{\phi}J)\frac{ds}{2\tilde{\lambda}}.$

Indeed, applying $\nabla$ to (3) and appealing to Leibniz’ rule in conjunction
with the identity555where the latter is understood as a composition of
endomorphisms $(\nabla\tilde{\phi}^{i})\otimes\partial_{i}d\phi$
$=\nabla\tilde{\phi}\nabla d\phi$, we obtain

(11) $\displaystyle d\nabla\tilde{\phi}=\varepsilon({\rm
id}+\nabla\tilde{\phi})\nabla
d\phi+\varepsilon\tilde{\phi}^{i}\nabla\partial_{i}d\phi.$

In view of (6), the main task is to identify the infinitesimal increment of
$\varepsilon\psi J\nabla\tilde{\phi}$; by Itô calculus this involves two
covariations of the drivers:

(12) $\displaystyle d\big{(}\varepsilon\psi
J\nabla\tilde{\phi}\big{)}=\varepsilon d\psi
J\nabla\tilde{\phi}+\varepsilon\psi J\nabla
d\tilde{\phi}+\varepsilon^{2}J({\rm id}+\nabla\tilde{\phi})[d\psi\,\nabla
d\phi]+\varepsilon^{2}\tilde{\phi}^{i}[d\psi\,\nabla\partial_{i}d\phi].$

These covariations of drivers are deterministic by the independence of
increments and constant in space by stationarity. Since $d\psi$ is invariant
in law under point reflection, $d\phi$ is (jointly) odd under the same
operation in view of (2), and thus also $\nabla\partial_{i}d\phi$. Hence the
last contribution to (12) vanishes, i. e.

(13) $\displaystyle[d\psi\,\nabla\partial_{i}d\phi]=0.$

Thus the task is to identify the (deterministic) endomorphism $B$ (of
cotangent space) in

(14) $[d\psi\,\nabla d\phi]=Bds.$

Since the scalar $d\psi$ is invariant in law under rotations, the endomorphism
$\nabla d\phi$ is invariant in law under conjugation with rotations. This
transmits to $B$, which amounts to

(15) $\displaystyle B=\;\mbox{scalar}\;\times\;\mbox{rotation}.$

In order to identify $B$, we derive local relations between the drivers
$d\phi$ and $d\psi$ from (2):

(16) $\displaystyle{\rm tr}\nabla
d\phi=0\quad\mbox{and}\quad\tilde{\lambda}{\rm tr}\nabla d\phi J=d\psi.$

To this purpose, we consider the $i$-th component of (2), which amounts to
applying this endomorphism identity to the cotangent vector $e^{i}$ (where
$\\{e^{i}\\}_{i=1,2}$ denotes the dual basis) and yields
$\tilde{\lambda}\nabla d\phi^{i}+d\psi Je^{i}$ $=J\nabla d\sigma^{i}$.
Applying the divergence666where for a cotangent vector $\xi$, $\xi\cdot$
denotes the corresponding tangent vector $\nabla\cdot=\partial_{j}e^{j}\cdot$
and using $\nabla\cdot\nabla=\triangle$ gives $\tilde{\lambda}\triangle
d\phi^{i}+\partial_{j}d\psi(e^{j}\cdot Je^{i})=0$. In view of
$d\phi^{i}e_{i}=d\phi$ and the skewness of $J$ in form of $(e^{j}\cdot
Je^{i})e_{i}=-(Je^{j}\cdot e^{i})e_{i}=-(Je^{j})\cdot$, we obtain the tangent-
vector identity

(17) $\displaystyle\tilde{\lambda}\triangle d\phi=(J\nabla d\psi)\cdot.$

Applying $\nabla=\partial_{i}e^{i}$, and noting that
$e^{i}\otimes\partial_{i}\nabla d\psi\cdot$ coincides with the Hessian
endomorphism $\nabla^{2}d\psi$ of cotangent space, we obtain once more by
skewness of $J$ the tensor-valued identity

(18) $\displaystyle\tilde{\lambda}\triangle\nabla d\phi=-\nabla^{2}d\psi J.$

Taking the trace and appealing to the symmetry of $\nabla^{2}d\psi$ on the one
hand, and multiplying with $J$ from the right, using $J^{2}=-{\rm id}$, and
then taking the trace and appealing to ${\rm tr}\nabla^{2}=\triangle$ on the
other hand gives

$\displaystyle\triangle{\rm tr}\nabla
d\phi=0\quad\mbox{and}\quad\triangle(\tilde{\lambda}{\rm tr}\nabla d\phi
J-d\psi)=0.$

Since both expressions under $\triangle$ are stationary fields of vanishing
expectation, this yields (16).

Because of the underlying white-noise type normalization

(19) $\displaystyle[d\psi\,d\psi]=ds,$

on the level of (14), and the locality777in differential geometric jargon of
the covariation seen as a bilinear form, (16) translates into

$\displaystyle{\rm tr}B=0\quad\mbox{and}\quad\tilde{\lambda}{\rm tr}BJ=1,$

which in view of (15) finally implies $B=-\frac{1}{2\tilde{\lambda}}J$. We
thus have identified (14):

(20) $\displaystyle[d\psi\,\nabla d\phi]=-J\frac{ds}{2\tilde{\lambda}}.$

Inserting (13) and (20) into (12) yields

$\displaystyle d\big{(}\varepsilon\psi J\nabla\tilde{\phi}\big{)}=\varepsilon
d\psi J\nabla\tilde{\phi}+\varepsilon\psi J\nabla
d\tilde{\phi}+\varepsilon^{2}({\rm
id}-J\nabla\tilde{\phi}J)\frac{ds}{2\tilde{\lambda}},$

which by (2) and (6) turns into (10).

Equipped with (10) we finally turn to (9). We now argue that the first part of
the drift term in (10) is balanced by (1) whereas the second part is
compensated by (4). Indeed, applying $\nabla$ to (4) we obtain because of
$\nabla\tilde{\sigma}^{i}\otimes\partial_{i}d\phi=\nabla\tilde{\sigma}\nabla
d\phi$ and $\nabla(J\tilde{\phi})=\nabla\tilde{\phi}J^{*}$, where $J^{*}$
denotes the adjoint endomorphism w. r. t. the Euclidean inner product on the
(co)tangent space,

(21) $\displaystyle dJ\nabla\tilde{\sigma}=\varepsilon(J\nabla
d\sigma+J\nabla\tilde{\sigma}\nabla
d\phi+\tilde{\sigma}^{i}J\nabla\partial_{i}d\phi+J\nabla
d\psi\otimes\tilde{\phi}+d\psi
J\nabla\tilde{\phi})-\varepsilon^{2}J\nabla\tilde{\phi}J\,\frac{ds}{2\tilde{\lambda}}.$

Applying $d$ to (8) and inserting (1), (11), (10), and (21) yields

$\displaystyle dr$ $\displaystyle=-\varepsilon\tilde{\lambda}(\nabla
d\phi)({\rm id}+\nabla\tilde{\phi})+\varepsilon J(\nabla
d\sigma)\nabla\tilde{\phi}+\varepsilon a({\rm id}+\nabla\tilde{\phi})\nabla
d\phi$
$\displaystyle+\varepsilon\tilde{\phi}^{i}a\nabla\partial_{i}d\phi-\varepsilon
J\nabla\tilde{\sigma}\nabla
d\phi-\varepsilon\tilde{\sigma}^{i}J\nabla\partial_{i}d\phi-\varepsilon
J\nabla d\psi\otimes\tilde{\phi}-\varepsilon d\psi J\nabla\tilde{\phi}-df,$

which by definition (2) of $(\nabla d\phi,\nabla d\sigma)$ in terms of $d\psi$
simplifies to

$\displaystyle dr$ $\displaystyle=-\varepsilon\tilde{\lambda}\nabla
d\phi+\varepsilon a({\rm id}+\nabla\tilde{\phi})\nabla d\phi$
$\displaystyle+\varepsilon\tilde{\phi}^{i}a\nabla\partial_{i}d\phi-\varepsilon
J\nabla\tilde{\sigma}\nabla d\phi-\varepsilon
J\tilde{\sigma}^{i}\nabla\partial_{i}d\phi-\varepsilon J\nabla
d\psi\otimes\tilde{\phi}-df.$

This motivates the form of (7) that compensates the terms that have an
additional derivative on the drivers $\nabla d\phi$ and $d\psi$:

$\displaystyle dr=\varepsilon(-\tilde{\lambda}{\rm id}+a({\rm
id}+\nabla\tilde{\phi})-J\nabla\tilde{\sigma}-f)\nabla d\phi,$

which by definition (8) of $r$ turns into the desired identity (9).

##  3\. Proof of estimates by Itô calculus

Lemma 2 from [2] provides an iteration formula for $f_{L+}$ that is used to
estimate $\mathbb{E}|f_{L}|^{2}$. In the continuum approach, based on the SDE
(7), we use Itô calculus to derive the differential inequality

(22)
$d\mathbb{E}|f|^{2}-\varepsilon^{2}\mathbb{E}|f|^{2}\frac{ds}{2\tilde{\lambda}^{2}}\leq\varepsilon^{2}\mathbb{E}|\tilde{\phi}\otimes
a-\tilde{\sigma}\otimes
J|^{2}\frac{ds}{2L^{2}\tilde{\lambda}^{2}}+\varepsilon^{2}\mathbb{E}|\tilde{\phi}|^{2}\frac{2ds}{L^{2}}.$

Here, $|f|^{2}$ stands for the squared Frobenius norm given by
$|f|^{2}:=\operatorname{\text{tr}}f^{*}f$. Likewise, $|\tilde{\phi}\otimes
a-\tilde{\sigma}\otimes J|^{2}$ denotes the squared Frobenius norm of that
3-tensor.

###### Proof of [2, Lemma 2] revisited..

Since $f$ is a martingale, its quadratic variation determines the evolution of
$\mathbb{E}|f|^{2}$, the expectation of a quadratic expression in $f$. Because
of the specific form of the latter, this involves a specific contraction of
the 4-tensor describing all covariations,

(23) $\displaystyle
d\mathbb{E}|f|^{2}=\mathbb{E}[\operatorname{\text{tr}}df^{*}df],$

which will be useful. Turning to (7), we note that the parity argument used
for (13) yields

$\displaystyle[\nabla
d\phi\,(\nabla\partial_{i}d\phi)^{*}]=0\quad\mbox{and}\quad[\nabla
d\phi\,(\nabla d\psi\otimes e_{i})^{*}]=0.$

Hence the r. h. s. of (23) splits

$\displaystyle
d\mathbb{E}|f|^{2}=\varepsilon^{2}\mathbb{E}[\operatorname{\text{tr}}(f\nabla
d\phi)^{*}(f\nabla d\phi)]$
$\displaystyle+\varepsilon^{2}\mathbb{E}[\operatorname{\text{tr}}((\tilde{\phi}^{i}a-\tilde{\sigma}^{i}J)\nabla\partial_{i}d\phi-{\tilde{\phi}^{i}}J\nabla
d\psi\otimes{e_{i}})^{*}((\tilde{\phi}^{j}a-\tilde{\sigma}^{j}J)\nabla\partial_{j}d\phi-{\tilde{\phi}^{j}}J\nabla
d\psi\otimes{e_{j}})].$

By Cauchy-Schwarz applied to the inner product $\mathbb{E}[{\rm tr}f^{*}g]$,
we may further split the last term:

$\displaystyle d\mathbb{E}|f|^{2}$
$\displaystyle\leq\varepsilon^{2}\mathbb{E}[\operatorname{\text{tr}}(f\nabla
d\phi)^{*}(f\nabla
d\phi)]+2\varepsilon^{2}\mathbb{E}[\operatorname{\text{tr}}((\tilde{\phi}^{i}a-\tilde{\sigma}^{i}J)\nabla\partial_{i}d\phi)^{*}((\tilde{\phi}^{j}a-\tilde{\sigma}^{j}J)\nabla\partial_{j}d\phi)]$
$\displaystyle+2\varepsilon^{2}\mathbb{E}[\operatorname{\text{tr}}(J\nabla
d\psi\otimes\tilde{\phi})^{*}(J\nabla d\psi\otimes\tilde{\phi})].$

Using the algebraic identities ${\rm tr}(ff^{\prime})^{*}ff^{\prime}={\rm
tr}f^{*}ff^{\prime}f^{\prime*}$, $(\alpha^{i}{\rm
id}+\beta^{i}J)^{*}(\alpha^{j}{\rm id}+\beta^{j}J)$
$=(\alpha^{i}\alpha^{j}+\beta^{i}\beta^{j}){\rm id}$
$+(\alpha^{i}\beta^{j}-\alpha^{j}\beta^{i})J$, and ${\rm
tr}(\xi\otimes\phi)^{*}(\xi\otimes\phi)$ $=|\xi|^{2}|\phi|^{2}$, we isolate
the quadratic variations of the drivers:

$\displaystyle d\mathbb{E}|f|^{2}$
$\displaystyle\leq\varepsilon^{2}\mathbb{E}\operatorname{\text{tr}}f^{*}f[\nabla
d\phi(\nabla d\phi)^{*}]$
$\displaystyle+2\varepsilon^{2}\mathbb{E}\big{(}\tilde{\phi}^{i}\tilde{\phi}^{j}+(\varepsilon\psi\tilde{\phi}^{i}-\tilde{\sigma}^{i})(\varepsilon\psi\tilde{\phi}^{j}-\tilde{\sigma}^{j})\big{)}{\rm
tr}[\nabla\partial_{j}d\phi(\nabla\partial_{i}d\phi)^{*}]+2\varepsilon^{2}\mathbb{E}|\tilde{\phi}|^{2}[\nabla
d\psi\cdot\nabla d\psi],$

where in the second r. h. s. term, the skew-symmetric part in $(i,j)$ drops
out because $[\nabla\partial_{j}d\phi(\nabla\partial_{i}d\phi)^{*}]$ is
symmetric in $(i,j)$ by stationarity888As a bilinear form the covariation
satisfies a Leibniz rule. Since the covariation under consideration is
deterministic and by stationarity constant in space, this implies that we can
integrate by parts.. By isotropy, the symmetric endomorphism $[\nabla
d\phi(\nabla d\phi)^{*}]$ is a multiple of the identity; likewise, ${\rm
tr}[\nabla\partial_{j}d\phi(\nabla\partial_{i}d\phi)^{*}]$ is a multiple of
$\delta_{ij}$. Therefore, and using

(24)
$2(|\tilde{\phi}|^{2}+|\varepsilon\psi\tilde{\phi}-\tilde{\sigma}|^{2})=|\tilde{\phi}\otimes
a-\tilde{\sigma}\otimes J|^{2},$

the above inequality collapses to

$\displaystyle d\mathbb{E}|f|^{2}$
$\displaystyle\leq\frac{\varepsilon^{2}}{2}\mathbb{E}|f|^{2}[{\rm tr}\nabla
d\phi(\nabla d\phi)^{*}]$
$\displaystyle+\frac{\varepsilon^{2}}{2}\mathbb{E}|\tilde{\phi}\otimes
a-\tilde{\sigma}\otimes J|^{2}{\textstyle\sum_{i=1}^{2}}[{\rm
tr}\nabla\partial_{i}d\phi(\nabla\partial_{i}d\phi)^{*}]+2\varepsilon^{2}\mathbb{E}|\tilde{\phi}|^{2}[\nabla
d\psi\cdot\nabla d\psi].$

Hence by (1) in form of

(25) $\displaystyle d\tilde{\lambda}^{2}=\varepsilon^{2}ds,$

(22) follows once we identify the following quadratic variations

(26) $\displaystyle\tilde{\lambda}^{2}[{\rm tr}\nabla d\phi\,(\nabla
d\phi)^{*}]$ $\displaystyle=ds,$ (27)
$\displaystyle\tilde{\lambda}^{2}L^{2}{\textstyle\sum_{i=1}^{2}}[{\rm
tr}\nabla\partial_{i}d\phi\,(\nabla\partial_{i}d\phi)^{*}]$
$\displaystyle=ds,$ (28) $\displaystyle L^{2}[\nabla d\psi\cdot\nabla d\psi]$
$\displaystyle=ds.$

We start by deriving (26) from (20). By stationarity we have $[{\rm tr}\nabla
d\phi(\nabla d\phi)^{*}]$ $=[d\phi\cdot(-\triangle)d\phi]$, into which we
insert (17) to the effect of999recall that $\xi.\phi$ denotes the canonical
pairing of a cotangent vector $\xi$ and a tangent vector $\phi$
$\tilde{\lambda}[{\rm tr}\nabla d\phi(\nabla d\phi)^{*}]$ $=[J\nabla
d\psi.d\phi]$. Once more by stationarity this turns into $\tilde{\lambda}[{\rm
tr}\nabla d\phi\,(\nabla d\phi)^{*}]$ $={\rm tr}J[d\psi\,\nabla d\phi]$, so
that it remains to appeal to (20) and $J^{2}=-{\rm id}$.

We now argue that (28) follows from (19). Indeed, by stationarity we have
$[\nabla d\psi\cdot\nabla d\psi]=[d\psi\,(-\triangle)d\psi]$. By definition of
the increment $d\psi$, it has values in the space of stationary fields which
are Fourier supported on $L|k|=1$, to the effect of
$L^{2}(-\triangle)d\psi=d\psi$. Hence we obtain the desired

(29) $\displaystyle L^{2}[\nabla d\psi\cdot\nabla d\psi]=[d\psi\,d\psi].$

Finally, we deduce (27) from (26) by a similar argument: By stationarity, we
have $\sum_{i=1}^{2}$ $[{\rm
tr}\nabla\partial_{i}d\phi(\nabla\partial_{i}d\phi)^{*}]$ $=[{\rm tr}\nabla
d\phi\,(-\triangle\nabla d\phi)^{*}]$, the Fourier support of $d\psi$
transmits to $\nabla d\phi$ via the defining equation (2). ∎

To estimate $\mathbb{E}|f|^{2}$, one needs the estimates on
$\mathbb{E}|\tilde{\phi}|^{2}$, $\mathbb{E}|\tilde{\sigma}|^{2}$ and
$\mathbb{E}|\psi\tilde{\phi}|^{2}$ that are stated in [2, Lemma 4]. We now
give their proofs in the continuum setting.

###### Proof of [2, Lemma 4] revisited..

Step 1: Evolution of $\mathbb{E}|\tilde{\phi}|^{2}$. We claim that

(30)
$d\mathbb{E}|\tilde{\phi}|^{2}-\varepsilon^{2}\mathbb{E}|\tilde{\phi}|^{2}\frac{ds}{2\tilde{\lambda}^{2}}=\varepsilon^{2}\frac{L^{2}ds}{\tilde{\lambda}^{2}},$

which integrates to

(31)
$\mathbb{E}|\tilde{\phi}|^{2}\lesssim\frac{\varepsilon^{2}L^{2}}{\tilde{\lambda}^{2}}.$

Let us first argue how (30) implies (31): We introduce the variables
$a:=\mathbb{E}|\tilde{\phi}|^{2}$ and $x:=\tilde{\lambda}^{2}$ so that
$L^{2}=\exp(\frac{2}{\varepsilon^{2}}(x-1))$,
$d\tilde{\lambda}^{2}=\varepsilon^{2}ds$ and (30) reads

$\sqrt{x}\frac{d}{dx}\frac{a}{\sqrt{x}}=\frac{da}{dx}-\frac{1}{2x}a=\frac{\exp(\frac{2}{\varepsilon^{2}}(x-1))}{x},$

which, since $a(x=1)=0$, implies

$a=\frac{\exp(\frac{2}{\varepsilon^{2}}(x-1))}{x}\int_{1}^{x}dy\,(\frac{x}{y})^{\frac{3}{2}}\exp(\frac{2}{\varepsilon^{2}}(y-x))\lesssim\varepsilon^{2}\frac{\exp(\frac{2}{\varepsilon^{2}}(x-1))}{x}.$

We now argue that (30) holds. Since $\tilde{\phi}$ is a martingale we have

$\displaystyle d\mathbb{E}|\tilde{\phi}|^{2}=\mathbb{E}[d\tilde{\phi}\cdot
d\tilde{\phi}].$

Inserting (3), the same parity argument that led to (13) yields the splitting

(32) $\displaystyle
d\mathbb{E}|\tilde{\phi}|^{2}=\varepsilon^{2}\mathbb{E}[d\phi\cdot
d\phi]+\varepsilon^{2}\mathbb{E}\tilde{\phi}^{i}\tilde{\phi}^{j}[\partial_{i}d\phi\cdot\partial_{j}d\phi].$

Turning to the quadratic variation of the drivers, we have by isotropy

(33) $\displaystyle
2\tilde{\lambda}^{2}[\partial_{i}d\phi\cdot\partial_{j}d\phi]=\delta_{ij}\tilde{\lambda}^{2}[{\rm
tr}\nabla d\phi(\nabla
d\phi)^{*}]\stackrel{{\scriptstyle(\ref{wr18})}}{{=}}\delta_{ij}ds.$

By the same argument that implied (29) one infers

(34) $\displaystyle\tilde{\lambda}^{2}[d\phi\cdot
d\phi]=L^{2}\tilde{\lambda}^{2}[{\rm tr}\nabla d\phi(\nabla
d\phi)^{*}]\stackrel{{\scriptstyle(\ref{wr18})}}{{=}}L^{2}ds.$

Inserting these two identities into (32) and appealing to (25) gives (30).

Step 2: Evolution of $\mathbb{E}|\tilde{\sigma}|^{2}$. We claim that

(35)
$d\mathbb{E}|\tilde{\sigma}|^{2}=\varepsilon^{2}\mathbb{E}|\tilde{\sigma}|^{2}\frac{ds}{2\tilde{\lambda}^{2}}+2\varepsilon^{2}\mathbb{E}\tilde{\sigma}\cdot
J\tilde{\phi}\frac{ds}{\tilde{\lambda}}+\varepsilon^{2}L^{2}ds+\varepsilon^{2}\mathbb{E}|\phi|^{2}ds,$

which upon integration becomes

(36) $\mathbb{E}|\tilde{\sigma}|^{2}\lesssim\varepsilon^{2}L^{2}.$

Indeed, using Young’s inequality and (31) on (35) implies

$d\mathbb{E}|\tilde{\sigma}|^{2}-\varepsilon^{2}\mathbb{E}|\tilde{\sigma}|^{2}\frac{ds}{\tilde{\lambda}^{2}}\lesssim\varepsilon^{2}L^{2}ds,$

which in terms of $b:=\mathbb{E}|\tilde{\sigma}|^{2}$ and (the previously
introduced) $x=\tilde{\lambda}^{2}$ may be rewritten as

$x\frac{d}{dx}\frac{b}{x}=\frac{db}{dx}-\frac{b}{x}\lesssim\exp(\frac{2}{\varepsilon^{2}}(x-1)),$

so that

$b\lesssim\exp(\frac{2}{\varepsilon^{2}}(x-1))\int_{1}^{x}dy\,\frac{x}{y}\exp(\frac{2}{\varepsilon^{2}}(y-x)).$

(36) follows from the last estimate.

Now comes the argument for (35). Note that by Itô’s formula, we have

$d|\tilde{\sigma}|^{2}=2\tilde{\sigma}\cdot
d\tilde{\sigma}+[d\tilde{\sigma}\cdot d\tilde{\sigma}]$

so that by (4), in which we substitute (25), we obtain

$d\mathbb{E}|\tilde{\sigma}|^{2}=2\varepsilon^{2}\mathbb{E}\tilde{\sigma}\cdot
J\tilde{\phi}\frac{ds}{\tilde{\lambda}}+\mathbb{E}[d\tilde{\sigma}\cdot
d\tilde{\sigma}].$

Let us observe that a parity argument similar to (13) implies

$[d\tilde{\sigma}\cdot
d\tilde{\sigma}]=\varepsilon^{2}\tilde{\sigma}^{i}\tilde{\sigma}^{j}[\partial_{i}d\phi\cdot\partial_{j}d\phi]+2\varepsilon^{2}\tilde{\sigma}^{i}\tilde{\phi}\cdot[d\psi\,\partial_{i}d\phi]+\varepsilon^{2}[(d\sigma+d\psi\,\tilde{\phi})\cdot(d\sigma+d\psi\,\tilde{\phi})].$

Using (20), (25) and (33) this equation simplifies to

$[d\tilde{\sigma}\cdot
d\tilde{\sigma}]=\varepsilon^{2}|\tilde{\sigma}|^{2}\frac{ds}{2\tilde{\lambda}^{2}}+\varepsilon^{2}\tilde{\sigma}\cdot
J\tilde{\phi}\frac{ds}{\tilde{\lambda}}+\varepsilon^{2}[d\sigma\cdot
d\sigma]+2\varepsilon^{2}\tilde{\phi}\cdot[d\psi\,d\sigma]+\varepsilon^{2}|\tilde{\phi}|^{2}[d\psi\,d\psi].$

Since $\tilde{\phi}$ has vanishing expectation, the mixed term
$\tilde{\phi}\cdot[d\psi\,d\sigma]$ has vanishing expectation, so that overall

$d\mathbb{E}|\tilde{\sigma}|^{2}=\varepsilon^{2}\mathbb{E}|\tilde{\sigma}|^{2}\frac{ds}{2\tilde{\lambda}^{2}}+2\varepsilon^{2}\mathbb{E}\tilde{\sigma}\cdot
J\tilde{\phi}\frac{ds}{\tilde{\lambda}}+\varepsilon^{2}\mathbb{E}[d\sigma\cdot
d\sigma]+\varepsilon^{2}\mathbb{E}|\tilde{\phi}|^{2}ds.$

Below we will argue that

(37) $\varepsilon^{2}[d\sigma\cdot d\sigma]=\varepsilon^{2}L^{2}ds,$

which together with the previous equation implies (35).

We now give the argument for (37). Via an integration by parts we obtain
$\mathbb{E}\nabla d\phi^{i}\cdot J\nabla d\sigma^{i}=0$ so that in particular
$\operatorname{\text{tr}}\mathbb{E}\nabla d\phi(J\nabla d\sigma)^{*}=0$, which
means that $\nabla d\phi$ and $J\nabla d\sigma$ are uncorrelated. Hence, (2)
provides an orthogonal decomposition of $d\psi J$. Therefore, and since all of
the following expressions are deterministic, we may conclude

$[\operatorname{\text{tr}}d\psi J(d\psi
J)^{*}]=\tilde{\lambda}^{2}[\operatorname{\text{tr}}\nabla d\phi(\nabla
d\phi)^{*}]+[\operatorname{\text{tr}}J\nabla d\sigma(J\nabla d\sigma)^{*}].$

In view of (19), and (26) this implies

$[\operatorname{\text{tr}}\nabla d\sigma(\nabla
d\sigma)^{*}]=2[d\psi\,d\psi]-\tilde{\lambda}^{2}[\operatorname{\text{tr}}\nabla
d\phi(\nabla d\phi)^{*}]=ds$

The argument that leads to (29) now implies (37).

Step 3: Evolution of $\mathbb{E}|\tilde{\phi}|^{4}$. We claim that

(38)
$d\mathbb{E}|\tilde{\phi}|^{4}\lesssim\varepsilon^{2}\mathbb{E}|\tilde{\phi}|^{4}\frac{ds}{\tilde{\lambda}^{2}}+\varepsilon^{2}\mathbb{E}|\tilde{\phi}|^{2}\frac{L^{2}ds}{\tilde{\lambda}^{2}},$

which integrates to

(39)
$\mathbb{E}|\tilde{\phi}|^{4}\lesssim\varepsilon^{4}\frac{L^{4}}{\tilde{\lambda}^{4}}.$

The integration follows along the lines of (36).

Itô’s formula implies

$d|\tilde{\phi}|^{2}=2\tilde{\phi}\cdot
d\tilde{\phi}+[\tilde{\phi}\cdot\tilde{\phi}],$

which we combine with (3) and (33), (34) on the quadratic variation to obtain

$d|\tilde{\phi}|^{2}=2\varepsilon\tilde{\phi}\cdot
d\phi+2\varepsilon\tilde{\phi}^{i}\tilde{\phi}\cdot\partial_{i}d\phi+\varepsilon^{2}\frac{L^{2}ds}{\tilde{\lambda}^{2}}+\varepsilon^{2}|\tilde{\phi}|^{2}\frac{ds}{2\tilde{\lambda}^{2}}.$

Since the covariation of $\phi$ and $\partial_{i}\phi$ vanishes, we obtain

$\displaystyle d|\tilde{\phi}|^{4}$
$\displaystyle=\text{martingale}+\varepsilon^{2}|\tilde{\phi}|^{4}\frac{ds}{\tilde{\lambda}^{2}}+2\varepsilon^{2}|\tilde{\phi}|^{2}\frac{L^{2}ds}{\tilde{\lambda}^{2}}$
$\displaystyle+4\varepsilon^{2}[(\tilde{\phi}\cdot d\phi)(\tilde{\phi}\cdot
d\phi)]+4\varepsilon^{2}[(\tilde{\phi}^{i}\tilde{\phi}\cdot\partial_{i}d\phi)(\tilde{\phi}^{j}\tilde{\phi}\cdot\partial_{j}d\phi)],$

which implies

$\displaystyle d\mathbb{E}|\tilde{\phi}|^{4}$
$\displaystyle=\varepsilon^{2}\mathbb{E}|\tilde{\phi}|^{4}\frac{ds}{\tilde{\lambda}^{2}}+2\varepsilon^{2}\mathbb{E}|\tilde{\phi}|^{2}\frac{L^{2}ds}{\tilde{\lambda}^{2}}$
$\displaystyle+4\varepsilon^{2}\mathbb{E}[(\tilde{\phi}\cdot
d\phi)(\tilde{\phi}\cdot
d\phi)]+4\varepsilon^{2}\mathbb{E}[(\tilde{\phi}^{i}\tilde{\phi}\cdot\partial_{i}d\phi)(\tilde{\phi}^{j}\tilde{\phi}\cdot\partial_{j}d\phi)].$

For the third term, we use (25) and (33), which yields

$4\varepsilon^{2}\mathbb{E}[(\tilde{\phi}\cdot d\phi)(\tilde{\phi}\cdot
d\phi)]=2\varepsilon^{2}\mathbb{E}|\tilde{\phi}|^{2}\frac{L^{2}ds}{\tilde{\lambda}^{2}}.$

For the fourth term, we note that the full covariation of first derivatives
can be estimated by Cauchy-Schwarz, i. e.

$\mathbb{E}[\partial_{i}d\phi^{k}\,\partial_{j}d\phi^{l}]\leq\mathbb{E}\operatorname{\text{tr}}[(\nabla
d\phi)^{*}(\nabla d\phi)].$

By a contraction we obtain the estimate

$\mathbb{E}[(\tilde{\phi}^{i}\tilde{\phi}\cdot\partial_{i}d\phi)(\tilde{\phi}^{j}\tilde{\phi}\cdot\partial_{j}d\phi)]\lesssim\mathbb{E}|\tilde{\phi}|^{4}\operatorname{\text{tr}}[(\nabla
d\phi)^{*}(\nabla d\phi)].$

By using (26) this yields (38). ∎

We collected all ingredients to estimate $\mathbb{E}|f|^{2}$, i. e. provide an
alternative proof of [2, Lemma 5].

###### Proof of Proposition 1.1..

We are now ready to prove

(40) $\mathbb{E}|f|^{2}\lesssim\varepsilon^{2}\tilde{\lambda}$

in the continuum setting. Indeed, (22), together with (24) and (31), implies

$\displaystyle
d\mathbb{E}|f|^{2}-\varepsilon^{2}\mathbb{E}|f|^{2}\frac{ds}{2\tilde{\lambda}^{2}}\lesssim\varepsilon^{2}\mathbb{E}|\varepsilon\psi\tilde{\phi}-\tilde{\sigma}|^{2}\frac{ds}{L^{2}\tilde{\lambda}^{2}}+\varepsilon^{4}\frac{ds}{\tilde{\lambda}^{2}}.$

By using Cauchy-Schwarz on the product $|\tilde{\phi}\psi|^{2}$ together with
(36) and (39) the above becomes

$d\mathbb{E}|f|^{2}-\varepsilon^{2}\mathbb{E}|f|^{2}\frac{ds}{2\tilde{\lambda}^{2}}\lesssim\varepsilon^{6}(\mathbb{E}\psi^{4})^{\frac{1}{2}}\frac{ds}{\tilde{\lambda}^{4}}+\varepsilon^{4}\frac{ds}{\tilde{\lambda}^{2}}.$

The Gaussianity of $\psi$ implies
$(\mathbb{E}\psi^{4})^{\frac{1}{2}}\sim\mathbb{E}\psi^{2}=\log L$, which
together with (1) in form of $\tilde{\lambda}^{2}=1+\varepsilon^{2}\log L$,
shows that the last term on the r. h. s. is the dominant term so that

$d\mathbb{E}|f|^{2}-\varepsilon^{2}\mathbb{E}|f|^{2}\frac{ds}{2\tilde{\lambda}^{2}}\lesssim\varepsilon^{4}\frac{ds}{\tilde{\lambda}^{2}}.$

In terms of $c:=\mathbb{E}|f|^{2}$ and $x=\tilde{\lambda}^{2}$ this may be
rewritten in the form

$\sqrt{x}\frac{d}{dx}\frac{c}{\sqrt{x}}\lesssim\varepsilon^{2}\frac{1}{x}$

so that

$c\lesssim\varepsilon^{2}\sqrt{x}\int_{1}^{x}dy\,\frac{1}{y^{\frac{3}{2}}}.$

This is nothing but (40) in terms of $c$ and $x$. ∎

## Acknowledgements

The authors thank Lihan Wang for fruitful discussions in earlier stages of
this work.

## References

  * [1] N. Barashkov and M. Gubinelli.  A variational method for $\Phi^{4}_{3}$.  _Duke Math. J._ 169 (17) 3339 - 3415, 15 November 2020.
  * [2] G. Chatzigeorgiou, P. Morfe, F. Otto, and L. Wang.  The Gaussian free-field as a stream function: asymptotics of effective diffusivity in infra-red cut-off.  _arXiv preprint_ arXiv:2212.14244 (2022).
  * [3] M. Gubinelli, and S.-J. Meyer.  The FBSDE approach to sine-Gordon up to $6\pi$.  _arXiv preprint_ arXiv:2401.13648 (2024).


# Development of a 127Xe calibration source for nEXO

B. G. Lenardo  C. A. Hardy  R. H. M. Tsang  J. C. Nzobadila Ondze  A. Piepke
S. Triambak  A. Jamil  G. Adhikari  S. Al Kharusi  E. Angelico  I. J. Arnquist
V. Belov  E. P. Bernard  A. Bhat  T. Bhatta  A. Bolotnikov  P. A. Breur  J. P.
Brodsky  E. Brown  T. Brunner  E. Caden  G. F. Cao 111Also at: University of
Chinese Academy of Sciences, Beijing, China L. Cao  B. Chana  S. A. Charlebois
D. Chernyak  M. Chiu  J. R. Cohen  R. Collister  J. Dalmasson  T. Daniels  L.
Darroch  R. DeVoe  M. L. di Vacri  Y. Y. Ding  M. J. Dolinski  J. Echevers  B.
Eckert  M. Elbeltagi  L. Fabris  D. Fairbank  W. Fairbank  J. Farine  Y. S. Fu
G. Gallina  P. Gautam  G. Giacomini  W. Gillis  C. Gingras  R. Gornea  G.
Gratta  K. Harouaka  M. Heffner  E. Hein  J. Hößl  A. House  A. Iverson  X. S.
Jiang  A. Karelin  L. J. Kaufman  R. Krücken  A. Kuchenkov  K. S. Kumar  A.
Larson  K. G. Leach  D. S. Leonard  G. Li  S. Li  Z. Li  C. Licciardi  R.
Lindsay  R. MacLellan  J. Masbou  K. McMichael  M. Medina Peregrina  B. Mong
D. C. Moore  K. Murray  J. Nattress  C. R. Natzke  X. E. Ngwadla  K. Ni  Z.
Ning  J. L. Orrell  G. S. Ortega  I. Ostrovskiy  C. T. Overman  A. Perna  T.
Pinto Franco  A. Pocar  J. F. Pratte  N. Priel  E. Raguzin  G. J. Ramonnye  H.
Rasiwala  K. Raymond  G. Richardson  M. Richman  J. Ringuette  P. C. Rowson
R. Saldanha  S. Sangiorgio  X. Shang  A. K. Soma  F. Spadoni  V. Stekhanov  X.
L. Sun  S. Thibado  A. Tidball  J. Todd  T. Totev  O. A. Tyuka  F. Vachon  V.
Veeraraghavan  S. Viel  K. Wamba  Y. Wang  Q. Wang  W. Wei  L. J. Wen  U.
Wichoski  S. Wilde  W. H. Wu  W. Yan  L. Yang  O. Zeldovich  J. Zhao  T.
Ziegler

###### Abstract

We study a possible calibration technique for the nEXO experiment using a
127Xe electron capture source. nEXO is a next-generation search for
neutrinoless double beta decay ($0\nu\beta\beta$) that will use a 5-tonne,
monolithic liquid xenon time projection chamber (TPC). The xenon, used both as
source and detection medium, will be enriched to 90% in 136Xe. To optimize the
event reconstruction and energy resolution, calibrations are needed to map the
position- and time-dependent detector response. The
$36.3\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}$ half-life of 127Xe and its small
$Q$-value compared to that of 136Xe $0\nu\beta\beta$ would allow a small
activity to be maintained continuously in the detector during normal
operations without introducing additional backgrounds, thereby enabling _in-
situ_ calibration and monitoring of the detector response. In this work we
describe a process for producing the source and preliminary experimental
tests. We then use simulations to project the precision with which such a
source could calibrate spatial corrections to the light and charge response of
the nEXO TPC.

##  1 Introduction

The nEXO experiment is a planned tonne-scale search for neutrinoless double
beta decay ($0\nu\beta\beta$) in 136Xe using a cylindrical liquid-phase time
projection chamber (TPC) [1]. The low intrinsic background, 3-dimensional
event reconstruction, and powerful self-shielding of the active liquid xenon
volume will enable nEXO to achieve the ultra-low backgrounds needed to reach a
sensitivity to $0\nu\beta\beta$ beyond a half-life of $10^{28}$ years [2]. The
TPC provides a dual-channel measurement of interactions in the liquid xenon
target: a scintillation signal, detected promptly via photosensors around the
barrel of the detector; and an ionization signal, detected by applying a
uniform electric field across the TPC to drift the charge to a collection
plane at the anode. The two signals provide complementary information that can
be combined to enable the reconstruction of both the 3-dimensional position of
each interaction vertex and the deposited energy. One of the performance
targets for nEXO is to achieve an energy resolution better than $\sigma/E=1$%
at the $0\nu\beta\beta$ Q-value (2.457 MeV), which contributes to the
rejection of backgrounds, in particular allowing for the separation of a
$0\nu\beta\beta$ signal from the endpoint of the two-neutrino double beta
decay ($2\nu\beta\beta$) spectrum.

The scintillation and ionization signals are strongly anticorrelated due to
large fluctuations in electron-ion recombination [3, 4], meaning the energy is
optimally reconstructed as a linear combination of the two signals. This can
be expressed as

$E=W\cdot(S_{0}+Q_{0}),$ (1.1)

where $S_{0}$ ($Q_{0}$) is the number of scintillation photons (ionization
electrons) released by the event after recombination, and $W$ is a
proportionality constant which represents the average energy required to
produce a single quantum of either light or charge. Assuming perfect linearity
and perfect anticorrelation between light and charge (that is, each electron-
ion pair which recombines produces one scintillation photon), $W$ is field-
and energy-independent; these assumptions are supported by measurements in
Ref. [5] and Ref. [6], respectively. In this picture, the energy resolution is
defined by the intrinsic fluctuations in $(S_{0}+Q_{0})$ combined with the
sources of fluctuations in the reconstruction of $S_{0}$ and $Q_{0}$ from
measured quantities.

An important source of resolution broadening stems from the position-dependent
detection efficiencies for both scintillation light and ionized charge. For
the scintillation signals, position dependence arises from a combination of
geometrical effects, surface reflectivities, and the detection efficiency of
the photosensors. $S_{0}$ is reconstructed as

$S_{0}=\frac{S_{\text{meas}}}{\epsilon_{QE}\,\times\,\epsilon_{{LM}}(x,y,z)},$
(1.2)

where $S_{\text{meas}}$ is the measured number of scintillation photons,
$\epsilon_{QE}$ is the photosensor quantum efficiency, and $\epsilon_{{LM}}$
is the “lightmap”, a 3-dimensional function that describes the photon
transport efficiency in the TPC. The lightmap is a function of detector
properties and may vary on timescales of several months [7]. For the
ionization signals, position dependence is driven primarily by electron
attachment on electronegative impurities in the liquid xenon, which can be
modeled as an exponential attenuation of the charge signal as a function of
the drift time $t$. Then $Q_{0}$ is reconstructed as

$Q_{0}=\frac{Q_{\text{meas}}}{e^{-t/\tau_{e}}},$ (1.3)

where $Q_{\text{meas}}$ is the measured number of ionization electrons, and
$\tau_{e}$ (known as the “electron lifetime”) represents the average time a
free electron can drift in the liquid before attaching to an impurity. By
design, the drift field in the TPC’s active volume is uniform, meaning $t$ is
proportional to the $z$ position of the event, i.e. $t=z/v_{d}$, where $v_{d}$
is the drift velocity.***Here we assume that the electric field is constant
throughout the active region of the TPC. Non-uniformities in the field would
create non-uniforimities in both the recombination fraction and the drift
velocity $v_{d}$, which would introduce additional second-order position
dependencies in the detected light and charge signals. The purity of the
liquid xenon is highly correlated with operation of the recirculation system
and can vary on timescales of $O(1)$ day [8]. Optimizing the energy
reconstruction relies on the ability to measure both $\epsilon_{LM}$ and
$\tau_{e}$ with appropriate frequency and accuracy using calibration data.

The simplest technique for measuring these two quantities is to use a source
of ionizing radiation to create events of fixed energy throughout the TPC. A
standard technique is to use $\gamma$-ray sources positioned next to the
detector. While this is the baseline method for nEXO, the powerful self-
shielding of liquid xenon means that, to achieve sufficient statistics in the
center of the TPC, the readout system needs to cope with high rates and pileup
at the edges and may require long calibration campaigns. An alternative
strategy is to use radioisotopes that can be mixed directly into the active
liquid xenon. Several such sources have been used previously for calibrating
liquid xenon detectors, namely 83mKr [9, 10, 11], tritiated methane [12], 37Ar
[13], the neutron-activated isomers 129mXe and 131mXe [14], and 220Rn [15,
16]. However, the first five produce signals at or below nEXO’s $\sim$200 keV
trigger threshold. While 220Rn is indeed under consideration for use in nEXO,
a) the short half-life of the decay chain (dominated by the
$10.6\text{\,}\mathrm{h}\mathrm{r}$ half-life of 212Pb) is comparable to the
xenon recirculation time and may affect the uniformity with which it can
distribute through the TPC, and b) 208Tl $\beta$-decay in the 220Rn decay
chain interferes with the $0\nu\beta\beta$ $Q$-value, limiting the frequency
with which that source could be used. It is therefore of interest to identify
other radioisotopes that can both mix throughout the TPC and produce
monoenergetic signals above nEXO’s trigger threshold.

Here we study the use of 127Xe as an injected calibration source for nEXO.
This isotope decays to 127I via electron capture (EC), predominantly releasing
a total of either 236 keV or 408 keV of ionizing radiation that can be used
for calibration. The full decay scheme is shown in Figure ‡ ‣ 1. While this
isotope has been used previously to characterize the scintillation and
ionization yields of liquid xenon [17, 18], here we specifically study its use
for calibrating position-dependent detection efficiencies in large-scale
detectors. The 36-day half-life will ensure that the source has sufficient
time to mix uniformly throughout the TPC, and the 662.3(20) keV $Q$-value
ensures that these events do not produce enough energy to interfere with the
$0\nu\beta\beta$ search. We consider a calibration strategy in which a
constant activity of $\sim$1 Bq is maintained continuously in the nEXO TPC via
frequent, controlled injections of 127Xe into the xenon recirculation loop.
This activity is similar to the expected overall background rate in nEXO (the
$2\nu\beta\beta$ decay of 136Xe alone produces 0.2 Bq), and will add
negligible dead time. Calibration data could then be taken concurrently with
physics data, providing quasi-real-time information on the detector response.

Figure 1: Decay scheme of 127Xe, from [19].‡‡‡Figure produced using the
Laraweb tool: http://www.nucleide.org/Laraweb/index.php127Xe decays through
electron capture to 127I, primarily populating the levels at either
$374.991(12)$ keV or $202.860(10)$ keV, which relax to the ground state via
$\gamma$-ray emission. Electron capture from the K-shell and L-shell have a
branching ratios of 84.2(8)% and 12.9(1)% and result in the emission of an
additional 33.045 keV or 5.185 keV in the form of x-rays and/or Auger
electrons. Most decays therefore deposit a total of either 408 keV or 236 keV
in the liquid xenon, with branching fractions of 40% and 44%, respectively.

First, we discuss production of the source via neutron activation at a
research reactor. Next, we demonstrate its use for measuring the electron
lifetime in a prototype liquid xenon detector. Finally, we describe
simulations of 127Xe decays with a detailed model of the nEXO TPC and estimate
the precision with which such a source could calibrate both the lightmap and
the electron lifetime for different integration times.

##  2 Source production

###  2.1 Procedure for neutron activation of ${}^{\text{nat}}$Xe

127Xe can be readily produced via neutron capture on 126Xe, which is present
in natural xenon at an isotopic abundance of $\sim$0.1%. To produce the
source, a 150 cm3 stainless steel (316L) sample cylinder was filled with
$69\text{\,}\mathrm{g}$ of natXe gas and shipped to the research reactor at
McClellan Nuclear Research Center (MNRC). The sample cylinder was placed in
the Neutron Transmutation Doping (NTD) void, for which the steady-state
neutron spectrum at 1 MW is given in Ref [20] and shown in Figure 2(a). The
irradiation was performed for 15 minutes at a power of
$\sim$$250\text{\,}\mathrm{k}\mathrm{W}$. The production of radioisotopes in
both the stainless steel cylinder and the natXe was calculated numerically by
folding the expected reactor spectrum with cross sections from standard
libraries. Neutron capture cross sections for most isotopes were taken from
ENDF/B-VII.0 [21]. For the two metastable isomers, 129mXe and 131mXe, we used
cross sections from the TENDL-2019 library [22], which are conveniently given
in terms of the total cross section for populating the desired state given a
specific reaction.

The neutron flux incident on the sample was calibrated using activation in the
stainless steel cylinder, which produces three long-lived products: 51Cr and
58Fe, which are produced primarily by thermal neutron capture, and 58Co, which
is produced primarily by fast neutron $(n,p)$ reactions. Two short ($\sim$30
min) radioassay measurements of the cylinder using high-purity germanium
(HPGe) counters were performed: the first was taken at MNRC 9 days after
irradiation; the second was taken at Stanford University 25 days after
irradiation. We estimate $\sim$50% systematic uncertainties in each
measurement due to uncertainties in the counting geometry. The measured
activities are listed in Table 1 and plotted in Figure 2(b). We infer the
overall neutron flux by scaling the neutron spectrum in the calculations
described above so that the predicted activities match the measured
activities. While the estimated fluxes from the measurements at 25 days are
systematically lower than those at 9 days by approximately 30%, both
measurements are consistent within the estimated uncertainties. Importantly,
the neutron fluxes inferred from the thermal-neutron-induced reactions are
consistent with those inferred from the fast neutron reactions, indicating
that the assumed spectral shape is adequate. Taking the average and standard
deviation of all inferred values gives an estimated neutron flux of $(2.1\pm
0.3)\times 10^{10}$ n/cm2/s. We note that, for unclear reasons, this flux is
lower by a factor of $\sim$5 than the expected steady-state flux at a reactor
power of $250\text{\,}\mathrm{k}\mathrm{W}$. We attribute this discrepancy to
operating the reactor in a transient mode during irradiation of our sample.

The predicted activities of radioisotopes produced in the natXe gas, given the
flux estimated above, are shown in Figure 2(c). The activity in the natXe gas
sample shortly after irradiation is dominated by short-lived isotopes. Of
these, the longest-lived are the metastable isomers 129mXe and 131mXe. Thermal
neutron capture is expected to be the dominant production mechanism for these
isotopes at a reactor, in contrast with the 252Cf-based activation scheme
reported in Ref. [14]. Despite the high initial activity, after a cool-off
period of $\sim$100 days the remaining activity is dominated by 127Xe. We note
that we predict a non-negligible amount of long-lived 137Cs produced by
neutron capture on 136Xe followed by beta decay of 137Xe ($T_{1/2}$ = 3.8
min), but Cs is expected to be easily removed from the Xe gas by standard
purification techniques prior to a deployment in nEXO.

(a)

(b)

(c)

Figure 2: The expected neutron spectrum and flux at MNRC [20] is shown in (a), along with neutron-induced activities of various radioisotopes in both the stainless steel cylinder (b) and the xenon gas (c). The overall scaling of the neutron spectrum is calibrated using two separate measurements of the stainless steel taken at MNRC (b, circles) and Stanford University (b, triangles), as explained in the text. The error bars are dominated by the estimated 50% systematic uncertainty in each measurement. This measured neutron flux is used to calculate the expected activity vs. time (solid lines).  Table 1: Radioisotopes produced by neutron activation of stainless steel, which we use to calibrate the neutron flux incident on the sample. The value of $\Delta T$ denotes the elapsed time between irradiation and when the measurements were made. Systematic uncertainties in the measurements are estimated to be $\pm$50%, due to uncertainties in the counting geometry. Statistical uncertainties are $O(1)$% and are omitted for clarity. Isotope | Half-life | Production mode |  $\gamma$-ray energy | $\Delta T$ | Meas. activity | Inferred flux
---|---|---|---|---|---|---
51Cr | 27.7 d | 50Cr $(n,\gamma)$ 51Cr | 320 keV | 9 d | 2960 kBq |  $2.3\times 10^{10}$ n/cm2/s
25 d | 1690 kBq |  $1.7\times 10^{10}$ n/cm2/s
59Fe | 44.5 d | 58Fe $(n,\gamma)$ 59Fe | 1099 & 1291 keV | 9 d | 61 kBq |  $2.7\times 10^{10}$ n/cm2/s
25 d | 40 kBq |  $2.3\times 10^{10}$ n/cm2/s
58Co | 70.9 d | 58Ni $(n,p)$ 58Co | 810 keV | 9 d | 11.3 kBq |  $2.3\times 10^{10}$ n/cm2/s
25 d | 8.1 kBq |  $1.8\times 10^{10}$ n/cm2/s

###  2.2 Low-background radioassay measurements of activated gas

High-precision radioassay measurements of the activated Xe gas were made by
transferring the gas into an unactivated cylinder and counting the sample with
a low-background HPGe detector at the University of Alabama [23]. A photo of
the experimental setup is shown in Figure 3(a). The purpose of these
measurements was twofold: first, we compare measured and predicted activities
in the gas to estimate the accuracy of our calculations, which requires more
sensitive measurements than those of the steel due to the smaller amount of
material; second, we use these data to search for any unexpected activation
products in the gas that could produce unwanted backgrounds in nEXO. Before
filling, the new cylinder was counted for approximately two weeks to obtain a
background spectrum. The cylinder was then filled and counted for two more
weeks. In each case, data were acquired in 4-hour intervals. The total energy
spectra, summed across the entire campaign, are shown in Figure 4.

We used a moving average window technique to detect peaks in the summed
spectrum, then fit them in each of the 4-hour time slices using a Gaussian
line shape plus a linear background. The fitted peak areas as a function of
time were then fitted to an exponential model to extract the mean lifetime of
the nuclides. The fitted energy and mean lifetime were compared to the
Evaluated Nuclear Structure Data Files (ENSDF, Ref. [24]) to identify a
candidate nuclide, then an activity of the nuclide in the $i$th time slice,
$A_{i}$, was determined by the following formula:

$A_{i}=\frac{C_{i}}{\varepsilon\cdot b\cdot\Delta t}$ (2.1)

where $C_{i}$ is the number of counts registered, $\Delta t=4$ hours, $b$ is
the tabulated $\gamma$-ray intensity from NNDC, and $\varepsilon$ is the
energy-dependent $\gamma$-ray detection efficiency. The efficiency estimation,
which includes modeling of the detector dead layer, was performed with
calibration sources and the Geant4-based GeSim package [23] (the GeSim
rendering of the counting setup is shown in Figure 3(b)). For point sources,
the systematic uncertainty in $\varepsilon$ is 9%. The dead-time and pileup in
the HPGe detector were negligible in these measurements. The measured
activities are shown in Table 2.

(a) Counting setup

(b) GEANT4 rendering

Figure 3: Activated Xe gas being counted at the University of Alabama. The
cylinder was placed directly below the HPGe detector, inside low-background
copper shielding. The GEANT4 model shown in (b) is used to calculate the
$\gamma$-ray detection efficiency of the counting geometry. Figure 4:
$\gamma$-ray energy spectrum of the activated Xe gas measured in a low-
background HPGe detector at the University of Alabama. The raw spectra (top)
consist of two measurements of an un-activated stainless steel sample
cylinder, the first under vacuum (black) and the second filled with the
activated Xe gas (blue). The former is subtracted from the latter to get the
$\gamma$-ray spectrum of the activated gas sample (bottom, red).  Table 2:
Measurements of long-lived ($T_{1/2}>1$ d) radioisotopes produced by neutron
activation of natXe. Activities are reported in terms of the activity
immediately following the irradiation campaign, as explained in the text. The
predicted activity uses the calibrated neutron flux obtained from measurements
of the stainless steel (Table 1). Statistical uncertainties on the measured
activities are $O(0.1)$% and are omitted for clarity. Systematic uncertainties
in the measurements are $O(10)$% and are discussed in the text (Section 2.2).

Isotope | Half-life | Production mode |  $\gamma$-ray energy | Measured activity | Predicted activity
---|---|---|---|---|---
127Xe | 36.3 d | 126Xe $(n,\gamma)$ 127Xe | 145 keV | 4.8 kBq | 6.5 kBq
172 keV | 5.0 kBq
203 keV | 5.5 kBq
317 keV* | -
375 keV | 6.3 kBq
406 keV* | -
578 keV* | -
618 keV | 6.9 kBq
129mXe | 8.88 d |  128Xe $(n,\gamma)$ 129mXe | 197 keV | 137 kBq | 36 kBq
129Xe $(n,n^{\prime})$ 129mXe
131mXe | 11.8 d |  130Xe $(n,\gamma)$ 131mXe | 164 keV | 187 kBq | 113 kBq
131Xe $(n,n^{\prime})$ 131mXe
133Xe | 5.25 d | 132Xe $(n,\gamma)$ 133Xe | 303 keV | 2170 kBq | 2190 kBq
384 keV | 2250 kBq
137Cs | 30.1 y |  136Xe $(n,\gamma)$ 137Xe, | 662 keV | $3.1\times 10^{-4}$ kBq | $3.4\times 10^{-2}$ kBq
137Xe $\rightarrow$ 137Cs

  * *

Features at these energies are produced by pile-up of lower-energy
$\gamma$-rays.

The measured activities are generally in agreement with the predictions, based
on the neutron flux evaluation described previously. For 127Xe, the measured
activity using different emission lines spans 4.8 – 6.9 kBq, in agreement with
the predicted value at the level of $\sim$30%. The measured values are
systematically lower than the predictions when the gamma energy goes below
$\sim$200 keV; we attribute this to additional systematic error in the
estimated $\gamma$-ray detection efficiency (beyond the expected 9%) due to
simplifications in the simulated geometry of the gas cylinder, which will
affect efficiency estimates for lower-energy $\gamma$-rays more strongly due
to their higher attenuation in materials. The activities derived from the
lines at 375 and 618 keV, for which the efficiencies are easier to estimate,
agree with the predicted value within 5%. Similar agreement is observed for
133Xe. Larger discrepancies are observed for the metastable isomers 129mXe and
131mXe, which we attribute to uncertainties in the evaluated reaction cross
sections for populating specific excited states; evaluations are extrapolated
from a single measurement of 129mXe and 131mXe activation by thermal neutrons
which carries 30–40% uncertainties [25]. The 137Cs activity exhibits the
largest discrepancy, with the measured activity a factor of $\sim$100 smaller
than the predicted value. We hypothesize that most of the Cs attaches to the
inner surface of the activated cylinder and is consequently not transferred to
the new cylinder used for counting. We take this as an encouraging sign that,
using dedicated purification techniques, the 137Cs contamination of such a
source could be reduced to negligible levels. We find no evidence for the
production of unexpected isotopes in the radioassay data. While we note that
additional measurements may be required to ensure that such a source meets the
stringent ultra-low-background requirements for nEXO, these results
demonstrate that neutron activation of natXe gas is a promising path to
producing a 127Xe calibration source.

##  3 Experimental demonstration

###  3.1 Stanford LXe TPC

The activated xenon was injected into the Stanford liquid xenon TPC to
demonstrate its use as a calibration source.

The TPC is housed in a cylindrical stainless steel chamber,
$20.3\text{\,}\mathrm{c}\mathrm{m}$ long by
$25.4\text{\,}\mathrm{c}\mathrm{m}$ diameter, maintained at
$165\text{\,}\mathrm{K}$ and filled with $27\text{\,}\mathrm{k}\mathrm{g}$ of
liquid xenon. The TPC itself, illustrated in Figure 5(a), consists of a
$13.5\text{\,}\mathrm{c}\mathrm{m}$ drift volume defined at the top by the
charge-sensing anode plane and at the bottom by a stainless steel cathode
grid. A uniform electric field is maintained by five field shaping rings
connected by $1\text{\,}\mathrm{G\SIUnitSymbolOhm}$ resistors. During the
measurements discussed here, the cathode was maintained at
$6\text{\,}\mathrm{k}\mathrm{V}$, producing an electric field of
$\sim$$400\text{\,}\mathrm{V}\mathrm{/}\mathrm{c}\mathrm{m}$ in the active
volume of the TPC. Scintillation light is detected by an array of 24 VUV-
sensitive FBK VUV-HD1 silicon photomultipliers (SiPMs) which are paired into
12 readout channels, located at the bottom of the chamber approximately
$4\text{\,}\mathrm{c}\mathrm{m}$ below the cathode. Ionization is detected by
a prototype nEXO charge tile, described in detail in Ref. [26]. The tile
consists of square gold pads deposited on a quartz substrate, connected into
strips of $9\text{\,}\mathrm{c}\mathrm{m}$ length and
$3\text{\,}\mathrm{m}\mathrm{m}$ pitch in both the $x$ and $y$ dimensions. The
strips are connected via feedthroughs to discrete preamplifiers (based on the
design in Ref. [27]), operating at $\sim$$165\text{\,}\mathrm{K}$ but located
outside the xenon space. Due to a limited number of feedthroughs in the xenon
chamber, some strips are ganged together into a single channel, as illustrated
in Figure 5(b). Both the light and the charge signals are digitized by Struck
SIS3316 digitizers at a rate of $62.5\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$.
Data acquisition is triggered by a two-fold coincidence requirement on the
SiPM channels, with the thresholds on each channel set at the mean pulse
height of single photoelectrons.

(a)

(b)

Figure 5: (a) CAD model of chamber and tile readout. The individual,
$3\text{\,}\mathrm{m}\mathrm{m}$-pitch strips on the charge-sensing anode are
ganged into channels as illustrated in (b).

###  3.2 Injection procedure

(a)

(b)

Figure 6: (a) Diagram of the source injection hardware. (b) Source pressure,
measured inside the 5.6 mL buffer volume during each fill/release cycle, and
detector trigger rate during the 127Xe injection test campaign. The shaded
grey region indicates a pause in data taking.

The activated xenon injection hardware is illustrated in Figure 6(a). The
cylinder filled with gas is connected through a valve (V1) to a tee, which
connects both to a high-pressure gauge and to the main xenon recirculation
loop through a second valve (V2). The 5.6 mL volume between the two valves
serves as a buffer volume which can be pressurized with activated gas from the
cylinder, then opened to the main xenon recirculation loop to inject the gas
into the system. During each fill/release cycle we measure the pressure in the
buffer volume to calculate the amount of gas injected.

The injection test was done in two stages, each of which consisted of $\sim$20
fill/release cycles. Xenon was continuously recirculated at approximately
$5\text{\,}\mathrm{S}\mathrm{L}\mathrm{P}\mathrm{M}$ to promote mixing. During
each injection, the SiPM trigger rate in the TPC was monitored to ensure that
the activity was reaching the detector. The results are shown in Figure 6(b).
A clear correlation between injected gas and trigger rate is observed,
indicating that at least some of the the activated gas mixed quickly into the
TPC. There is also evidence of increasing activity after the injection was
stopped, suggesting that the distribution of 127Xe throughout the TPC was not
immediate, and that mixing continued for some time.

###  3.3 Electron lifetime measurement

After injection, the xenon was recirculated for several days to ensure that
the 127Xe was distributed uniformly. Recirculation was then stopped and data
were taken in four separate acquisitions over the course of two days to
measure the electron lifetime in the TPC by mesuring the detected ionization
signal as a function of the drift time.

Charge collection signals on individual channels are included in our analysis
if they are above 3 times the RMS of the baseline noise. The charge drift time
is computed as the time difference between the 90% rise time of the waveform
and the scintillation trigger. Signals are first grouped into “clusters” if
their reconstructed drift times are within 3 $\mu$s of each other. We then
apply an event selection cut that requires each cluster to be reconstructed on
at least one $x$ and $y$ channel to ensure that events are fully reconstructed
in all three spatial dimensions. Finally, we select events for which the
charge-weighted average position in the $x$-$y$ plane falls in the central
$\pm$15 mm of the TPC, to select the region of the TPC with the finest-grained
position resolution. The total ionization energy of the event is reconstructed
by summing the charge from all the detected signals.

The data are then divided into 5 $\mu$s bins along the drift time axis, as
shown in Figure 7(a). An example bin is illustrated by Figure 7(b). The two
peaks expected from the 236 keV and 408 keV decay branches are clearly
observed. For each drift time bin, the charge energy spectrum is fitted with
two Gaussian distributions plus an exponentially-decaying background,

$f(x)=a_{1}\exp\left[-\frac{(x-\mu_{1})^{2}}{2\sigma_{1}^{2}}\right]+a_{2}\exp\left[-\frac{(x-\mu_{2})^{2}}{2\sigma_{2}^{2}}\right]+c\exp(-bx),$
(3.1)

where $a_{i}$, $\mu_{i}$ and $\sigma_{i}$ are the normalization constants, the
centroids and the standard deviations, respectively.

Because the charge sensors are bare strips, the charge induced in the sensors
by positive ions must be accounted for. Positive ions, which are produced
along with the electrons in an event, induce a surface charge on the strips of
opposite sign to the collected charge, leading to an effective reduction in
the detected charge relative to the number of electrons collected. On the
timescale of a single event, the ions are effectively stationary and their
induced charge can be calculated analytically as described in Ref. [26]. The
effective reduction of the detected charge is a function of the distance to
the readout strips and the number of strips over which the charge signal is
summed (the latter of which affects the surface area over which the ion-
induced charge is integrated). This effect is stronger for events closer to
the anode, i.e. at shorter drift times, so that its trend is qualitatively
opposite to that deriving from the finite electron lifetime. In this analysis,
we analytically calculate the strength of the induced charge for a single
strip and scale it by the average number of channels above threshold in each
event. For the higher-energy (lower-energy) peak, this number is 2.3 (2.0)
channels.

(a)

(b)

Figure 7:  Measurement of the electron lifetime in the Stanford TPC using
127Xe. (a) shows the distribution of all events in the charge vs. drift time
plane (blue points). The orange and green points show the fitted centroids of
the high- and low-energy peaks, respectively. The dashed line shows the best-
fit model using Equation 3.2. (b) shows the distribution for a single slice in
drift time. Two peaks corresponding to the 236 keV and 408 keV decay branches
are clearly observable. Figure 8: Fitted values of the electron lifetime for
both the low-energy (top) and high-energy (bottom) peaks. Data were taken in
four separate acquisitions; the $x$ error bars show the duration of the
individual acquisition, while the $y$ error bars show the uncertainty in the
fitted peak position. The dashed line and error band show the best-fit
electron lifetime and uncertainty from the combined dataset.

To extract the electron lifetime, we fit the centroid values for each peak to
the function:

$Q_{e}(t)=I(t)\times Q_{0}\,\,e^{(-t/\tau_{e})}$ (3.2)

where $I(t)$ is the analytically-calculated positive ion induction signal and
$Q_{0}$ is the initial ionized charge in the event. Both $Q_{0}$ and
$\tau_{e}$ are left as free parameters. Drift times beyond 50 $\mu$s are
subject to detector threshold effects and are not included in the fit. We fit
the data from each of the four acquisitions independently, then, having
confirmed their consistency, combine them into a single large dataset which is
fitted separately.

The best-fit values of $\tau_{e}$ are shown in Fig. 8, for each dataset
individually (data points) and for the full concatenated dataset (dashed
line). The concatenated dataset with the best-fit curve from Eq. 3.2 is shown
in Figure 7(a). We measure an electron lifetime of approximately
$45\text{\,}\mathrm{\SIUnitSymbolMicro s}$, which is likely limited by
outgassing from PVC-insulated wire used in the liquid xenon volume during this
particular run. We note that we obtain consistent results for both the 236 keV
and 408 keV peaks.

##  4 Projections of 127Xe calibrations in nEXO

To model 127Xe decay events in nEXO, we use the Geant4-based nexo-offline
simulation package, which has been used to estimate nEXO’s sensitivity to
$0\nu\beta\beta$ [2]. The software contains a detailed geometry of the nEXO
experiment and uses NEST [28], tuned to match the measurements in Ref. [6], to
model the production of ionization electrons and scintillation photons in the
liquid xenon medium. It then models the drift and diffusion of charge through
the TPC to the anode plane, properly accounting for the signal development and
noise in the charge readout electronics at the level of individual readout
channels [29]. In this work we start from the reconstruction algorithm used in
Ref. [2] and extend it with simplified models of the readout noise to generate
the charge and light detection signals.

###  4.1 Simulation of 127Xe calibration datasets

Figure 9: 127Xe decay events simulated in the active volume of the nEXO TPC.
We show the Monte-Carlo-truth distributions of ionization and scintillation
produced by the events, as well as their projections onto each observable
individually.

To generate calibration datasets, we simulate 127Xe decays distributed
uniformly throughout the liquid xenon volume in the nEXO TPC. The MC-truth
distribution of scintillation light versus charge is shown in Figure 9, before
noise and detection efficiencies are taken into account. As a result of
recombination fluctuations, the projections of these distributions into either
charge or light appear as two peaks. These two peaks will then be broadened by
the detection efficiencies and readout noise of the nEXO detector systems.

In the ionization channel, broadening is introduced both by fluctuations in
the charge detection efficiency (due to the finite electron lifetime) and by
noise in the readout electronics. The former is included directly in the
simulations; the MC-truth $z$ position is used to calculate an expected
attenuation probability from Equation 1.3, then the “detected” charge is drawn
from a binomial distribution to model fluctuations in electron attachment to
electronegative impurities in the liquid xenon. The readout noise is modeled
for each channel using a normal distribution with a $1\sigma$ width of 600
electrons, which is based on conservative estimates of the noise from
simulated waveforms analyzed with a trapezoidal filtering scheme. A similar
strategy was used in Ref. [2].

In the scintillation channel, broadening is introduced by fluctuations in the
light collection efficiency and the noise due to correlated avalanches in the
SiPMs. To model the former, the MC-truth lightmap $\epsilon_{\text{LM}}$ for
the nEXO TPC is calculated from a high-statistics light propagation simulation
using $5\times 10^{8}$ point-like photon sources distributed uniformly
throughout the TPC. This is performed using Chroma [30], which uses CUDA-
enabled GPUs for fast, high statistics photon transport simulations. The nEXO
detector geometry was imported directly into the simulation from CAD software.
Details of the optical properties included in the simulation can be found in
Ref. [2]. The resulting position-dependent efficiency values are binned in $r$
and $z$ with a bin size of 0.25 mm by 0.25 mm. A Gaussian blur with a width of
1 bin is then applied to the histogram to produce a smooth and continuous
lightmap.

The $xy$-positions of simulated 127Xe events are determined using an average
position of the charge signals for each energy deposition, weighted by the
detected charge. The $z$ position is calculated from the drift time and anode
position using an electron drift velocity of 0.171 cm/$\mu$s, appropriate for
the design field of 400 V/cm [1]. For each event, the number of detected
photons is given by

$\displaystyle n_{\mathrm{hit}}$
$\displaystyle=\mathcal{B}\left(S_{0},\epsilon_{QE}\cdot\epsilon_{LM}\left(x,y,z\right)\right)$
(4.1) $\displaystyle n_{\mathrm{av}}$
$\displaystyle=\mathcal{P}(n_{\mathrm{hit}}\cdot\Lambda)$ (4.2) $\displaystyle
n_{\mathrm{det}}$
$\displaystyle=\frac{n_{\mathrm{hit}}+n_{\mathrm{av}}}{\left(1+\Lambda\right)}$
(4.3)

where $\mathcal{B}$ and $\mathcal{P}$ represent Binomial and Poisson random
variables, respectively, $n_{\mathrm{hit}}$ is the number of scintillation
photons that create a signal in a SiPM, $\epsilon_{LM}(x,y,z)$ is the value of
the MC-truth lightmap at the reconstructed event position, $n_{\mathrm{av}}$
is the number of photons resulting from correlated avalanches, $\Lambda$ is
the correlated avalanche fraction, and $n_{\mathrm{det}}$ is the reconstructed
number of detected photons. In this study, we use $\epsilon_{QE}=0.186$ and
$\Lambda=0.2$, the current projections for nEXO [2]. Furthermore, we
conservatively assume that the $236\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$
peak from 127Xe will fall below the nEXO trigger threshold, and therefore
estimate calibration capabilities using only events from the
$408\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ peak.

###  4.2 Electron lifetime calibration

In a large detector like nEXO, the diffusion of an electron cloud as it drifts
across the TPC spreads the charge across more channels, leading to more
channels collecting charge for events with longer drift times. This in turn
causes more charge to be "lost" below threshold, which results in a systematic
bias in the measured charge that mimics a shorter electron lifetime. To avoid
this issue, we adopt a "no-threshold" analysis on the charge signals via the
following procedure.

First, a temporary 1200 e- threshold is applied to find channels that collect
large amounts of charge. Next, the charge-weighted position of the event is
computed using these channels. Finally, the charge is reconstructed by adding
the integrated signal on the five nearest channels in each direction (20 in
total, counting $\pm\,x$ and $\pm\,y$), regardless of the amount of detected
charge on each. As the events tend to be highly localized, this enables us to
include the charge collected on channels that do not pass the initial
threshold cut, at the cost of introducing additional noise into the
reconstructed charge. In addition, the effect of induced charge from the
positive ions introduces a bias for events at short drift times similar to
that observed in Section 3. This effect is exacerbated by the summing of many
channels, which increases the induction component to the signal, and leads to
a strong rolloff in the reconstructed charge for events closer to the anode.
The effect is illustrated by the blue points in Figure 10(a). To correct for
this effect, we again use an analytical model of induced charge to calculate
the expected bias, making the simplifying assumption that all of the charge is
located at the reconstructed charge-weighted position. After applying the
correction and performing the "no-threshold" analysis, the biases from
diffusion and positive ion induction are removed and the data are well-
described by an exponential attenuation, as illustrated by the red points in
Figure 10(a).

We evaluate the precision of electron lifetime calibrations as a function of
both $\tau_{e}$ and $N$, where $N$ is the number of events in a given
calibration dataset. For each value of $\tau_{e}$, $10^{7}$ events are
simulated in the active volume of the TPC. Smaller datasets are created by
randomly sampling subsets of these simulations, with replacement. The electron
lifetimes considered here range between $\tau_{e}=$1 – 10 ms, and the sizes of
datasets range from $N=10^{3}$ – $10^{6}$ events. For each ($\tau_{e}$,$N$)
pair, we analyze seven sample datasets; we run the fits and extract $\lambda$
(where $\lambda=1/\tau_{e}$) and $\sigma_{\lambda}$ for each, then take the
mean and standard deviation of the seven values to estimate the average
performance and the size of expected statistical fluctuations. The result for
datasets with $\tau_{e}=10$ ms, as a function of $N$, is shown in Figure
11(a). The average estimated uncertainty, shown by the solid red line, is
found to be proportional to $1/\sqrt{N}$, indicating that statistical
fluctuations are the dominant source of error in our reconstruction of
$\tau_{e}$. For $\tau_{e}$ significantly greater than the average drift time
(i.e., $\tau_{e}\geq\sim$1 ms) the absolute uncertainty in $\lambda$ is
approximately independent of the true value of $\tau_{e}$, meaning these
conclusions hold for any $\tau_{e}$ considered here.

(a) Figure 10: (a) Reconstructed charge peak position as a function of the
drift time, for a simulated dataset with 30,000 127Xe decays in the active
volume of the nEXO TPC. The uncorrected data (blue) show the positive ion
induction effect on the charge measurement. The data are corrected (red) using
the estimated event position as described in the text. The corrected data are
fitted to an exponential to determine the electron lifetime.

(a)

(b)

Figure 11: (a) Uncertainty in the calibration of the inverse electron lifetime
$\lambda$ as a function of the number of simulated events in the active
volume. The central values and error bands represent the mean and standard
deviation of the fitted $\lambda$s over seven simulated datasets. An
uncertainty of $10^{-2}$ ms-1 corresponds to a relative uncertainty of 10%.
(b) The expected energy resolution of nEXO at $Q_{\beta\beta}=2.457$ MeV as a
function of the integration time, assuming a 1 Bq steady-state activity in the
liquid xenon. The expected energy resolution assuming $\tau=10$ ms and no
uncertainty on the electron lifetime is shown by the dashed red line.

Using these results, we then calculate the impact of the electron lifetime
calibration uncertainty on the total energy resolution of nEXO. The simulated
number of events in the TPC is converted into an integration time, assuming a
127Xe activity of $1\text{\,}\mathrm{B}\mathrm{q}$ distributed throughout the
entire xenon system. The results are shown in Figure 11(b). When fewer
electrons are detected, charge noise makes up a larger fraction of the total
charge signal. This results in an asymptotically worse energy resolution at
shorter electron lifetimes, despite the independence of the absolute
uncertainty in $\lambda$ on the value of $\tau_{e}$. We find that, for the
benchmark case where $\tau_{e}=10$ ms, the uncertainty in the electron
lifetime calibration after a 24 hr integration period introduces only a 0.03%
absolute broadening of the energy resolution, indicating that a
$1\text{\,}\mathrm{B}\mathrm{q}$ source is able to calibrate the electron
lifetime on a daily basis with sufficient precision for nEXO.

###  4.3 Lightmap calibration

We evaluate the lightmap reconstruction capability by determining the degree
to which a lightmap reconstructed from a set of simulated 127Xe calibration
data deviates from the MC-truth lightmap. Figure 12(a) shows the uncorrected
calibration data, with variations in the photon transport efficiency smearing
the light distribution significantly. The photon transport efficiency
associated with each calibration event is given by

$\epsilon_{C}=\frac{n_{\mathrm{det}}}{\left<S\right>\times\epsilon_{QE}}$
(4.4)

where $\left<S\right>=18114$ is the expected number of photons produced for an
event corresponding to the 408 keV peak, as determined by NEST. Figure 12(b)
shows the MC-truth scintillation peaks for a sample calibration dataset, along
with the scintillation peaks reconstructed both with and without a lightmap
correction. From this figure, it should be noted that the width of the
uncorrected peaks is dominated by variations in the photon transport
efficiency across the detector, while statistical variations in the production
and detection of photons contribute less significantly.

(a)

(b)

Figure 12: Monte Carlo simulated data from a 127Xe calibration campaign. In
(a), the two peaks are clearly separable with a linear cut in the space of
detected electrons and photons. Both peaks appear smeared toward low numbers
of scintillation photons as a result of events originating in regions of lower
photon transport efficiency. In (b), the effect of a lightmap correction on
the scintillation peak width is shown. The green histogram shows the true
number of scintillation photons produced, the purple shows how this can be
reconstructed assuming the lightmap is calibrated perfectly (i.e. corrected by
the same lightmap from which the detected photons were sampled), and the red
shows the uncorrected data scaled by a constant.

There are many techniques that can be used to extract a lightmap from
calibration data. The simplest – choosing an appropriate bin size and binning
the data in 3 dimensions – is not well-suited to nEXO due to the large
detector volume. To properly capture the spatial variation in the lightmap,
the binning must be sufficiently fine; however this requires a large number of
events to avoid empty bins and significant statistical fluctuations between
adjacent bins. Techniques that avoid this drawback by interpolating or
smoothing over fewer events require the choice of a length parameter based on
the data, while the optimal choice of such a parameter often varies by region
throughout the TPC. A neural network model avoids these problems.
Specifically, it requires no choice of bin size or length scale, it is
adaptive to regions in the TPC where the efficiency varies at differing rates,
and the computation time scales well with increasing dataset sizes.

The neural net can be trained on the 3-dimensional position coordinates and
associated efficiencies of a calibration dataset, giving a continuous
efficiency function defined across the full TPC volume. We use a basic feed-
forward neural network architecture with five hidden layers. Although the nEXO
TPC is designed to be approximately cylindrically symmetric, we reconstruct
the lightmap in all three spatial dimensions in order to allow for the
possibility of unexpected asymmetries in the real experiment. Neural net
hyperparameters were selected to minimize the error in the reconstructed
lightmap. Once a satisfactory neural net architecture was chosen, it was
trained repeatedly on calibration datasets of varying sizes. The loss function
evaluated on both the training and validation datasets was monitored to ensure
the neural net was not overfitting to training data. Figure 13 shows the full
sequence of simulating calibration data from the MC-truth lightmap and then
reconstructing the lightmap using the neural net.

Figure 13: Left: the MC-truth lightmap from which the photon transport
efficiency is sampled during the simulation of calibration data. Middle:
Points representing the reconstructed positions of $10^{4}$ calibration
events, colored by the photon transport efficiency as calculated from the
number of detected photons and the expected number predicted by NEST. Right:
the reconstructed lightmap resulting from training a neural net on the
calibration dataset.

Subsets of the simulated events of varying sizes were fed to the neural net to
understand the dependence of the reconstructed lightmap accuracy on the number
of 127Xe decays in the active volume. We look specifically at four dataset
sizes ranging from $10^{3}$ to $10^{6}$ events. For each of these four dataset
sizes, 25 calibration datasets were sampled, with replacement, from the full
set of $10^{7}$ events.

Figure 14: Reconstructed lightmap accuracy at the chosen angular slice of the
3-dimensional lightmap for four sample datasets ranging from $10^{3}$ to
$10^{6}$ events. The red dotted line shows the boundary that defines the inner
2 tonnes, while the grey band around the detector edges indicates the regions
removed by the 20 mm standoff cut.

To compare the reconstructed lightmap to the MC-truth lightmap, the trained
neural net is passed a uniformly-spaced grid of points along an arbitrary
angular slice of the TPC. Dividing point-by-point by the corresponding points
in the MC-truth lightmap gives the spatially-dependent lightmap reconstruction
accuracy in $r$ and $z$ for the chosen slice. This is plotted in Fig. 14 for a
representative dataset for each of the four dataset sizes. The lightmap error
is given by the standard deviation of this accuracy distribution.

In the absence of a sufficient density of calibration events to capture the
spatial dependence of the true lightmap, the neural net converges to a largely
uniform lightmap centered around the average photon transport efficiency. The
effect of this overly-uniform lightmap is to introduce a systematic bias in
regions with large gradients in the lightmap. With larger calibration
datasets, the reconstructed lightmaps converge toward the true lightmap, with
the residual errors receding to the edges of the TPC. With increasing
statistics, the reduction in residual error continues to the point at which
further improvements are limited by both the distance between $\gamma$-ray
energy depositions in a single event and the accuracy of the position
reconstruction. When the average spacing between calibration events approaches
the uncertainty in the event position, further improvements cannot be achieved
with larger datasets. Using this lightmap calibration technique, it will be
possible to achieve a lightmap error of less than 1% in the full volume after
calibrating for 16 days with a 127Xe activity maintained at 1 Bq. In the inner
two tonnes, this error will fall below 0.5%. These results are summarized in
Figure 15.

Figure 15: The lightmap error achievable when reconstructing the lightmap from
a calibration dataset with a give number of decays in the TPC volume. The
solid lines and shaded regions represent the mean and standard deviation of
the lightmap errors calculated from all 25 datasets at each point.

##  5 Conclusions and prospects for nEXO

In this work, we have demonstrated the viability of using 127Xe produced via
the neutron activation of natXe gas as an internal calibration source for
liquid xenon TPCs. Such a source was procured, assayed for radioactive
contaminants, and used to measure the electron lifetime in a liquid xenon TPC
using prototype instrumentation for the nEXO experiment. Simulations show that
this source can be used to calibrate both the electron lifetime and the
lightmap in nEXO without requiring detector downtime for dedicated calibration
campaigns.

In this work we have assumed an activity of 1 Bq is maintained in the liquid
xenon continuously. In practice, this could be accomplished by an initial
injection of 1 Bq followed every two weeks by injections of $\sim$0.25 Bq to
maintain a near-constant activity. Figure 16 shows the amount of xenon gas
corresponding to an activity of 1 Bq for two different scenarios: a source
with the same activity as the one used in this work, and a source with an
activity that is larger by a factor of four. In the latter case, we find that,
even up to one year after activation, the source can supply 0.25 Bq injections
with less than 1 g of activated xenon gas, resulting in negligible dilution of
the enriched xenon in nEXO. Metering the injected gas could be accomplished by
a simple volume-sharing scheme such as that used in Section 3, but with a
larger ratio of initial volume to expansion volume to permit higher precision.

Figure 16: The inverse specific activity for the source used in Section 3
(blue) and a source with a factor of four higher initial activity (green). In
the latter case, a 1 Bq activity corresponds to less than 1 g of gas up to 300
days after activation.

## Acknowledgments

The authors thank Harrison Redman of the Health Physics Group at Stanford for
assistance with HPGe radioassay measurements, and Wesley Frey at MNRC for help
with source development. The authors gratefully acknowledge support for nEXO
from the Office of Nuclear Physics within DOE’s Office of Science, and NSF in
the United States; from the National Research Foundataion (NRF) of South
Africa; from NSERC, CFI, FRQNT, NRC, and the McDonald Institute (CFREF) in
Canada; from IBS in Korea; from RFBR in Russia; and from CAS and NSFC in
China.

## References

  * [1] nEXO Collaboration, S. A. Kharusi, A. Alamre, J. B. Albert, M. Alfaris, G. Anton et al., _nEXO Pre-Conceptual Design Report_ , 1805.11142.
  * [2] nEXO Collaboration, G. Adhikari, S. A. Kharusi, E. Angelico, G. Anton, I. J. Arnquist et al., _nEXO: neutrinoless double beta decay search beyond $10^{28}$ year half-life sensitivity_, _Journal of Physics G: Nuclear and Particle Physics_ 49 (2021) 015104, [2106.16243].
  * [3] EXO-200 collaboration, E. Conti et al., _Correlated fluctuations between luminescence and ionization in liquid xenon_ , _Phys. Rev. B_ 68 (2003) 054201, [hep-ex/0303008].
  * [4] E. Aprile, K. L. Giboni, P. Majewski, K. Ni and M. Yamashita, _Observation of Anti-correlation between Scintillation and Ionization for MeV Gamma-Rays in Liquid Xenon_ , _Phys. Rev. B_ 76 (2007) 014115, [0704.1118].
  * [5] XENON collaboration, E. Aprile et al., _Energy resolution and linearity of XENON1T in the MeV energy range_ , _Eur. Phys. J. C_ 80 (2020) 785, [2003.03825].
  * [6] EXO-200 collaboration, G. Anton et al., _Measurement of the scintillation and ionization response of liquid xenon at MeV energies in the EXO-200 experiment_ , _Phys. Rev. C_ 101 (2020) 065501, [1908.04128].
  * [7] EXO-200 collaboration, C. G. Davis et al., _An Optimal Energy Estimator to Reduce Correlated Noise for the EXO-200 Light Readout_ , _JINST_ 11 (2016) P07015, [1605.06552].
  * [8] N. Ackerman et al., _The EXO-200 detector, part II: Auxiliary Systems_ , 2107.06007.
  * [9] A. Manalaysay et al., _Spatially uniform calibration of a liquid xenon detector at low energies using 83m-Kr_ , _Rev. Sci. Instrum._ 81 (2010) 073303, [0908.0616].
  * [10] L. W. Kastens, S. B. Cahn, A. Manzur and D. N. McKinsey, _Calibration of a Liquid Xenon Detector with Kr-83m_ , _Phys. Rev. C_ 80 (2009) 045809, [0905.1766].
  * [11] LUX collaboration, D. S. Akerib et al., _${}^{83\textrm{m}}$ Kr calibration of the 2013 LUX dark matter search_, _Phys. Rev. D_ 96 (2017) 112009, [1708.02566].
  * [12] LUX collaboration, D. S. Akerib, H. M. Araújo, X. Bai, A. J. Bailey, J. Balajthy, P. Beltrame et al., _Tritium calibration of the LUX dark matter experiment_ , _Phys. Rev. D_ 93 (Apr, 2016) 072009.
  * [13] E. M. Boulton et al., _Calibration of a two-phase xenon time projection chamber with a 37Ar source_, _JINST_ 12 (2017) P08004, [1705.08958].
  * [14] K. Ni, R. Hasty, T. M. Wongjirad, L. Kastens, A. Manzur and D. N. McKinsey, _Preparation of Neutron-activated Xenon for Liquid Xenon Detector Calibration_ , _Nucl. Instrum. Meth. A_ 582 (2007) 569–574, [0708.1976].
  * [15] R. Lang, A. Brown, E. Brown, M. Cervantes, S. Macmullin, D. Masson et al., _A 220Rn source for the calibration of low-background experiments_, _Journal of Instrumentation_ 11 (2016) P04004–P04004.
  * [16] XENON collaboration, E. Aprile et al., _Results from a Calibration of XENON100 Using a Source of Dissolved Radon-220_ , _Phys. Rev. D_ 95 (2017) 072008, [1611.03585].
  * [17] D. S. Akerib, S. Alsum, H. M. Araújo, X. Bai, A. J. Bailey, J. Balajthy et al., _Ultralow energy calibration of LUX detector using ${}^{127}\mathrm{Xe}$ electron capture_, _Phys. Rev. D_ 96 (Dec, 2017) 112011.
  * [18] D. J. Temples, J. McLaughlin, J. Bargemann, D. Baxter, A. Cottle, C. E. Dahl et al., _Measurement of charge and light yields for ${}^{127}\mathrm{Xe}$ $l$-shell electron captures in liquid xenon_, _Phys. Rev. D_ 104 (Dec, 2021) 112001.
  * [19] M.-M. Bé, V. Chisté, C. Dulieu, M. Kellett, X. Mougeot, A. Arinc et al., _Table of Radionuclides_ , vol. 8 of _Monographie BIPM-5_.  Bureau International des Poids et Mesures, Pavillon de Breteuil, F-92310 Sèvres, France, 2016.
  * [20] McClellan Nuclear Research Center, _MNRC User’s Data - CIF PTS NTD Spectra_ , 2015.
  * [21] M. Chadwick, M. Herman, P. Obložinský, M. Dunn, Y. Danon, A. Kahler et al., _ENDF/B-VII.1 Nuclear Data for Science and Technology: Cross Sections, Covariances, Fission Product Yields and Decay Data_ , _Nuclear Data Sheets_ 112 (2011) 2887 – 2996.
  * [22] A. Koning, D. Rochman, J.-C. Sublet, N. Dzysiuk, M. Fleming and S. van der Marck, _TENDL: Complete Nuclear Data Library for Innovative Nuclear Science and Technology_ , _Nuclear Data Sheets_ 155 (2019) 1–55.
  * [23] R. H. M. Tsang, A. Piepke, D. J. Auty, B. Cleveland, S. Delaquis, T. Didberidze et al., _GEANT4 models of HPGe detectors for radioassay_ , _Nucl. Instrum. Meth. A_ 935 (2019) 75–82, [1902.06847].
  * [24] B. N. L. National Nuclear Data Center, “Nudat-2.8 (nuclear structure and decay data).”
  * [25] E. Kondaiah, N. Ranakumar and R. Fink, _Thermal neutron activation cross sections for kr and xe isotopes_ , _Nuclear Physics A_ 120 (1968) 329–336.
  * [26] nEXO Collaboration, M. Jewell et al., _Characterization of an Ionization Readout Tile for nEXO_ , _JINST_ 13 (2018) P01006, [1710.05109].
  * [27] L. Fabris, N. W. Madden and H. Yaver, _A fast, compact solution for low noise charge preamplifiers_ , _Nucl. Instr. Meth._ A424 (Mar, 1999) 545–551.
  * [28] M. Szydagis, J. Balajthy, J. Brodsky, J. Cutter, J. Huang, E. Kozlova et al., _Noble Element Simulation Technique v2.0_ , July, 2018.  10.5281/zenodo.1314669.
  * [29] nEXO Collaboration, Z. Li et al., _Simulation of charge readout with segmented tiles in nEXO_ , _JINST_ 14 (2019) P09020, [1907.07512].
  * [30] B. Land, S. Seibert and A. LaTorre, “Chroma.” https://github.com/BenLand100/chroma, 2021.