e1e-mail<EMAIL_ADDRESS>
11institutetext: Please look at the end of the paper for the author list

# Measurement of isotopic separation of argon with the prototype of the
cryogenic distillation plant Aria for dark matter searches

The DarkSide-20k Collaboration${}^{\text{a,1}}$

###### Abstract

The Aria cryogenic distillation plant, located in Sardinia, Italy, is a key
component of the DarkSide-20k experimental program for WIMP dark matter
searches at the INFN Laboratori Nazionali del Gran Sasso, Italy. Aria is
designed to purify the argon, extracted from underground wells in Colorado,
USA, and used as the DarkSide-20k target material, to detector-grade quality.

In this paper, we report the first measurement of argon isotopic separation by
distillation with the 26 m tall Aria prototype. We discuss the measurement of
the operating parameters of the column and the observation of the simultaneous
separation of the three stable argon isotopes: ^36Ar, ^38Ar, and ^40Ar. We
also provide a detailed comparison of the experimental results with commercial
process simulation software.

This measurement of isotopic separation of argon is a significant achievement
for the project, building on the success of the initial demonstration of
isotopic separation of nitrogen using the same equipment in 2019.

††journal: Eur. Phys. J. C

##  1 Introduction

Aria is an industrial-scale plant comprising a $350$ m cryogenic distillation
column for isotopic separation Darkside:2021 , the tallest ever built, and is
currently being installed in a mine shaft at Carbosulcis S.p.A., Nuraxi-Figus
(SU), Italy. Aria is designed to reduce the isotopic abundance of ^39Ar in the
low-radioactivity argon, or undergorund argon (UAr), extracted from
underground sources and used for dark matter searches.

A full description of the Aria plant and of the column structure and
performance can be found in Ref. Darkside:2021 .

In recent years, interest in UAr has grown to include its use in searches for
neutrino-less double-beta decay with the LEGEND experiment Legend ; in the
measurement of coherent elastic neutrino-nucleus scattering with the COHERENT
experiment Coherent ); and in measuring low-energy processes with one of the
planned modules of the DUNE experiment (PNNL-SA-171088). The projected demand
for UAr varies greatly, from one tonne for the COHERENT experiment to several
thousand tonnes for a detector of similar size to DUNE. The required ^39Ar
depletion factor with respect to atmospheric argon levels spans from a factor
$1000$ for COHERENT and DarkSide-20k to more than $10000$ for future dark
matter searches with, e.g., a tonne-scale DarkSide-LowMass experiment
(arXiv:2209.01177).

The first application of the Aria plant will be to purify 120 t of UAr through
cryogenic distillation, which will be used in the DarkSide-20k experiment
Aalseth:2018gq . However, the high depletion factors required for future
projects also demonstrate the need for isotopic separation using Aria.

The initial phase of the Aria program includes assessing the efficiency of the
distillation process and determining the thermodynamic parameters of the
prototype distillation plant. A successful nitrogen isotope separation
campaign was conducted with the same plant in 2019 and described in Ref.
Darkside:2021 . The prototype plant is a shortened version of the Aria column,
consisting of a top condenser, a bottom reboiler, and one central module (out
of $28$). It is equipped with all the auxiliary instrumentation of the full
column and located in a surface building at Carbosulcis S.p.A. The results
obtained using nitrogen were very promising. However, argon has distinct
thermodynamic properties, such as its boiling point in relation to pressure
and its latent heat at the boiling point, compared to nitrogen. Therefore, a
separate characterization of the plant’s performance is necessary when using
argon.

This article presents the results of a distillation run for the separation of
argon isotopes conducted in 2021. Whereas ^39Ar is only present in traces, the
two stable isotopes other than ^40Ar in atmospheric argon, ^36Ar and ^38Ar,
have a non-negligible isotopic abundance of 0.334% and 0.063%, respectively
Lee ; Bohlke , and can be used to characterize the distillation performance of
the plant.

With respect to Ref. Darkside:2021 , this paper presents a more in-depth
characterization of the plant performance, including hydraulic parameters,
details of the distillation process, and a thorough comparison between the
actual performance and that obtained using a commercial process simulation
software. We have also demonstrated through experimentation that multi-
component distillation modeling is not required for argon distillation,
despite the presence of more than two isotopes in the input feed. This was
already hypothesized in Ref. Darkside:2021  without experimental proof.

Here we briefly reproduce a few basic formulae of distillation columns. For a
two-component distillation at total reflux, i.e. without a significant amount
of argon entering or exiting the column during distillation, the separation
between isotopes $i$ and $j$ between the top (T) and the bottom (B) of the
column is given by:

$S^{TB}_{i-j}=(\alpha_{i-j})^{N},$ (1)

where $\alpha_{i-j}$ is the relative volatility between the lighter, more
volatile isotope $i$ and the heavier, less volatile isotope $j$, and $N$ is
the number of theoretical stages. The separation can be measured from:

$S^{TB}_{i-j}=\frac{(R_{i-j})_{T}}{(R_{i-j})_{B}},$ (2)

where $(R_{i-j})_{T(B)}$ is the isotopic ratio, i.e. the relative abundance
between the lighter and the heavier isotope, measured at the top (bottom) of
the column.

The logarithm of the relative volatility, also known as the Vapor Pressure
Isotope Effect (VPIE), is given by:

$\mathrm{VPIE}=\ln{\alpha_{i-j}}\simeq\ln\frac{{P_{i}}}{{P_{j}}}\simeq\ln\frac{{f_{c}}}{{f_{g}}},$
(3)

with $P_{i}$ and $P_{j}$ the vapor pressures of the two isotopes $i$ and $j$
and $f_{c}$ and $f_{g}$ the reduced partition functions in the condensed phase
and in the ideal gas.

In a column with structured packing like Aria, the distillation performance is
related to the height equivalent to a theoretical plate or stage, HETP, and is
given by:

$\mathrm{HETP}=\frac{h}{N}$ (4)

where $h$ is the cumulative height of the packing material of the column.
Combining Eq. 1 and Eq. 4, it follows that:

$\mathrm{HETP}={\frac{h\times\ln(\alpha_{i-j})}{\ln(S^{TB}_{i-j})}}$ (5)

The rest of the paper is organized as follows. In Sect. 2 we discuss the
simulation programs used in this work. In Sect. 3 we present the experimental
setup. In Sect. 4 we discuss the pressure drop along the column and the liquid
holdup. In Sect. 5 and 6 we discuss the distillation measurements and in Sect.
7 we compare them with the simulations. In Sect. 8 we conclude and give some
perspectives. In A we derive an estimate of the relative volatilities of ^40Ar
with respect to ^36Ar and ^38Ar from existing measurements.

##  2 The simulation programs

We use the Aspen HYSYS (©2022 Aspen Technology Inc) software, a powerful tool
for industrial process simulations, to model our system. HYSYS performs
rigorous distillation calculations, i.e. solve the equations for mass,
component, energy, and equilibrium balance for each stage of the column.
Rigorous calculations also allow the simulation of multi-component
distillation processes, with three or more components modeled simultaneously.
The standard HYSYS library of chemicals does not include isotopes. For this
simulation, we added the argon isotopes to the HYSYS library, together with
their molecular weight and Antoine equation parameters.

The results for binary distillation obtained with HYSYS were compared with
those from a calculation with the graphical McCabe-Thiele (MCT) method.

We also used the Sulcol^T^M $3.5$ (Sulzer) software package for structured and
random packing design to model the hydraulic parameters of the distillation
column. Taking input properties, such as the density and viscosity of the gas
and liquid phases, the surface tension, and the mass flow rates of the gas and
the liquid inside the column, Sulcol calculates column parameters including
the pressure drop per unit length and the liquid hold-up of the column.

##  3 The experimental setup and measurements

Before starting operations, the column and all the process lines were
evacuated through a scroll pump. The column was then filled with argon gas.
The total amount of argon filling the column was about $250$ kg.

Plant operation followed a procedure similar to our previous run with nitrogen
Darkside:2021 . Automatic process control was implemented for most of the
system. For this run, it took $12$ days of operation to stabilize the plant
and start the distillation measurements reported below, which lasted about 4
days. As for the past run, nitrogen was used as the refrigerant fluid.

The distillation was performed in total reflux mode. About 1 L h-1 of gas was
extracted for the sampling, corresponding to an overall loss during the run of
$<$0.1% of the loaded argon mass. Therefore, we expect Eq. 1 to hold.

The operating conditions of the argon distillation runs are presented in Fig.
1. The top panel (Fig. 1) displays the vapor mass flow rate of nitrogen in the
auxiliary system (red); the bottom panel (Fig. 1) shows the average pressure
$p$ inside the column (blue) as a function of time. The time averages of the
same quantities are summarized in Table 1, together with the maximum pressure
variation during the run, $\delta_{p}$. The collected data is divided in three
runs (A, B, and C) corresponding to changing operating conditions. At the
beginning of run A, the vapor mass flow rate of nitrogen was about 10% lower
than the plant design value of 550 kg h-1. In run B, we kept the vapor mass
flow rate unchanged and lowered the pressure inside the column by increasing
the amount of nitrogen introduced into the auxiliary system. In run C we
increased the vapor mass flow rate to a value about 10% higher than the design
value. The lowest average pressure achieved inside the column was about 18%
larger than the design value described in Ref. Darkside:2021 .

Figure 1: Vapor mass flow rate of nitrogen in the auxiliary system (red, 1)
and average pressure inside the column $p$ (blue, 1) as a function of time.
The shaded regions correspond to periods with changing conditions in the
column and are not included in the data analysis.

With argon distillation, we have achieved a ten-fold lower pressure variation
during the run, than during the nitrogen distillation run discussed in Ref.
Darkside:2021 .

Table 1: Time averaged operational parameters during runs A, B and C: pressure inside the column, maximum pressure variation $\delta_{p}$ during the run, and vapor mass flow rate of nitrogen in the auxiliary system. |  $p$ (bar) |  $\delta_{p}$ (bar) | vapor flow rate (kg h-1)
---|---|---|---
A | $1.824\pm 0.001$ | $\pm 0.02$ | $491\pm 1$
B | $1.544\pm 0.001$ | $\pm 0.02$ | $487\pm 1$
C | $1.953\pm 0.001$ | $\pm 0.04$ | $624\pm 1$

Table 2 shows three derived operational parameters: the time-averaged
saturation temperature, $T$, obtained from the Antoine equation using the
measured pressures from Table 1; the argon vapor mass flow rate, V, derived
from the energy balance in the reboiler; and the thermal power of the process,
given by the product of nitrogen latent heat and nitrogen mass flow rate.

Table 2: Derived operational parameters: time-averaged saturation temperature, $T$, argon vapor mass flow rate, $V$, and the thermal power of the process, $Q$. |  $T$ (K) |  V (kg h-1) |  $Q$ (kW)
---|---|---|---
A |  $93.3\pm 0.1$ |  $543\pm 1$ |  $23.0\pm 0.1$
B |  $91.5\pm 0.1$ |  $544\pm 1$ |  $23.8\pm 0.1$
C |  $94.0\pm 0.1$ |  $689\pm 1$ |  $29.7\pm 0.1$

###  3.1 The gas sampling system

The instrument used to measure the isotopic distillation performance is an MKS
Instruments, Cirrus™ 3-XD quadrupole mass spectrometer. During the run, the
mass spectrometer was continuously in operation to sample the gas coming from
the top (Top) and bottom (Bottom) of the column, and from the gas bottles
(Feed), as shown in Fig. 2.

Figure 2: Schematic of the gas sampling system.

The gas from the column flows through 1/8" copper lines between 4 m (Bottom)
and 30 m (Top) in length. These lines are split before the mass spectrometer.
Most of the gas flow, controlled via mass flow controllers set around $50$ mL
min-1, purge the lines and is vented to air. The gas along the other branch
flows through, respectively, a 1/16" ($1$ mm ID) stainless steel (SS), $7$ cm
long, capillary line; a multi-port valve (VALCO Instruments); another 1/16"
($1$ mm ID), $15$ cm long SS capillary line; and a manual pressure regulator,
to limit the gas pressure to a maximum of $1$ bar, thus protecting the input
of the mass spectrometer. The gas then passes through a 1/16" ($1$ mm ID), $1$
m long SS line, and a $19$ cm ($0.07$ mm ID) long SS capillary, to the mass
spectrometer. The inlet flow of the mass spectrometer is $20$ mL min-1. Before
the run, the response time of the sampling system was determined by connecting
a bottle of argon and one of nitrogen to the multi-inlet valve. It took about
$10$ min to see a change in the measured gas composition and 40-45 min to see
it stabilize. This lag time is caused by the combination of the relatively
large buffer volume of the manual pressure regulator and the low gas inlet
flow. During gas sampling of argon from the column, the inlet valve was
programmed to switch among the three input capillaries every hour. This is a
compromise between the response time of the sampling system and the need for
quickly detecting changes in isotopic mass fraction in the column when the
operating parameters are changed. In Fig. 3 we show an example of a background
subtracted mass spectrum from the spectrometer during the run, zoomed in the
region of interest for argon isotopes

Figure 3: Example of a background subtracted mass spectrum measured during the
run, zoomed in the region of interest for argon isotopes.

and featuring well-separated peaks.

We obtain the isotopic ratios $R_{36-40}$ and $R_{38-40}$ as ratios of the
^36Ar (^38Ar) and ^40Ar peak heights. $R_{36-40}$ is displayed vs time in Fig.
4, where each point corresponds to a measurement derived from a spectrum like
the one in Fig. 3.

Figure 4: Isotopic ratio, $R_{36-40}$, as a function of time, for a time slice
of run C. Every hour, the inlet valve is switched from Top (red) to Bottom
(blue) and Feed (green). The isotopic ratios were calculated at the time when
the inlet valve is switched to the Bottom sampling (e.g. at times marked by
dashed vertical lines), and only considering the last 15 min before each
switching of the inlet valve.

Fig. 4 clearly shows that the measurement of the isotopic ratio settles to a
plateau in approximately 45 min. Because of this delayed response, we only
considered the last 15 min of each 1-hour sampling period during the
measurement of the distillation performance,. The isotopic ratios were
calculated with the inlet valve switched to the Bottom sampling (e.g. at times
marked by dashed vertical lines). We estimated a systematic uncertainty of the
order of 2% on the measured isotopic ratio due to this choice, by comparing it
with other possible choices of times for calculating the isotopic ratios.

##  4 Pressure drop along the column and liquid holdup

The Aria column is filled with structured stainless steel packing (Sulzer CY
gauze) sections $2.56\text{\,}\mathrm{m}$ long, four per module, interleaved
with a liquid distributor to optimize flow uniformity across the column cross-
section, whose structure is discussed below. The purpose of the packing is to
support the liquid phase, for optimal thermodynamic contact between the rising
vapor and the descending liquid. During operations, the liquid is dispersed as
a film coating the surfaces of the packing and immersed in the vapor phase.
This configuration is referred to as irrigated bed.

Two very relevant operational parameters of the packed distillation column are
the pressure drop per unit length and the liquid holdup Kister:1992 ,
described hereafter.

###  4.1 Pressure drop

The pressure drop is due to friction and directional changes in the gas flow
inside the column caused by the packing and the distributors. Therefore, the
total pressure drop in the column, $\Delta p_{tot}$, is given by the sum of
the pressure drop in the distributors, $\Delta p_{\mathrm{distributors}}$, and
in the wetted packing elements, $\Delta p_{\mathrm{packing,irr}}$ as:

$\Delta p_{tot}=\Delta p_{\mathrm{distributors}}+\Delta
p_{\mathrm{packing,irr}}$ (6)

####  4.1.1 Calculation of the pressure drop due to the distributors

The distributors, whose picture can be found in Ref. Darkside:2021 , are made
by horizontal plates, intercepting the downward liquid argon flow along the
column. The plates are perforated, with 10 cm hollow vertical pipes or
chimneys with a top cap at the hole locations on the upper side, which are
uniformly distributed on the plate surface. The vapor flows up through these
chimneys. The liquid formed on the distributor plate is streamed, through 0.3
cm holes located at 3 cm, 4 cm, and 5 cm height in the chimneys, to the
packing section below.

In our calculations, the distributed pressure drop due to the short pipes can
be neglected and the distributor is, therefore, approximated as a perforated
sheet with holes. The pressure drop from the four distributors, $\Delta
p_{\mathrm{distributors}}$, is determined following Ref. caleffi28 :

$\Delta
p_{\mathrm{distributors}}=4\times\zeta\;\rho{{}_{\mathrm{G}}}\;\frac{v_{\mathrm{G}}^{2}}{2},$
(7)

where $\zeta$ is the concentrated loss coefficient, and
$\rho{{}_{\mathrm{G}}}$ and $v_{\mathrm{G}}$ are the density and the velocity
of the argon gas, respectively. The value of $\zeta$ is evaluated by
interpolating the data reported on page 18 of Ref. caleffi28  for $\zeta$ vs.
$A^{*}$/$A$, and using parameters from Table 3 as input. We find $\zeta=224\pm
100$, with the uncertainty resulting from different parametrizations used to
interpolate the data.

Table 3: Input parameters for the calculation of the concentrated loss coefficient in the distributors, $\zeta$ Gas passage channel cross-section | $a^{*}$ | 0.00138 m2
---|---|---
Number of gas passage channels per distributor | $n_{P}$ | $7$
Column cross-section | $A$ |  $0.079$ m2
$A^{*}$/$A$= $n_{P}\times a^{*}$/$A$ |  | $0.12$

Table 4 shows the results obtained for the three runs, with the dominant
uncertainty coming from the determination of $\zeta$.

Table 4: Calculated velocity and density of the gas passing through the distributors, and pressure drop due to the distributors for runs A, B, and C from Eq. 7 |  $v_{\mathrm{G}}$(m s-1) |  $\rho_{\mathrm{G}}$(kg m-3) |  $\Delta p_{\mathrm{distributors}}$(mbar)
---|---|---|---
A | $0.202\pm 0.001$ | $9.40\pm 0.01$ | $1.7\pm 0.7$
B | $0.235\pm 0.001$ | $8.09\pm 0.01$ | $2.0\pm 0.9$
C | $0.242\pm 0.001$ | $9.97\pm 0.01$ | $2.6\pm 1.1$

####  4.1.2 Calculation of the pressure drop in the packing

The Sulcol and HYSYS programs are used to calculate the pressure drop per unit
length of the packing for the dry bed, i.e. in the absence of liquid, $\Delta
p_{\mathrm{packing,dry}}/\Delta z$ (the $z$ axis is directed downwards along
the column). These calculations were cross-checked with a model for pressure
drop calculation for packed columns described in Ref. stilch , using the input
parameters related to the packing reported in Table 14-14 of Ref. Perry . The
results are summarized in Table 5 and show that the three methods agree within
30%. The pressure drop in the packing increases with the liquid load since the
liquid impedes the passage of gas.

Table 5: Pressure drop per unit length of the packing for dry packing
calculated for runs A, B, and C. $\Delta p_{\mathrm{packing,dry}}/\Delta z$
(mbar m-1)
---
|  Sulcol |  HYSYS |  method of Ref. stilch
A |  $0.45$ |  $0.49$ |  $0.38$
B |  $0.52$ |  $0.56$ |  $0.44$
C |  $0.68$ |  $0.73$ |  $0.55$

The authors of Ref. stilch  also provide a method to estimate the pressure
drop of the irrigated bed. Table 6 summarizes the values obtained for the
pressure drop per unit length of the irrigated bed, $\Delta
p_{\mathrm{packing,irr}}/\Delta z$, which results about three times larger
than that for the dry bed. We assigned an uncertainty of 30% on this number
based on the comparisons of Table 5.

Table 6: Pressure drop per unit length for irrigated packing calculated for runs A, B and C using the method of Ref. stilch . |  $\Delta p_{\mathrm{packing,irr}}/\Delta z$ (mbar m-1)
---|---
A |  $1.1\pm 0.3$
B |  $1.2\pm 0.4$
C |  $1.6\pm 0.4$

The total pressure drop, $\Delta p_{\mathrm{tot,calc}}$, is given by the sum
of the pressure drop of the distributors, from Table 4, and that of the
irrigated packing, from Table 6 multiplied by the packing height of
$10.24\text{\,}\mathrm{m}$, and is shown in Table 7.

Table 7: Total experimental, $\Delta p_{\mathrm{tot,exp}}$, and calculated, $\Delta p_{\mathrm{tot,calc}}$, pressure drop for tests A, B, and C.  |  $\Delta p_{\mathrm{tot,exp}}$ (mbar) |  $\Delta p_{\mathrm{tot,calc}}$ (mbar)
---|---|---
A |  $16.3\pm 0.1$ |  $12.9\pm 4.0$
B |  $14.9\pm 0.2$ |  $14.3\pm 4.5$
C |  $19.6\pm 0.2$ |  $19.0\pm 6.0$

Table 7 also shows the measured pressure drop of the column, $\Delta
p_{\mathrm{tot,exp}}$, i.e., the difference between the top and bottom
recorded pressures (SIEMENS SITRANS P, DS III). Data and simulation show
agreement within an uncertainty of the order of 30%. We will assign a
corresponding uncertainty to the values of $\Delta
p_{\mathrm{packing,dry}}/\Delta z$, which will be used for comparison with
Sulzer data in Fig. 7. in relation to Table 7 and Table 5.

###  4.2 Liquid holdup

The liquid holdup (or retention), $\mathrm{h_{L}}$, is the fraction of the
packing volume occupied by the liquid. The holdup consists of two parts:
static hold-up and dynamic hold-up. The dominant one is the dynamic hold-up,
which depends on the liquid load, and, to some extent, on the gas load.

The holdup was evaluated using both HYSYS and Sulcol based on the operating
conditions of the column, and with a direct calculation based on measured and
estimated liquid volumes in the various elements of the plant, as described
below. Given the quantity of gas initially loaded, about $250$ kg, and the
operating conditions of the column, we calculate that there are about $10$ kg
of gas and $240$ kg of liquid, corresponding to $178$ L. The liquid is
distributed both in the packed sections and in the processing circuit. There
are about 80 L at the bottom of the column, as measured by a level indicator,
about 40 L in the reboiler, about 8 L in the condenser (the amount of liquid
argon is related to the amount of nitrogen in the auxiliary circuit, as
measured by a level indicator), about 11 L in the distributors (the height of
the liquid level in the distributors during distillation must be at least 3.5
cm, which is right above the height of the first hole in the distributor
pipes), and about 2 L in the piping. By subtraction, the amount of liquid in
the packed sections is approximately 37 L, about $20\%$ of the total liquid
present in the column. The stainless steel structure of the packing occupies
about 15% of the volume. The estimated liquid holdup is, therefore, between
5.5% and 6.0%, as shown in Table 8. This relatively low value compares well
with those obtained with HYSYS and Sulcol and agrees with measurements with
the same packing by others leveque .

Table 8: Liquid hold-up, $\mathrm{h_{L}}$, calculations for runs A, B, and C.
$\mathrm{h_{L}}(\%)$
---
|  Sulcol |  HYSYS |  our estimate
A |  $6.0$ |  $5.5$ |  $5.5$
B |  $6.0$ |  $5.5$ |  $5.5$
C |  $6.5$ |  $5.5$ |  $6.0$

##  5 Measurement of distillation parameters

As is well known, to achieve correct results from mass spectrometry in the
measurements of isotopic ratios, it is necessary to get calibration factors.
Therefore, we applied a time-dependent correction factor to the measured
isotopic ratios $(R_{36-40})_{T}$ and $(R_{36-40})_{B}$, given by the ratio of
the measured isotopic ratio in the feed, $(R_{36-40})_{F}$, and the known
value of the natural isotopic ratio of 0.334%. Figure 5 shows the corrected
isotopic ratios $(R_{36-40})_{T}$ and $(R_{36-40})_{B}$ as a function of time.
The correction also compensates for time variations in the response of the
sampling system and reduces them by about a factor of ten. A systematic
uncertainty on the isotopic ratios was evaluated, with reference to Fig. 3, by
comparing the measurement of the isotopic ratio using the peak height value
and the integral of the mass distributions. We find a $\leq 5$% discrepancy.

Figure 5: The isotopic ratios $(R_{36-40})_{T}$ (red) and $(R_{36-40})_{B}$
(blue) as a function of time, after correction using the measured isotopic
ratio in the feed, $(R_{36-40})_{F}$ (see text). The dashed-dotted line
corresponds to the natural isotopic ratio. The shaded regions correspond to
changing conditions in the column (see Fig. 1) and are not included in the
data analysis.

The measured separations $S^{TB}_{36-40}$ and $S^{TB}_{38-40}$ vs. time,
defined in Eq. 2 are shown for the three runs A, B, and C in Fig. 6. Time-
dependent effects in this ratio of isotopic ratios largely cancel out and lead
to negligible systematic uncertainty. The comparison of the measurement of
separation values with the peak height value and the integral of the mass
distributions has a maximum discrepancy of about 2%, of the same size of the
statistical uncertainty.

We measure one separation value every three hours, with occasional data points
not available due to connectivity problems with the mass spectrometer. In runs
A and B we observed a few % stability of $S^{TB}_{36-40}$, even though with a
limited number of measurements and duration of the run. In run C $\sim~{}7$%
fluctuations are observed, possibly due to small pressure variations in the
column during the run. We evaluated the impact of these fluctuations on
$\ln{\alpha_{36-40}}$ in A and found it negligible. Therefore, in the
following, we use the weighted average of all measured separation values
within a run.

Figure 6: Separation $S^{TB}_{36-40}$ (red) and $S^{TB}_{38-40}$ (green) as a
function of time. The shaded regions correspond to changing conditions in the
column (see Fig. 1) and are not included in the data analysis.

The time-averaged separations for each run, ${\cal S}^{TB}_{36-40}$ and ${\cal
S}^{TB}_{38-40}$, are shown in Table 9 and 10.

Table 9: Time-averaged separation ${\cal S}^{TB}_{36-40}$ during runs A, B, and C, logarithm of the relative volatility $\ln{\alpha_{36-40}}$, equivalent number of theoretical stages $N$ and HETP.  | ${\cal S}^{TB}_{36-40}$ | $\ln{\alpha_{36-40}}$ | $N$ | HETP (cm)
---|---|---|---|---
A | $1.49\pm 0.03$ | $(5.0\pm 0.4)\times 10^{-3}$ | $80\pm 7$ | $13\pm 1$
B | $1.51\pm 0.03$ | $(5.2\pm 0.4)\times 10^{-3}$ | $79\pm 7$ | $13\pm 1$
C | $1.48\pm 0.03$ | $(4.9\pm 0.4)\times 10^{-3}$ | $80\pm 8$ | $13\pm 1$
Table 10: Time-averaged separation ${\cal S}^{TB}_{38-40}$ during runs A, B, and C, logarithm of relative volatility $\ln{\alpha_{38-40}}$, equivalent number of theoretical stages $N$ and HETP. | ${\cal S}^{TB}_{38-40}$ | $\ln{\alpha_{38-40}}$ | $N$ | HETP (cm)
---|---|---|---|---
A | $1.20\pm 0.04$ | $(2.4\pm 0.3)\times 10^{-3}$ | $76\pm 17$ | $14\pm 3$
B | $1.21\pm 0.03$ | $(2.5\pm 0.3)\times 10^{-3}$ | $76\pm 13$ | $14\pm 2$
C | $1.20\pm 0.03$ | $(2.3\pm 0.3)\times 10^{-3}$ | $79\pm 15$ | $13\pm 2$

The corresponding relative volatility values were determined as discussed in
A, yielding an uncertainty of about 7% for $\ln{\alpha_{36-40}}$ and about 12%
for $\ln{\alpha_{38-40}}$. The equivalent number of theoretical stages $N$ and
the HETP were calculated from Eq. 1 and Eq. 5. The values of $N$ and HETP
obtained with the mass 36 and 38 isotopes and reported in the two tables agree
within measurement uncertainties. Given their rather strong correlation, in
the following, we only consider the measurement in Table 9.

Figure 7 shows the measured HETP (top) from Table 9 and the calculated $\Delta
p_{packing,dry}/\Delta z$ (bottom) (from Table 5) vs. the sizing parameter
$F_{G}$. $F_{G}$ is the F-factor or vapor load and is defined as $V_{G}$
$\times$ $\sqrt{\rho_{G}}$.

Figure 7:  Measured HETP (top), from Table 9 (red, green and violet dots), and
$\Delta p_{packing,dry}/\Delta z$ (bottom), from the third column of Table 5
assuming 30% uncertainty, vs. sizing parameter, $F_{G}$, for the three runs
with argon (red). Blue (cyan) line: measurements with
chlorobenzene/ethylbenzene mixtures at $p_{C}=0.96$ bar ($0.40$ bar), Sulzer
CY Gauze Packing, partial data-set extracted from the Structured Packing
brochure of Sulzer Ltd.. A perfect match between our data and those from
Sulzer is not expected due to the different thermodynamical parameters of the
fluids.

The three measured values of HETP are 20-30% larger than the expectation based
on the Sulzer data. This is possibly due to the different thermodynamical
parameters between the argon and the fluid used by Sulzer.

The measured separations and the derived HETP are conservative estimates due
to the limited data taking time (from 15 h to two days per run) at the same
conditions and the somewhat larger average gas pressure in the column compared
to the design value of 1.3 bar presented in Ref. Darkside:2021 .

The measured pressure drop per unit length, which we assume for the Sulzer
data to be the pressure drop per unit length for the dry bed, follows quite
well the curve measured by Sulzer, within uncertainties.

##  6 Multi-component distillation

In the paper on nitrogen distillation Darkside:2021 , we argued that the
distillation of ^39Ar would not be affected by the presence of the more
abundant ^36Ar and ^38Ar, based on the argument that for a gas mixture with
one dominant component and several other components with mass fractions below
a few 0.1% each, the distillation of these is essentially independent of the
presence of the others and is well described as a binary distillation with
respect to ^40Ar. Therefore, a multi-component approach for distillation of
argon Kister:1992  is not needed. This argument does not hold, e.g., for xenon
Back . This assumption was directly tested with this run.

In total reflux condition, if the distillation ^36Ar and ^38Ar proceed
independently of each other, we expect

$S^{TB}_{36-40}=(\alpha_{36-40})^{N}$ (8)

and

$S^{TB}_{38-40}=(\alpha_{38-40})^{N}$ (9)

From Eq.13 one infers that:

$k_{36-38}\coloneqq\frac{\mathrm{ln}(S^{TB}_{36-40})}{\mathrm{ln}(S^{TB}_{38-40})}\times\frac{1}{2.1}=1.$
(10)

Figure 8: Measured $k_{36-38}$ (as defined in Eq. 10) as a function of time.
The shaded regions correspond to changing conditions in the column (see Fig.
1) and are not included in the data analysis.

Fig. 8 shows the measured value of $k_{36-38}$ vs time. The time average of
these values is $0.998\pm 0.030$, consistent with 1, and demonstrates that no
multi-component approach is needed for argon distillation.

##  7 Comparison with simulation

The measured distillation data for run B were compared to simulations with
HYSYS. We used the conditions of run B as input parameters for the simulation.
These are the mean values of $p$ from Table 1, $V$ from Table 2, HETP from
Table 9, and ${\cal S}^{FB}_{36-40}$ from Table 11. The comparison of measured
and simulated separations is shown in the first two columns of Table 11. Good
agreement is observed within the uncertainties. The simulation was also
performed for binary systems using both the rigorous and the MCT methods. The
latter is included for completeness and because of its use in Ref.
Darkside:2021 . We find good agreement across the board.

Table 11: Comparison of separations between measurements and simulation, related to run B conditions: data (Exp.) from Table 9 and Table 10, HYSYS rigorous simulation with three components (Rig. 3), with two components (Rig. 2) and McCabe-Thiele (MCT) simulation. | Exp. | Rig. 3 | Rig. 2 | MCT
---|---|---|---|---
${\cal S}^{TB}_{36-40}$ | $1.51\pm 0.03$ | $1.5\pm 0.1$ | $1.5\pm 0.1$ | $1.5\pm 0.1$
${\cal S}^{TB}_{38-40}$ | $1.21\pm 0.03$ | $1.20\pm 0.04$ |  |

##  8 Conclusions and outlook

We validated the performance of the Aria plant with argon running the
prototype plant for a few days in total reflux. We significantly improved our
understanding of the plant performance compared to the previous nitrogen
distillation run, achieving a deeper understanding of operating
characteristics such as the hydraulic parameters and the behavior of multi-
component distillation. We also successfully modeled the Aria prototype
performance using commercial process simulation software. This study is a
validation of both the plant behavior and the software itself. Originally
designed for distilling organic compounds at room temperature, it is also
suitable for the simulation of cryogenic isotopic distillation.

The average measured HETP is about 30% larger than that measured in a
different column with the same packing, using organic mixtures at room
temperature by Sulzer. The discrepancy could be due to the different
thermodynamical parameters of the fluids used in the two cases or to the
somewhat larger average gas pressure in the column compared to the design
value of 1.3 bar presented in Ref. Darkside:2021 .

The measured pressure drop per unit length is in agreement, within
uncertainties, with that measured with the organic mixture.

## Acknowledgements

This report is based upon work supported by FSC 2014-2020 - Patto per lo
Sviluppo, Regione Sardegna, Italy, the U. S. National Science Foundation (NSF)
(Grants No. PHY-0919363, No. PHY-1004054, No. PHY-1004072, No. PHY-1242585,
No. PHY-1314483, No. PHY- 1314507, associated collaborative grants, No.
PHY-1211308, No. PHY-1314501, and No. PHY-1455351, as well as Major Research
Instrumentation Grant No. MRI-1429544), the Italian Istituto Nazionale di
Fisica Nucleare (Grants from Italian Ministero dell’Istruzione, Università, e
Ricerca Progetto Premiale 2013 and Commissione Scientific Nazionale II), the
Natural Sciences and Engineering Research Council of Canada, SNOLAB, and the
Arthur B. McDonald Canadian Astroparticle Physics Research Institute. We
acknowledge the financial support by LabEx UnivEarthS (ANR-10-LABX-0023 and
ANR18-IDEX-0001), the São Paulo Research Foundation (Grant
FAPESP-2017/26238-4), Chinese Academy of Sciences (113111KYSB20210030) and
National Natural Science Foundation of China (12020101004). The authors were
also supported by the Spanish Ministry of Science and Innovation (MICINN)
through the grant PID2019-109374GB-I00, the “Atraccion de Talento” Grant
2018-T2/ TIC-10494, the Polish NCN, Grant No. UMO- 2019/ 33/ B/ ST2/ 02884,
the Polish Ministry of Science and Higher Education, MNiSW, grant number
6811/IA/SP/2018, the International Research Agenda Programme AstroCeNT, Grant
No. MAB/2018/7, funded by the Foundation for Polish Science from the European
Regional Development Fund, the European Union’s Horizon 2020 research and
innovation program under grant agreement No 952480 (DarkWave), the Science and
Technology Facilities Council, part of the United Kingdom Research and
Innovation, and The Royal Society (United Kingdom), and IN2P3-COPIN consortium
(Grant No. 20-152). I.F.M.A is supported in part by Conselho Nacional de
Desenvolvimento Científico e Tecnológico (CNPq). We also wish to acknowledge
the support from Pacific Northwest National Laboratory, which is operated by
Battelle for the U.S. Department of Energy under Contract No. DE-
AC05-76RL01830. This research was supported by the Fermi National Accelerator
Laboratory (Fermilab), a U.S. Department of Energy, Office of Science, HEP
User Facility. Fermilab is managed by Fermi Research Alliance, LLC (FRA),
acting under Contract No. DE-AC02-07CH11359.

We acknowledge the professional contribution of the Mine and Electrical
Maintenance staff of Carbosulcis S.p.A. to this activity. We thank Polaris
S.r.L. and in particular M. Masetto and E.V.Canesi for their continuous
support during the preparation and execution of the run. We thank Fondazione
Aria for its contribution to the project and, in particular, for the help
during the data-taking. We thank S.Nisi of Istituto Nazionale di Fisica
Nucleare, Laboratori Nazionali del Gran Sasso, for his contribution in the
commissioning of the sampling system and M.Guetti of Istituto Nazionale di
Fisica Nucleare, Laboratori Nazionali del Gran Sasso, for the help in plant
maintenance. We thank M. Arba and M. Tuveri of the Cagliari Division of
Istituto Nazionale di Fisica Nucleare for their support during the preparation
of the run. We acknowledge the contribution of eng. S. Tosti, T. Pinna,
D.Dongiovanni, A.Santucci of ENEA for the safety analyses for Aria.

##  Appendix A Evaluation of relative volatilities

In this section, we give an update of our best estimate of the relative
volatilities $\alpha_{36-40}$ and $\alpha_{38-40}$, based on the scrutiny of
the literature, and we update the result derived in Appendix A of Ref.
Darkside:2021 . For $\mathrm{ln\alpha_{36-40}}$, we included measurements from
Ref. Bigeleisen:1972  and from Ref. Lee:1970a , corresponding to the second
definition of relative volatility of Eq. 3.

In Fig. 9 we report the measured dependence of $\mathrm{ln\alpha_{36-40}}$ on
temperature. The function $\mathrm{ln{\alpha}_{36-40}}=A\times 1/T^{2}+B$ is
fitted to the data of Ref. Boato1961  and Ref. Bigeleisen:1972 .

Figure 9: Measurements of $\mathrm{ln\alpha_{36-40}}$ vs. temperature. The
data are taken from Boato1961  (blue squares), Bigeleisen:1972  (violet
stars), and Lee:1970a  (green and red dots). The function
$\mathrm{ln{\alpha}_{36-40}}=A\times 1/T^{2}+B$ is fitted (continuous lines)
to the data of Ref. Boato1961  and Ref. Bigeleisen:1972 . The corresponding
one-sided standard deviation curves (dashed lines)are also shown. The derived
values of $\mathrm{ln\alpha_{36-40}}$ corresponding to the three runs A, B,
and C are also shown (black dots).

The choice of this parametrization follows the theoretical considerations of
Ref. Boato1962b . The errors on the single measurement were set all equal in
the fit and determined in retrospect requiring the reduced $\chi^{2}$ to be
one. Applying error propagation for the estimate of
$\mathrm{ln{\hat{\alpha}}_{36-40}}=\hat{A}\times 1/T^{2}+\hat{B}$ and of the
uncertainty $\sigma_{\mathrm{ln\hat{\alpha}}}$ as:

$\sigma_{\mathrm{ln\hat{\alpha}}}=\sqrt{V_{00}\times(1/T^{2})^{2}+V_{11}+2/T^{2}\times
V_{22}}$ (11)

with $V_{ij}$ being the elements of the covariance matrix, one obtains the
mean (continuous lines) and one-sided standard deviation curves (dashed
lines). To take into account all existing measurements in our estimate, we
assumed that the one standard deviation value of $\mathrm{ln\alpha_{36-40}}$
at a given temperature lies inside the two dashed curves. The values of
$\mathrm{ln\alpha_{36-40}}$ corresponding to the temperatures of the three
runs A, B, and C are also shown and their values are reported in Table 9. The
horizontal error bars correspond to the uncertainty on the temperature related
to the pressure variations during the runs of Table 1. They give a negligible
contribution to the uncertainty on $\mathrm{ln\alpha_{36-40}}$.

The relative volatility $\ln{\alpha_{38-40}}$, at given temperature, was
derived from $\ln{\alpha_{36-40}}$ using the theoretical formula from Ref.
CanongiaLopes:2003ju , Eq. (25):

$\ln{\alpha_{A-40}}\propto\frac{40-A}{A^{3}},$ (12)

with $A$ the atomic mass of an isotope, which entails:

$\ln{\alpha_{36-40}}\sim\ln{\alpha_{38-40}}\times 2.1.$ (13)

This dependence is also consistent, within errors, with the measurements of
Ref. Alamre:2020 , of $\ln{\alpha_{36-40}}$ = $\ln{\alpha_{38-40}}$
$\times(2.3\pm 0.2)$ at $87.3$ K. The relative uncertainty from these
measurements was propagated on the final estimate of the quantity
$\ln{\alpha_{38-40}}$.

## References

  * (1) P. Agnes, S. Albergo, I. F. M. Albuquerque, T. Alexander, A. Alici, A. K. Alton, P. Amaudruz, M. Arba, P. Arpaia, S. Arcelli. et al. (DarkSide-20k Collaboration), Eur. Phys. J. C 81, 4, 359 (2021).
  * (2) C.E. Aalseth, F. Acerbi, P. Agnes, I. F. M. Albuquerque, T. Alexander, A. Alici, A. K. Alton, P. Antonioli, S. Arcelli, R. Ardito et al. (DarkSide-20k Collaboration), Eur. Phys. J. Plus 133, 131 (2018)
  * (3) N. Abgrall, A. Abramov, N. Abrosimov, I. Abt, M. Agostini, M. Agartioglu, A. Ajjaq, S. I. Alvis, F. T. Avignone III, X. Bai, et al., AIP Conf. Proc. 1894, 020027 (2017)
  * (4) B. Abi, R. Acciarri, M. A. Acero, G. Adamov, D. Adams, M. Adinolfi, Z. Ahmad, J. Ahmed, T. Alion, S. Alonso Monsalve et al. Eur. Phys. J. C 80, 1 (2020)
  * (5) D. Akimov, P. An, C. Awe, P. S. Barbeau, B. Becker, V. Belov, M. A. Blackston, A. Bolozdynya, B. Cabrera-Palmer, N. Chen, et al, Phys. Rev. Lett. 126, 012002 (2021)
  * (6) J.-Y. Lee, K. Marti, J. P. Severinghaus, K. Kawamura H.-S. Yoo, J. Bok Lee, J. Seog Kim, Geochim. Cosmochim. Acta 70, 4507 (2006)
  * (7) , J.K. Böhlke, Pure Appl. Chem., 86 1421 (2014)
  * (8) H. Z. Kister, J. R. Haas, D. R. Hart, D. R. Gill. _Distillation Design_ , (McGraw-Hill, New York, 1992)
  * (9) M. Doninelli, M. Doninelli, Idraulica 28, 3 (2005)
  * (10) J. Stilchlmair, J. L. Bravo, J. R. Fair, Gas Sep. $\&$ Purif. 3, 19 (1989).
  * (11) D. W. Green, M.Z.Southard, _Perry’s chemical engineer’s handbook, 9th edition_ (McGraw-Hill, New York, 2019)
  * (12) J. Lévêque, D. Rouzineau, M. Prévost, M. Meyer, Chem. Eng. Sci. 64, 2607 (2009).
  * (13) H.O. Back, D.R. Bottenus, C. Clayton, D. Stephenson, W. TeGrotenhuis, J. Instrum. 12, 09033 (2017)
  * (14) G. Boato, G. Casanova, G. Scoles, M. E. Vallauri, Nuovo Cim. 20, 88 (1961)
  * (15) J. T Phillips, C. U. Linderstrom-Lang, J.Bigeleisen, J. Chem. Phys. 56, 5053 (1972)
  * (16) M. W. Lee, S. Fuks, J. Bigeleisen, J. Chem. Phys. 53, 4066 (1970)
  * (17) G. Boato, G. Casanova, A. Levi, J. Chem. Phys. 37, 201 (1962)
  * (18) J .N. Canongia Lopes, A.A.H. Pádua, L.P.N. Rebelo, J. Bigeleisen, J. Chem. Phys. 118, 5028 (2003)
  * (19) A. Alamre, I. Badhress, B. Death, C. Licciardi, D. Sinclair, ACS omega. 5, 28977 (2020)

The DarkSide-20k Collaboration

E. Aaron49, P. Agnes30, I. Ahmad66, S. Albergo2 ,3 , I. F. M. Albuquerque4, T.
Alexander5, A. K. Alton8, P. Amaudruz9, M. Atzori Corona28 ,21 , M. Ave4, I.
Ch. Avetisov10, O. Azzolini11, H. O. Back5, Z. Balmforth12, A. Barrado
Olmedo14, P. Barrillon15, A. Basco16, G. Batignani17 ,18 , V. Bocci41, W. M.
Bonivento21, B. Bottino22 ,23 , M. G. Boulay24, J. Busto15, M. Cadeddu21, A.
Caminata23, N. Canci16, A. Capra9, S. Caprioli23, M. Caravati21, N. Cargioli28
,21 , M. Carlini29, P. Castello32 ,21 , P. Cavalcante29, S. Cavuoti33 ,16 ,35
, S. Cebrian36, J. M. Cela Ruiz14, S. Chashin13, A. Chepurnov13, E
Chyhyrynets11, L. Cifarelli6 ,7 , D. Cintas36, M. Citterio38, B. Cleveland68
,69 , V. Cocco21, E. Conde Vilda14, L. Consiglio29, S. Copello23 ,22 , G.
Covone33 ,16 , M. Czubak37, M. D’Aniello34, S. D’Auria38, M. D. Da Rocha
Rolo39, S. Davini23, S. De Cecco41 ,42 , G. De Guido45, D. De Gruttola43 ,44 ,
S. De Pasquale43 ,44 , G. De Rosa33 ,16 , G. Dellacasa39, A. V. Derbin46, A.
Devoto28 ,21 , F. Di Capua33 ,16 , L. Di Noto23, P. Di Stefano73, G.
Dolganov47, F. Dordei21, E. Ellingwood73, T. Erjavec49, S. Farenzena95**$*$Not
a member of the DarkSide-20k Collaboration, M. Fernandez Diaz14, G. Fiorillo33
,16 , P. Franchini74 ,12 , D. Franco51, N. Funicello43 ,44 , F. Gabriele21, D.
Gahan28 ,21 , C. Galbiati52 ,29 ,30 , G. Gallina52, G. Gallus21 ,32 , M.
Garbini31 ,7 , P. Garcia Abia14, A. Gendotti53, C. Ghiano29, C. Giganti40, G.
K. Giovanetti55, V. Goicoechea Casanueva56, A. Gola57 ,58 , G. Grauso16, G.
Grilli di Cortona50, A. Grobov47 ,60 , M. Gromov13 ,61 , M. Guan62, M.
Guerzoni7, M. Gulino63 ,64 , C. Guo62, B. R. Hackett5, A. L. Hallin65, A.
Hamer92 ,12 , M. Haranczyk37, T. Hessel51, S. Hill12, S. Horikawa94 ,29 , F.
Hubaut15, J. Hucker73, T. Hugues66 ,51 , An. Ianni52 ,29 , V. Ippolito41, C.
Jillings68 ,69 , S. Jois12, P. Kachru30 ,29 , A. A. Kemp73, C. L. Kendziora67,
M. Kimura66, I. Kochanek29, K. Kondo29, G. Korga12, S. Koulosousas12, A.
Kubankin70, M. Kuss17, M. Kuźniak66, M. La Commara71 ,16 , M. Lai28 ,21 , N.
Lami95**footnotemark: $*$, E. Le Guirriec15, E. Leason12, A. Leoni29, L.
Lidey5, F. Lippi95**footnotemark: $*$, M. Lissia21, L. Luzzi14, O.
Lychagina61, N. Maccioni95**footnotemark: $*$, O. Macfadyen12, I. N.
Machulin47 ,60 , S. Manecki68 ,69 , I. Manthos89, L. Mapelli52, A. Margotti7,
S. M. Mari26 ,27 , C. Mariani80, J. Maricic56, A. Marini22 ,23 , M. Martínez36
,72 , C. J. Martoff91, M. Mascia84, A. Masoni21, G. Matteucci33 ,16 , K.
Mavrokoridis88, C. Maxia95**footnotemark: $*$, A. B. McDonald73, A. Messina41
,42 , R. Milincic56, A. Mitra93, A. Moharana30 ,29 , S. Moioli45, J. Monroe12,
E. Moretti57, M. Morrocchi17 ,18 , T. Mróz37, V. N. Muratova46, C. Muscas32
,21 , P. Musico23, R. Nania7, M. Nessi29, K. Nikolopoulos89, J. Nowak74, K.
Olchansky9, A. Oleinik70, V. Oleynikov19 ,20 , P. Organtini52 ,29 , A. Ortiz
de Solórzano36, L. Pagani49, M. Pallavicini22 ,23 , L. Pandola64, E. Pantic49,
E. Paoloni17 ,18 , G. Paternoster57 ,58 , P. A. Pegoraro32 ,21 , K. Pelczar37,
L. A. Pellegrini45, C. Pellegrino7, V. Pesudo14, S. Piacentini42 ,41 , L.
Pietrofaccia29, N. Pino2 ,3 , A. Pocar48, D. M. Poehlmann49, S. Pordes67, P.
Pralavorio15, D. Price76, F. Ragusa77 ,38 , Y. Ramachers93, M. Razeti21, A. L.
Renshaw1, M. Rescigno41, F. Retiere9, L. P. Rignanese7 ,6 , C. Ripoli44 ,43 ,
A. Rivetti39, A. Roberts88, C. Roberts88, J. Rode40 ,51 , G. Rogers89, L.
Romero14, M. Rossi23 ,22 , A. Rubbia53, M. A. Sabia42 ,41 , G. M.
Sabiu95**footnotemark: $*$, P. Salomone42 ,41 , E. Sandford76, S.
Sanfilippo64, D. Santone12, R. Santorelli14, C. Savarese52, E. Scapparone7, G.
Schillaci64, F. Schukman73, G. Scioli6 ,7 , M. Simeone79 ,16 , P. Skensved73,
M. D. Skorokhvatov47 ,60 , O. Smirnov61, T. Smirnova47, B. Smith9, F.
Spadoni5, M. Spangenberg93, R. Stefanizzi28 ,21 , A. Steri21, V. Stornelli94
,29 , S. Stracka17, M. Stringer73, S. Sulis32 ,21 , A. Sung52, Y. Suvorov33
,16 ,47 , A. M. Szelc92, R. Tartaglia29, A. Taylor88, J. Taylor88, S.
Tedesco39 ,54 , G. Testera23, K. Thieme56, T. N. Thorpe81, A. Tonazzo51, A.
Tricomi2 ,3 , E. V. Unzhakov46, T. Vallivilayil John30 ,29 , M. Van Uffelen15,
T. Viant53, S. Viel24, R. B. Vogelaar80, J. Vossebeld88, M. Wada66 ,28 , M. B.
Walczak66, H. Wang81, Y. Wang62 ,85 , S. Westerdale87 ,52 , L. Williams82, I.
Wingerter-Seez15, R. Wojaczyński66, Ma. M. Wojcik37, T. Wright80, Y. Xie62 ,85
, C. Yang62 ,85 , A. Zabihi66, P. Zakhary66, A. Zani38, A. Zichichi6 ,7 , G.
Zuzel37, M. P. Zykova10

  1. 1

Department of Physics, University of Houston, Houston, TX 77204, USA

  2. 2

INFN Catania, Catania 95121, Italy

  3. 3

Università of Catania, Catania 95124, Italy

  4. 4

Instituto de Física, Universidade de São Paulo, São Paulo 05508-090, Brazil

  5. 5

Pacific Northwest National Laboratory, Richland, WA 99352, USA

  6. 6

Department of Physics and Astronomy, Università degli Studi di Bologna,
Bologna 40126, Italy

  7. 7

INFN Bologna, Bologna 40126, Italy

  8. 8

Physics Department, Augustana University, Sioux Falls, SD 57197, USA

  9. 9

TRIUMF, 4004 Wesbrook Mall, Vancouver, BC V6T 2A3, Canada

  10. 10

Mendeleev University of Chemical Technology, Moscow 125047, Russia

  11. 11

INFN Laboratori Nazionali di Legnaro, Legnaro (Padova) 35020, Italy

  12. 12

Department of Physics, Royal Holloway University of London, Egham TW20 0EX, UK

  13. 13

Skobeltsyn Institute of Nuclear Physics, Lomonosov Moscow State University,
Moscow 119234, Russia

  14. 14

CIEMAT, Centro de Investigaciones Energéticas, Medioambientales y
Tecnológicas, Madrid 28040, Spain

  15. 15

Centre de Physique des Particules de Marseille, Aix Marseille Univ,
CNRS/IN2P3, CPPM, Marseille, France

  16. 16

INFN Napoli, Napoli 80126, Italy

  17. 17

INFN Pisa, Pisa 56127, Italy

  18. 18

Physics Department, Università degli Studi di Pisa, Pisa 56127, Italy

  19. 19

Budker Institute of Nuclear Physics, Novosibirsk 630090, Russia

  20. 20

Novosibirsk State University, Novosibirsk 630090, Russia

  21. 21

INFN Cagliari, Cagliari 09042, Italy

  22. 22

Physics Department, Università degli Studi di Genova, Genova 16146, Italy

  23. 23

INFN Genova, Genova 16146, Italy

  24. 24

Department of Physics, Carleton University, Ottawa, ON K1S 5B6, Canada

  25. 25

CERN, European Organization for Nuclear Research 1211 Geneve 23, Switzerland,
CERN

  26. 26

INFN Roma Tre, Roma 00146, Italy

  27. 27

Mathematics and Physics Department, Università degli Studi Roma Tre, Roma
00146, Italy

  28. 28

Physics Department, Università degli Studi di Cagliari, Cagliari 09042, Italy

  29. 29

INFN Laboratori Nazionali del Gran Sasso, Assergi (AQ) 67100, Italy

  30. 30

Gran Sasso Science Institute, L’Aquila 67100, Italy

  31. 31

Museo Storico della Fisica e Centro Studi e Ricerche Enrico Fermi, Roma 00184,
Italy

  32. 32

Department of Electrical and Electronic Engineering, Università degli Studi di
Cagliari, Cagliari 09123, Italy

  33. 33

Physics Department, Università degli Studi “Federico II” di Napoli, Napoli
80126, Italy

  34. 34

Department of Strutture per l’Ingegneria e l’Architettura, Università degli
Studi “Federico II” di Napoli, Napoli 80131, Italy

  35. 35

INAF Osservatorio Astronomico di Capodimonte, 80131 Napoli, Italy

  36. 36

Centro de Astropartículas y Física de Altas Energías, Universidad de Zaragoza,
Zaragoza 50009, Spain

  37. 37

M. Smoluchowski Institute of Physics, Jagiellonian University, 30-348 Krakow,
Poland

  38. 38

INFN Milano, Milano 20133, Italy

  39. 39

INFN Torino, Torino 10125, Italy

  40. 40

LPNHE, CNRS/IN2P3, Sorbonne Université, Université Paris Diderot, Paris 75252,
France

  41. 41

INFN Sezione di Roma, Roma 00185, Italy

  42. 42

Physics Department, Sapienza Università di Roma, Roma 00185, Italy

  43. 43

Physics Department, Università degli Studi di Salerno, Salerno 84084, Italy

  44. 44

INFN Gruppo Collegato di Salerno, Sezione di Napoli, Salerno 84084, Italy

  45. 45

Chemistry, Materials and Chemical Engineering Department “G. Natta",
Politecnico di Milano, Milano 20133, Italy

  46. 46

Saint Petersburg Nuclear Physics Institute, Gatchina 188350, Russia

  47. 47

National Research Centre Kurchatov Institute, Moscow 123182, Russia

  48. 48

Amherst Center for Fundamental Interactions and Physics Department, University
of Massachusetts, Amherst, MA 01003, USA

  49. 49

Department of Physics, University of California, Davis, CA 95616, USA

  50. 50

INFN Laboratori Nazionali di Frascati, Frascati 00044, Italy

  51. 51

APC, Université de Paris, CNRS, Astroparticule et Cosmologie, Paris F-75013,
France

  52. 52

Physics Department, Princeton University, Princeton, NJ 08544, USA

  53. 53

Institute for Particle Physics, ETH Zürich, Zürich 8093, Switzerland

  54. 54

Department of Electronics and Communications, Politecnico di Torino, Torino
10129, Italy

  55. 55

Williams College, Physics Department, Williamstown, MA 01267 USA

  56. 56

Department of Physics and Astronomy, University of Hawai’i, Honolulu, HI
96822, USA

  57. 57

Fondazione Bruno Kessler, Povo 38123, Italy

  58. 58

Trento Institute for Fundamental Physics and Applications, Povo 38123, Italy

  59. 59

Universiatat de Barcelona, Barcelona E-08028, Catalonia, Spain

  60. 60

National Research Nuclear University MEPhI, Moscow 115409, Russia

  61. 61

Joint Institute for Nuclear Research, Dubna 141980, Russia

  62. 62

Institute of High Energy Physics, Beijing 100049, China

  63. 63

Engineering and Architecture Faculty, Università di Enna Kore, Enna 94100,
Italy

  64. 64

INFN Laboratori Nazionali del Sud, Catania 95123, Italy

  65. 65

Department of Physics, University of Alberta, Edmonton, AB T6G 2R3, Canada

  66. 66

AstroCeNT, Nicolaus Copernicus Astronomical Center of the Polish Academy of
Sciences, 00-614 Warsaw, Poland

  67. 67

Fermi National Accelerator Laboratory, Batavia, IL 60510, USA

  68. 68

SNOLAB, Lively, ON P3Y 1N2, Canada

  69. 69

Department of Physics and Astronomy, Laurentian University, Sudbury, ON P3E
2C6, Canada

  70. 70

Radiation Physics Laboratory, Belgorod National Research University, Belgorod
308007, Russia

  71. 71

Pharmacy Department, Università degli Studi “Federico II” di Napoli, Napoli
80131, Italy

  72. 72

Fundación ARAID, Universidad de Zaragoza, Zaragoza 50009, Spain

  73. 73

Department of Physics, Engineering Physics and Astronomy, Queen’s University,
Kingston, ON K7L 3N6, Canada

  74. 74

Physics Department, Lancaster University, Lancaster LA1 4YB, UK

  75. 75

Civil and Environmental Engineering Department, Politecnico di Milano, Milano
20133, Italy

  76. 76

Department of Physics and Astronomy, The University of Manchester, Manchester
M13 9PL, UK

  77. 77

Physics Department, Università degli Studi di Milano, Milano 20133, Italy

  78. 78

Brookhaven National Laboratory, Upton, NY 11973, USA

  79. 79

Chemical, Materials, and Industrial Production Engineering Department,
Università degli Studi “Federico II” di Napoli, Napoli 80126, Italy

  80. 80

Virginia Tech, Blacksburg, VA 24061, USA

  81. 81

Physics and Astronomy Department, University of California, Los Angeles, CA
90095, USA

  82. 82

Department of Physics and Engineering, Fort Lewis College, Durango, CO 81301,
USA

  83. 83

Institute of Applied Radiation Chemistry, Lodz University of Technology,
93-590 Lodz, Poland

  84. 84

Department of Mechanical, Chemical, and Materials Engineering, Università
degli Studi, Cagliari 09042, Italy

  85. 85

University of Chinese Academy of Sciences, Beijing 100049, China

  86. 86

Institut für Kernphysik, Johannes Gutenberg-Universität Mainz, D-55128 Mainz,
Germany

  87. 87

Department of Physics and Astronomy, University of California, Riverside, CA
92507, USA

  88. 88

Department of Physics, University of Liverpool, The Oliver Lodge Laboratory,
Liverpool L69 7ZE, UK

  89. 89

School of Physics and Astronomy, University of Birmingham, Edgbaston, B15 2TT,
Birmingham, UK

  90. 90

Science & Technology Facilities Council (STFC), Rutherford Appleton
Laboratory, Technology, Harwell Oxford, Didcot OX11 0QX, UK

  91. 91

Physics Department, Temple University, Philadelphia, PA 19122, USA

  92. 92

School of Physics and Astronomy, University of Edinburgh, Edinburgh EH9 3FD,
UK

  93. 93

University of Warwick, Department of Physics, Coventry CV47AL, UK

  94. 94

Università degli Studi dell’Aquila, L’Aquila 67100, Italy

  95. 95

CarboSulcis S.p.A. - Miniera Monte Sinni, 09010 Gonnesa, Italy


We now prove the limiting distribution of
$(N_{[rt]}+\Upsilon_{[rt]})/\sqrt{t}$ in following steps.

We first check $s_{t}^{-2}V_{t}^{2}\stackrel{{\scriptstyle
p}}{{\to}}\eta^{2}>0$ $a.s.$. By (Li et al.,, 2022) (page 55),
$\frac{1}{t}\sum_{i=B}^{t}\mathbbm{E}(\bm{\nu}_{i}\bm{\nu}_{i}^{T}|\mathcal{F}_{i-1})\stackrel{{\scriptstyle
p}}{{\to}}S+\mathbbm{E}(\Xi_{\bm{x}}(\bm{\theta}^{*})\Lambda\Xi_{\bm{x}}(\bm{\theta}^{*}))$,
then

$\displaystyle\frac{1}{t}V_{t}^{2}=\frac{1}{t}\sum_{i=1}^{t}\mathbbm{E}(X_{i}^{2}|\mathcal{F}_{i-1})=\bm{w}^{T}\left[\frac{1}{t}\sum_{i=B}^{t}\mathbbm{E}(\bm{\nu}_{i}\bm{\nu}_{i}^{T}|\mathcal{F}_{i-1})\right]\bm{w}$
$\displaystyle\qquad\stackrel{{\scriptstyle
d}}{{\to}}\bm{w}^{T}(S+\mathbbm{E}(\Xi_{\bm{x}}(\bm{\theta}^{*})\Lambda\Xi_{\bm{x}}(\bm{\theta}^{*})))\bm{w},$
$\displaystyle\frac{1}{t}s_{t}^{2}=\mathbbm{E}\left[\frac{1}{t}V_{t}^{2}\right]\to\bm{w}^{T}(S+\mathbbm{E}(\Xi_{\bm{x}}(\bm{\theta}^{*})\Lambda\Xi_{\bm{x}}(\bm{\theta}^{*})))\bm{w}.$

Since
$C:=\bm{w}^{T}(S+\mathbbm{E}(\Xi_{\bm{x}}(\bm{\theta}^{*})\Lambda\Xi_{\bm{x}}(\bm{\theta}^{*})))\bm{w}>0$
is a constant, we also have $\frac{1}{t}V_{t}^{2}$ converges in probability
implied by comvergence in distribution

$\frac{1}{t}V_{t}^{2}\stackrel{{\scriptstyle
p}}{{\to}}\bm{w}^{T}(S+\mathbbm{E}(\Xi_{\bm{x}}(\bm{\theta}^{*})\Lambda\Xi_{\bm{x}}(\bm{\theta}^{*})))\bm{w}.$

Thus

$s_{t}^{-2}V_{t}^{2}\stackrel{{\scriptstyle p}}{{\to}}1.$

We then check the conditional Lindeberg condition. By (Li et al.,, 2022) (page
56), one have that as $t\to\infty$,

$\displaystyle\forall\epsilon^{\prime}>0,R_{t}(\epsilon^{\prime}):=t^{-1}\sum_{i=1}^{t}\mathbbm{E}\left[\|\nu_{i}\|_{2}^{2}I(\|\nu_{i}\|_{2}>\epsilon^{\prime}\sqrt{t})|\mathcal{F}_{i-1}\right]\stackrel{{\scriptstyle
p}}{{\to}}0.$ (53)

To check the Lindeberg condition for $X_{i}$ sequence for any $\epsilon>0$, we
note that

$\displaystyle
s_{t}^{-2}\sum_{i=1}^{t}\mathbbm{E}\left[X_{i}^{2}I(|X_{i}|>\epsilon
s_{t})|\mathcal{F}_{i-1}\right]\to\frac{1}{tC}\sum_{i=1}^{t}\mathbbm{E}\left[X_{i}^{2}I(|X_{i}|>\epsilon\sqrt{Ct})|\mathcal{F}_{i-1}\right],$

where

$\displaystyle\frac{1}{t}\sum_{i=1}^{t}\mathbbm{E}\left[X_{i}^{2}I(|X_{i}|>\epsilon\sqrt{Ct})|\mathcal{F}_{i-1}\right]$
$\displaystyle\leq$
$\displaystyle\frac{1}{t}\sum_{i=1}^{t}\mathbbm{E}\left[\|\bm{\nu}_{i}\|_{2}^{2}\|\bm{w}\|_{2}^{2}I(\|\bm{\nu}_{i}\|_{2}>\epsilon\sqrt{Ct}/\|\bm{w}\|_{2})|\mathcal{F}_{i-1}\right]$
$\displaystyle\to$ $\displaystyle 0\left(\text{By taking
}\epsilon^{\prime}=\frac{\epsilon C}{\|\bm{w}\|_{2}}\text{ in
\eqref{eq:lindeberg_rootsgd}}\right).$

Thus the Lindeberg condition holds.

We now have all the conditions in Lemma 7 hold, thus we have $\zeta_{t}(r)$
satisfies that

$\zeta_{t}(r)\stackrel{{\scriptstyle d}}{{\to}}W_{r},$

where
$\zeta_{t}(r)=V_{t}^{-1}\\{S_{i}+(V_{i+1}^{2}-V_{i}^{2})^{-1}(rV_{t}^{2}-V_{i}^{2})X_{i+1}\\}$
and $i$ is such that $V_{i}^{2}\leq rV_{t}^{2}<V_{i+1}^{2}$.

On the other hand, ${V_{[rt]}^{2}}/{V_{t}^{2}}\to r$. Thus
$\zeta_{t}\left({V_{[rt]}^{2}}/{V_{t}^{2}}\right)\to\zeta_{t}(r)$. That is,
$V_{t}^{-1}S_{[rt]}=\zeta_{t}\left({V_{[rt]}^{2}}/{V_{t}^{2}}\right)\stackrel{{\scriptstyle
d}}{{\to}}W_{r}$, then by Slutsky’s theorem
$S_{[rt]}/\sqrt{t}\stackrel{{\scriptstyle d}}{{\to}}\sqrt{C}W_{r}$. In this
way,

$\bm{w}^{T}(N_{[rt]}+\Upsilon_{[rt]})/\sqrt{t}\to
S_{[rt]}/\sqrt{t}\stackrel{{\scriptstyle
d}}{{\to}}\sqrt{\bm{w}^{T}(S+\mathbbm{E}(\Xi_{\bm{x}}(\bm{\theta}^{*})\Lambda\Xi_{\bm{x}}(\bm{\theta}^{*})))\bm{w}}W_{r}.$

Recall that the above holds for any $\bm{w}\neq\bm{0}$, so we have

$(N_{[rt]}+\Upsilon_{[rt]})/\sqrt{t}\stackrel{{\scriptstyle
d}}{{\to}}(S+\mathbbm{E}(\Xi_{\bm{x}}(\bm{\theta}^{*})\Lambda\Xi_{\bm{x}}(\bm{\theta}^{*})))^{1/2}\bm{W}_{r},$

which is the limiting distribution as claimed. ∎

##  Appendix F Proof for Theorem 2

In this section, we prove Theorem 2 using Proposition 2.

###### Proof.

By Proposition 2, we have that the random function
$C_{t}(r)=r\sqrt{t}\bm{w}^{T}(\widehat{\bm{\theta}}_{[rt]}-\bm{\theta}^{*})$
satisfies

$C_{t}(r)\stackrel{{\scriptstyle
d}}{{\to}}(\bm{w}^{T}\Sigma\bm{w})^{1/2}W_{r}.$ (54)

Our statistic is

$\displaystyle\frac{\sqrt{t}\bm{w}^{T}(\widehat{\bm{\theta}}_{t}-\bm{\theta}^{*})}{\sqrt{\bm{w}^{T}V_{t}\bm{w}}}$
$\displaystyle=$
$\displaystyle\frac{C_{t}(1)}{\sqrt{\frac{1}{t}\sum_{i=1}^{t}\frac{i^{2}}{t}(\bm{w}^{T}(\widehat{\bm{\theta}}_{i}-\widehat{\bm{\theta}}_{t}))^{2}}}$
$\displaystyle\stackrel{{\scriptstyle t\to\infty}}{{\to}}$
$\displaystyle\frac{C_{t}(1)}{\sqrt{\int_{0}^{1}r^{2}t(\bm{w}^{T}(\widehat{\bm{\theta}}_{[rt]}-\widehat{\bm{\theta}}_{t}))^{2}dr}}$
$\displaystyle=$
$\displaystyle\frac{C_{t}(1)}{\sqrt{\int_{0}^{1}(C_{t}(r)-rC_{t}(1))^{2}dr}}.$

Now $\frac{C_{t}(1)}{\sqrt{\int_{0}^{1}(C_{t}(r)-rC_{t}(1))^{2}dr}}$ is a
continuous function of $C_{t}(\cdot)$, so by continuous mapping theorem
(Theorem 18.11 of (Van der Vaart,, 2000)) we have that

$\frac{C_{t}(1)}{\sqrt{\int_{0}^{1}(C_{t}(r)-rC_{t}(1))^{2}dr}}\stackrel{{\scriptstyle
d}}{{\to}}\frac{(\bm{w}^{T}\Sigma\bm{w})^{1/2}W_{1}}{(\bm{w}^{T}\Sigma\bm{w})^{1/2}\sqrt{\int_{0}^{1}(W_{r}-rW_{1})^{2}dr}}=\frac{W_{1}}{\sqrt{\int_{0}^{1}(W_{r}-rW_{1})^{2}dr}}.$

This gives us the theorem claim. ∎

##  Appendix G Comparison of Our Covariance Estimators with Those for the ASGD
Algorithm

In this section, we compare our plug-in covariance estimator and the random-
scaling estimator with their counterparts for the ASDG algorithm.

###  G.1 Plug-in Estimator Comparison

We compare our plug-in estimator with the plug-in estimator for ASGD in (Chen
et al.,, 2020) from two aspects: the convergence rate and the computational
burden.

For the convergence rate comparison, our plug-in estimator converges faster to
the true asymptotic covariance of ROOT-SGD compared to the ASGD counterpart.
Recall that for ASGD, the plug-in covariance estimator
$\widehat{\Sigma}_{t,ASGD}$ converges to the true asymptotic covariance
${\Sigma}_{ASGD}$ as
$\mathbbm{E}\|\Sigma_{ASGD}-\widehat{\Sigma}_{t,ASGD}\|_{2}\lesssim
t^{-\alpha/2}$. Since $\alpha\in(1/2,1)$, the convergence rate is strictly
slower than $\mathcal{O}(t^{-1/2})$. Our plug-in covariance estimator
$\widehat{\Sigma}_{t}$ converges to the true asymptotic covariance ${\Sigma}$
of ROOT-SGD as $\mathbbm{E}\|\Sigma-\widehat{\Sigma}_{t}\|_{2}\lesssim
t^{-1/2}$. So the convergence speed of our estimator is strictly faster than
that of the plug-in estimator in (Chen et al.,, 2020). Moreover, this
$\mathcal{O}(t^{-1/2})$ rate matches the optimal statistical rate in such a
random sampling scheme.

For the computational aspect of the plug-in estimator, both our estimator and
that in (Chen et al.,, 2020) can be computed in a fully online fashion.
Suppose the algorithm is updated for $t$ steps, and the parameter is of
dimension $p$. Our plug-in estimator takes $\mathcal{O}(p^{6}+tp^{4})$ total
arithmetic computation, while their estimator takes
$\mathcal{O}(p^{3}+tp^{2})$ total arithmetic computation. Our plug-in
estimator takes more computation. This is due to the intrinsic structure of
the asymptotic covariance of the ROOT-SGD estimator: the asymptotic covariance
depends on the quantity

$\mathbbm{E}(\nabla^{2}f(\bm{\theta}^{*};\bm{x})\Lambda\nabla^{2}f(\bm{\theta}^{*};\bm{x}));$
(55)

moreover, one needs to solve an equation that contains (55) in $\Lambda$.
Thus, in our plug-in estimator, we need to evaluate the empirical counter-part
of (55) for an unknown $\Lambda$, which is
$\widehat{P}(\Lambda)=\frac{1}{t}\sum_{i=1}^{t}\nabla^{2}f(\widehat{\bm{\theta}}_{i-1};\bm{x}_{i})\Lambda\nabla^{2}f(\widehat{\bm{\theta}}_{i-1};\bm{x}_{i})$.
To achieve this, we keep an online update of
$\widehat{P}=\frac{1}{t}\sum_{i=1}^{t}\nabla^{2}f(\widehat{\bm{\theta}}_{i-1};\bm{x}_{i})\otimes\nabla^{2}f(\widehat{\bm{\theta}}_{i-1};\bm{x}_{i})$.
We further invert $\widehat{P}$ when solving the empirical counterpart of
$\Lambda$. Computing $\widehat{P}\in\mathcal{R}^{p^{2}\times p^{2}}$ and its
inverse then needs $\mathcal{O}(tp^{4}+p^{6})$ operations, which dominates the
computation. As a comparison, the asymptotic covariance of the ASGD estimator
does not contain any $\Lambda$ term; thus, its plug-in estimator is computed
with fewer computations.

The computation burden of our plug-in estimator for ROOT-SGD is the same as
that for ASGD for the special case of
$\nabla^{2}f(\bm{\theta}^{*};\bm{x})\equiv\nabla^{2}F(\bm{\theta}^{*})$. Such
a special case holds for the exponential family model. In this case, the
asymptotic covariance of ROOT-SGD reduces to $A^{-1}SA^{-1}$, and the plug-in
covariance of ROOT-SGD is then $\widehat{A}^{-1}\widehat{S}\widehat{A}^{-1}$.
Then the computational complexity of our plug-in estimator becomes
$\mathcal{O}(p^{3}+tp^{2})$.

###  G.2 Hessian-free Estimator Comparison

We compare the computation of our random-scaling estimator with the random-
scaling estimator for ASGD in (Lee et al.,, 2021).

Both estimators can be computed fully online. For $t$ steps of the algorithm
on a problem of dimension $p$, our random-scaling estimator takes
$\mathcal{O}(tp^{2})$ arithmetic computations, which is the same as the
random-scaling estimator for ASGD and is less than our plug-in estimator.

As for the asymptotic convergence of the random-scaling estimators, both our
estimator and that in (Lee et al.,, 2021) are asymptotically consistent.
Unfortunately, there is no convergence rate result in either our work or (Lee
et al.,, 2021). Though the convergence rates of the random-scaling estimators
are not guaranteed, considering that ROOT-SGD converges faster than SGD, there
is still an advantage to use ROOT-SGD and our random-scaling estimator as
compared to the ASGD counterpart. And we do see such an advantage in our
experiment.