# High-precision mass measurement of 24Si and a refined determination of the
$rp$ process at the $A=22$ waiting point

D. Puentes<EMAIL_ADDRESS>Department of Physics and Astronomy, Michigan
State University, East Lansing, Michigan 48824, USA  National Superconducting
Cyclotron Laboratory, East Lansing, Michigan 48824, USA  Z. Meisel  Institute
of Nuclear & Particle Physics, Department of Physics and Astronomy, Ohio
University, Athens, Ohio 45701, USA  Joint Institute for Nuclear Astrophysics
– Center for the Evolution of the Elements, Michigan State University, East
Lansing, Michigan 48824, USA  G. Bollen  Department of Physics and Astronomy,
Michigan State University, East Lansing, Michigan 48824, USA  Facility for
Rare Isotope Beams, East Lansing, Michigan 48824, USA  A. Hamaker  Department
of Physics and Astronomy, Michigan State University, East Lansing, Michigan
48824, USA  National Superconducting Cyclotron Laboratory, East Lansing,
Michigan 48824, USA  C. Langer  University of Applied Sciences Aachen, 52066
Aachen, Germany  E. Leistenschneider  National Superconducting Cyclotron
Laboratory, East Lansing, Michigan 48824, USA  C. Nicoloff  Department of
Physics and Astronomy, Michigan State University, East Lansing, Michigan
48824, USA  National Superconducting Cyclotron Laboratory, East Lansing,
Michigan 48824, USA  W.-J. Ong  Nuclear and Chemical Sciences Division,
Lawrence Livermore National Laboratory, Livermore, California 94550, USA
Joint Institute for Nuclear Astrophysics – Center for the Evolution of the
Elements, Michigan State University, East Lansing, Michigan 48824, USA  M.
Redshaw  National Superconducting Cyclotron Laboratory, East Lansing, Michigan
48824, USA  Department of Physics, Central Michigan University, Mount
Pleasant, Michigan 48859, USA  R. Ringle  National Superconducting Cyclotron
Laboratory, East Lansing, Michigan 48824, USA  C. S. Sumithrarachchi  National
Superconducting Cyclotron Laboratory, East Lansing, Michigan 48824, USA  J.
Surbrook  Department of Physics and Astronomy, Michigan State University, East
Lansing, Michigan 48824, USA  National Superconducting Cyclotron Laboratory,
East Lansing, Michigan 48824, USA  A. A. Valverde  Department of Physics and
Astronomy, University of Manitoba, Winnipeg, MB R3T 2N2, Canada  I. T. Yandow
Department of Physics and Astronomy, Michigan State University, East Lansing,
Michigan 48824, USA  National Superconducting Cyclotron Laboratory, East
Lansing, Michigan 48824, USA

###### Abstract

We report a high precision mass measurement of ${}^{24}{\rm Si}$, performed
with the LEBIT facility at the National Superconducting Cyclotron Laboratory.
The atomic mass excess, $10\;753.8$(37) keV, is a factor of 5 more precise
than previous results. This substantially reduces the uncertainty of the
${}^{23}{\rm Al}(p,\gamma)^{24}{\rm Si}$ reaction rate, which is a key part of
the rapid proton capture ($rp$) process powering Type I X-ray bursts. The
updated rate constrains the onset temperature of the $(\alpha,p)$ process at
the ${}^{22}{\rm Mg}$ waiting-point to a precision of 9%.

Type I X-ray bursts occur at astrophysical sites consisting of a binary star
system where a neutron star accretes H/He-rich matter from a companion star,
leading to nuclear burning on the neutron star surface [1, 2]. As the
temperature increases, the hot CNO cycle proceeds until the triple-$\alpha$
process is initiated with an ensuing thermonuclear runaway. The thermonuclear
runaway explosion consists of multiple competing reaction sequences including
the rapid proton ($rp$) [3] and ($\alpha,p$) [4] processes. The rate at which
the runaway occurs as well as the reaction sequence competition depends on
factors such as accretion rate and fuel composition. The only observable from
these explosions, however, is the X-ray burst light curve, which can be
recorded with telescopes like the Neutron Star Interior Composition Explorer
(NICER) [5].

In addition to the astrophysical conditions, the main determinant of the shape
of the light curve is the nuclear physics involved [6]. The nuclear physics
input consists of the fast ($p,\gamma$), ($\alpha,\gamma$) and ($\alpha,p$)
reaction rates and slow $\beta^{+}$ decays. Single-zone [7] and multi-zone [8,
9] X-ray burst models have been used to understand the impact of the
uncertainties derived from these nuclear physics inputs [10, 11]. Through
these analyses, it was shown that only a few reactions play a large role in
the uncertainties of the simulated light curves [10]. The dramatic variations
in the simulated light curve due to these reaction rate uncertainties can make
it difficult to extract astrophysical parameters from comparisons with
observed bursting events.

One of those influential reactions is the 23Al($p,\gamma$)24Si reaction.
Variations within the uncertainty of this reaction rate lead to significant
shifts in the simulated X-ray light curve. The importance of this reaction
stems from its direct connection to the 22Mg waiting point along the $rp$
process [12]. The $Q$ value for ${}^{22}{\rm Mg}(p,\gamma)$ (0.141 MeV) [13]
is small enough to establish a ($\gamma,p$)-($p,\gamma$) equilibrium [14],
which makes it difficult for the $rp$ process to proceed in this mass region.
There is evidence [14] that the 22Mg waiting point can be bypassed via the
22Mg($\alpha,p$) reaction. However, it is also possible for the ${}^{23}{\rm
Al}(p,\gamma)$ reaction to undergo resonant capture to only a few low lying
states in 24Si when certain conditions are assumed for X-ray bursts. This
network competition in the reaction flow has been investigated in previous
experiments, providing constraints on when 22Mg($\alpha,p$) becomes
significant during X-ray burst events [15]. However, uncertainties in the
23Al($p,\gamma$) reaction need to be resolved in order to confirm the flow of
the $rp$ process at this waiting point. A recent experiment measured $\gamma$
rays from ${}^{23}{\rm Al}(d,n)$ [16] to identify unbound states in 24Si,
making the uncertainty in the mass of 24Si (19 keV [13]) the dominant source
of uncertainty for the ${}^{23}{\rm Al}(p,\gamma)$ reaction rate. A high-
precision mass measurement of 24Si would constrain the $rp$ process in this
mass region by precisely quantifying the bypass of the ${}^{22}{\rm Mg}$
waiting point. In this Letter, we present a new high-precision mass
measurement of 24Si using Penning trap mass spectrometry.

A beam of 24Si was produced at the National Superconducting Cyclotron
Laboratory (NSCL) [17] via projectile fragmentation of a 28Si14+ beam, at an
energy of 160 A MeV, with a 1530 mg/cm2 thin beryllium foil target. The
resulting cocktail beam was purified using the A1900 Fragment Separator [17].
Afterwards, the beam was sent through a momentum compression beamline, passing
through aluminum degraders and an Al wedge to degrade the beam energy before
entering the gas stopper [18]. Interactions between the 24Si ions and helium
gas in the stopper slowed the rare isotope ions, which were then extracted
using a combination of direct current and radio frequency (RF) fields. After
extraction from the gas cell, the rare isotope ions were injected into a
Radio-Frequency Quadrupole (RFQ), accelerated to an energy of 30 keV and
purified with a magnetic dipole mass separator based on their mass-to-charge
ratio (A/Q). Following the beam stopping facility, 24Si was delivered as a
molecular ion, [24SiO2H]+ at $A=57$. The beam was then transported to the Low
Energy Beam and Ion Trap (LEBIT) [19] facility where the mass measurements
were performed.

The first component of the LEBIT facility is a beam cooler and buncher [20].
The ions were cooled, collected, and ejected as a pulsed beam. The pulsed ion
beam was subsequently injected into the 9.4 T Penning trap mass spectrometer
[21] and directed off-axis relative to the trap center using Lorentz steerers
[22], capturing the ions in an initial magnetron orbit. During trapping,
contaminants were purged using a sum of dipole RF excitations at the
identified contaminants’ reduced cyclotron frequencies and a broadband SWIFT
excitation [23, 24] to remove contaminant ions. The major contaminants
identified during the experiment were [C3H2F]+, [C2HS]+ and [C3H5O]+. The
Time-of-Flight Ion Cyclotron Resonance (TOF-ICR) technique was then used to
measure the mass of the ion of interest. In this method, a quadrupole
excitation is performed for a range of frequencies near the cyclotron
frequency, $\nu_{c}=qB/2\pi m$, of the ion of interest [25, 26]. As the RF
excitation, $\nu_{RF}$, approaches $\nu_{c}$, the magnetron motion is
converted to reduced cyclotron motion. The conversion is maximized when
$\nu_{RF}\approx\nu_{c}$. Following the excitation, the ions are ejected from
the trap to a microchannel plate (MCP) detector, where the time of flight is
recorded. During the conversion, the ions gain energy, which is observed as a
minimum in the TOF response curve. The theoretical line shape can be fit to
the data to extract $\nu_{c}$ [26].

Figure 1: An example of a 75-ms [24SiO2H]+ time-of-flight Ramsey cyclotron
resonance curve from which $\nu_{c}$ was determined. The line shape is
described in Ref. [27].

Another excitation scheme, known as a Ramsey excitation [27], was also
performed in this experiment. This generates a TOF pattern visible in FIG. 1
and described by Ref [27], which allows superior precision to the simple RF
excitation scheme for the same ion trapping times. 18 measurements of
[24SiO2H]+ were performed, consisting of 6 normal quadrupole excitation
measurements and 12 Ramsey excitation measurements. The normal quadrupole
excitation measurements were performed with 50 ms and 75 ms excitations. The
Ramsey measurements were performed with 25 ms, 50 ms, 75 ms, 100 ms and 150 ms
excitations.

Cyclotron frequency measurements of [24SiO2H]+ were bracketed between two
reference measurements of [12C3H2F]+ to track variations in the magnetic field
strength. The weighted average of the ratio for
$\nu_{c}$([24SiO2H]+)/$\nu_{c}$([12C3H2F]+) was determined to be
$\bar{R}=1.000085151(71)$ with an associated Birge Ratio [28] of 1.03. Most
systematic uncertainties in $\bar{R}$ scale linearly with the mass difference
between the ion of interest and the reference ion [25]. This mass-dependent
shift has been measured to be $\Delta\bar{R}=2\times 10^{-10}$/u [29], which
is negligible for isobaric species. Other systematic uncertainties include
relativistic effects, ion-ion interactions, and non-linear magnetic field
fluctuations. Relativistic effects are negligible compared to the statistical
uncertainty [30]. Ion-ion interactions were limited by excluding shots with
more than 5 ions. The non-linear fluctuations in the magnetic field are on the
order of $\Delta\bar{R}=10^{-9}$/hour [31]. The longest single measurement was
2.5 hours making these fluctuations negligible as well. As a result, the
uncertainty of our averaged frequency measurement is statistical.

Figure 2: The reaction rates for 23Al($p,\gamma$) (Wolf et al. 2019 [16],
denoted by the superscript 1, and this work, denoted by the superscript 2),
22Mg($p,\gamma$) (Iliadis et al. 2010, [32]) and 22Mg($\alpha,p$) (Randhawa et
al. 2020, [15]) from 0.1 to 2 GK. Our calculations show a dramatic reduction
in the 1-$\sigma$ reaction rate uncertainty for 23Al($p,\gamma$) as a result
of our new mass measurement and the corresponding resonance energies.

The mass excess of 24Si is extracted directly from $\bar{R}$, and was
determined to be $10\;753.8$(37) keV. Our mass is in good agreement with the
AME2020 (10 745(19) keV), and a factor of 5 more precise. The AME2020 value is
primarily based on a 1980 pion double-charge exchange reaction measurement
[33]. With our measurement, a new proton separation energy for 24Si was
calculated to be $3\;283.2$(37) keV.

With a precisely known $Q$ value, the flow following the 22Mg waiting point in
the $rp$ process is strictly constrained based on nuclear physics information.
FIG. 2 displays the impact of the new $Q$ value on the uncertainty of the
23Al($p,\gamma$) reaction rate. Also shown in FIG. 2 are the 22Mg($p,\gamma$)
reaction rate, which is represented by the black band, based on information
from Ref. [32], and the reaction rate experimentally determined for the
22Mg($\alpha,p$) reaction represented by the blue band from Ref. [15]. We
estimate the ${}^{22}{\rm Mg}(\alpha,p)$ reaction rate uncertainty to be
$\pm$80% based on the uncertainty for the lowest-energy cross section
measurement from Ref. [15].

Table 1: Resonance levels used in the narrow-resonance approximation for the reaction rate calculation. Experimental information is taken from Ref. [16]. Other levels were calculated with the NuShellX@MSU package with their corresponding spectroscopic factors, radiative widths and proton widths [34], as denoted by *. Excitation energies from NuShellX@MSU calculations were defined with uncertainties of 150 keV [35]. The reaction rate was calculated with a Monte Carlo program to determine a robust rate uncertainty. A complete description of the reaction rate calculation can be found in the text. $E_{x}$ (MeV) | $J$ | $l$ | $\text{C}^{2}\text{S}$ |  $\Gamma_{\gamma}$ (eV) |  $\Gamma_{p}$ (eV)
---|---|---|---|---|---
3.449(5) | 2 | 0 | 0.7(4) | $1.9\times 10^{-2}$ | $1.0\times 10^{-4}$
|  | 2 | 0.002(1) |  |
|  | 2 | 0.3(2) |  |
3.471(6) | 0 | 2 | 0.8(4) | $1.6\times 10^{-3}$ | $6.2\times 10^{-5}$
4.256(150)* | 3 | 0 | 0.59 | $1.3\times 10^{-2}$ | $9.0\times 10^{3}$
|  | 2 | 0.17 |  |
5.353(150)* | 3 | 0 | 0.0012 | $2.8\times 10^{-2}$ | $3.6\times 10^{3}$
|  | 2 | 0.11 |  |
5.504(150)* | 2 | 0 | 0.044 | $2.2\times 10^{-1}$ | $2.8\times 10^{4}$
|  | 2 | 0.068 |  |
5.564(150)* | 4 | 2 | 0.048 | $2.2\times 10^{-2}$ | $2.2\times 10^{3}$
6.004(150)* | 4 | 2 | 0.28 | $6.5\times 10^{-3}$ | $2.8\times 10^{4}$
6.056(150)* | 0 | 2 | 0.053 | $3.4\times 10^{-3}$ | $5.6\times 10^{3}$
6.072(150)* | 2 | 0 | 0.012 | $5.0\times 10^{-2}$ | $2.4\times 10^{4}$
|  | 2 | 0.093 |  |

The ${}^{23}{\rm Al}$($p,\gamma$) reaction rate was calculated using the
narrow resonance approximation [36]. The information that contributed to this
calculation can be found in Table 1, based on experimental information [16],
and, for higher lying states, on calculations with NuShellX@MSU [34] that were
included in the reaction rate calculation of Ref. [16]. Based on a previous
measurement of longitudinal momentum distributions of ${}^{24}{\rm Si}$ states
[37], we assigned the 3471 state with J=0 and used the information from Ref.
[16].

The reaction rate properties varied using a Monte Carlo random sampling
approach. The excitation energies in the unbound states of ${}^{24}{\rm Si}$
and the $Q$ value for the ${}^{23}{\rm Al}(p,\gamma)$ reaction were varied,
assuming a normal distribution. The impact of the proton penetration factor
was also included in the reaction rate calculation [38]. The excitation
energies derived from NuShellX@MSU calculations were defined with
uncertainties of 150 keV [35]. The radiative widths and spectroscopic factors
were varied with a Log-normal distribution [39]. All radiative partial widths
assumed a factor uncertainty of 2 [40]. The same factor uncertainty was used
for all theoretical spectroscopic factors. Experimental spectroscopic factors
and uncertainties were varied using relations found in Ref. [39]. The proton
partial widths were scaled, based on the impact of the spectroscopic factor
variation.

For the direct capture component of the reaction rate calculation, the
spectroscopic factors for each astrophysical S-factor reported in Ref. [41]
were varied with a Log-normal distribution as well. The experimental ground
state spectroscopic factor was varied with its experimental uncertainty while
the three theoretical values calculated for the first excited state assumed a
factor uncertainty of 2. After varying the spectroscopic factor, the
corresponding astrophysical S-factors $S(E_{0})$ were scaled to calculate the
new S-factor before summing all four contributions in the direct capture
component calculation. The results for the calculated rate for the
${}^{23}{\rm Al}$($p,\gamma$) rate can be found in Table 2.

Table 2: The recommended reaction rate based on information from Ref. [16] and this work. Also included are the lower (16th percentile) and upper (64th percentile) uncertainties of this reaction. $T_{9}$ |  |  $N_{A}<\sigma v>$ ($\text{cm}^{3}/\text{s}/\text{mole}$) |
---|---|---|---
| Recommended | Lower | Upper
0.1 | $8.195\times 10^{-7}$ | $4.777\times 10^{-7}$ | $1.402\times 10^{-6}$
0.15 | $2.754\times 10^{-4}$ | $1.640\times 10^{-4}$ | $4.642\times 10^{-4}$
0.2 | $4.481\times 10^{-3}$ | $2.602\times 10^{-3}$ | $7.761\times 10^{-3}$
0.3 | $6.168\times 10^{-2}$ | $3.432\times 10^{-2}$ | $1.114\times 10^{-1}$
0.4 | $2.022\times 10^{-1}$ | $1.099\times 10^{-1}$ | $3.746\times 10^{-1}$
0.5 | $3.834\times 10^{-1}$ | $2.052\times 10^{-1}$ | $7.212\times 10^{-1}$
0.6 | $5.594\times 10^{-1}$ | $2.967\times 10^{-1}$ | 1.064
0.7 | $7.096\times 10^{-1}$ | $3.745\times 10^{-1}$ | 1.360
0.8 | $8.328\times 10^{-1}$ | $4.392\times 10^{-1}$ | 1.601
0.9 | $9.381\times 10^{-1}$ | $4.950\times 10^{-1}$ | 1.801
1.0 | 1.034 | $5.482\times 10^{-1}$ | 1.979
1.5 | 1.707 | $9.329\times 10^{-1}$ | 3.303
2.0 | 3.222 | 2.010 | 6.553
Figure 3: Flow through 22Mg($\alpha,p$) for the temperature range of relevance
for X-ray bursts, where the solid black lines assume the ${}^{23}{\rm
Al}(p,\gamma)$ reaction rate uncertainty of Ref. [16] and the dashed red lines
correspond to the uncertainty from this work. The inset shows the uncertainty
spread in the ($\alpha,p$) flow, which peaks at temperatures reached during
the light curve rise in X-ray burst model calculations.

We assessed the impact of the new ${}^{24}{\rm Si}$ mass on the onset of the
$\alpha p$-process by calculating the flow at the ${}^{22}{\rm Mg}$
bottleneck. To first order, at this point in the reaction network, there is a
competition between ${}^{23}{\rm Al}(p,\gamma)$ and ${}^{22}{\rm
Mg}(\alpha,p)$, given the rather long half-life of ${}^{22}{\rm Mg}$, 3.88 s
[42], and low proton separation energy $S_{p}$ of ${}^{23}{\rm Al}$, 0.141 MeV
[13]. The flow into the $\alpha p$-process is defined as
$\lambda_{\alpha,p}/\left(\lambda_{p,\gamma}+\lambda_{\alpha,p}\right)$. The
rates are $\lambda_{i}=W\left(X_{i}/A_{i}\right)\rho N_{\rm A}\langle\sigma
v\rangle$, where $X_{i}$ and $A_{i}$ are the mass fraction and mass number of
hydrogen and helium for the ${}^{23}{\rm Al}(p,\gamma)$ and ${}^{22}{\rm
Mg}(\alpha,p)$ reactions, respectively [43]. $W$ is a weight factor, which is
the equilibrium abundance of the nuclide as determined by the Saha equation
[6, 44]. Since the ${}^{22}{\rm Mg}$ $\beta^{+}$-decay rate is too slow to
compete, only strong and electromagnetic interactions are involved in the
flow, and the density dependence cancels. For $X_{i}$, we take the conditions
determined by Ref. [45] for the X-ray burst ignition region, where $X_{\rm
H}=0.03$ and $X_{\rm He}=0.31$.

Results are shown in FIG. 3 over the range of temperatures of relevance for
X-ray bursts. For the following discussion, we adopt 10% of flow through
${}^{22}{\rm Mg}(\alpha,p)$ as the point when the $\alpha p$-process becomes
significant. Previously, using the ${}^{23}{\rm Al}(p,\gamma)$ reaction rate
uncertainty from Ref. [16], the uncertainty band for when the $\alpha
p$-process becomes significant is $T_{9}=0.94\pm 0.15$ GK. With the new
${}^{24}{\rm Si}$ mass, this uncertainty band becomes $0.93\pm 0.08$ GK,
precisely pinpointing the onset of the $\alpha p$-process in X-ray bursts.
While adopting other ignition conditions shifts the temperature for
significant $\alpha p$-process flow, e.g. $0.99\pm 0.09$ GK for the $X_{\rm
H}=0.06$ $X_{\rm He}=0.19$ ignition conditions from the calculations of Ref.
[43], the precision remains 9%. Therefore, we firmly establish the relevance
of the $\alpha p$-process in powering X-ray burst light curves by determining
under which conditions there is $(\alpha,p)$ breakout from the ${}^{22}{\rm
Mg}$ waiting-point.

This Letter presents the first high-precision mass measurement of ${}^{24}{\rm
Si}$ using Penning trap mass spectrometry, resulting in a mass excess of
$10\;753.8$(37) keV. The achieved precision reduces the uncertainty of the
${}^{23}{\rm Al}(p,\gamma)$ reaction rate, pinpointing the temperature when
the the $\alpha p$-process turns-on in X-ray bursts with ${}^{22}{\rm
Mg}(\alpha,p)$ to 9%. Further refinement would require more precise resonance
strength determinations for ${}^{23}{\rm Al}(p,\gamma)$, as was recently
submitted for publication [46], and a measurement of the ${}^{22}{\rm
Mg}(\alpha,p)$ cross section to lower center-of-mass energies.

###### Acknowledgements.

This material is based upon work supported by the U.S. National Science
Foundation through Grants No. PHY-1565546, PHY-2111185 (The Precision Frontier
at FRIB: Masses, Radii, Moments, and Fundamental Interactions), PHY-1430152
(Joint Institute for Nuclear Astrophysics – Center for the Evolution of the
Elements), and No. PHY-191355, as well as the U.S. Department of Energy,
Office of Science, Office of Nuclear Physics under Grants No. DE-
FG02-88ER40387, DE-SC0019042, DE-NA0003909, and DE-SC0015927. A.A.V.
acknowledges support from NSERC (Canada) under Contract No. SAPPJ2018-00028.
Part of this work was performed under the auspices of the U.S. Department of
Energy by Lawrence Livermore National Laboratory under Contract DE-
AC52-07NA27344.

## References

  * Lewin  _et al._ [1993] W. H. Lewin, J. Van Paradijs, and R. E. Taam, Space Sci. Rev. 62, 223 (1993).
  * Schatz and Rehm [2006] H. Schatz and K. Rehm, Nucl. Phys. A 777, 601 (2006).
  * Schatz  _et al._ [1998] H. Schatz, A. Aprahamian, J. Görres, M. Wiescher, T. Rauscher, J. Rembges, F.-K. Thielemann, B. Pfeiffer, P. Möller, K.-L. Kratz, H. Herndl, B. Brown, and H. Rebel, Phys. Rep. 294, 167 (1998).
  * Wallace and Woosley [1981] R. Wallace and S. E. Woosley, Astrophys. J. Suppl. Ser. 45, 389 (1981).
  * Bult  _et al._ [2021] P. Bult, D. Altamirano, Z. Arzoumanian, A. V. Bilous, D. Chakrabarty, K. C. Gendreau, T. Güver, G. K. Jaisawal, E. Kuulkers, C. Malacaria, M. Ng, A. Sanna, and T. E. Strohmayer, Astrophys. J. 907, 79 (2021).
  * Meisel  _et al._ [2018] Z. Meisel, A. Deibel, L. Keek, P. Shternin, and J. Elfritz, J. Phys. G 45, 093001 (2018).
  * Schatz  _et al._ [2001] H. Schatz, A. Aprahamian, V. Barnard, L. Bildsten, A. Cumming, M. Ouellette, T. Rauscher, F.-K. Thielemann, and M. Wiescher, Phys. Rev. Lett. 86, 3471 (2001).
  * Heger  _et al._ [2007] A. Heger, A. Cumming, D. K. Galloway, and S. E. Woosley, Astrophys. J Lett. 671, L141 (2007).
  * Zamfir  _et al._ [2012] M. Zamfir, A. Cumming, and D. K. Galloway, Astrophys. J. 749, 69 (2012).
  * Cyburt  _et al._ [2016] R. H. Cyburt, A. M. Amthor, A. Heger, E. Johnson, L. Keek, Z. Meisel, H. Schatz, and K. Smith, Astrophys. J. 830, 55 (2016).
  * Parikh  _et al._ [2008] A. Parikh, J. José, F. Moreno, and C. Iliadis, Astrophys. J. Suppl. Ser. 178, 110 (2008).
  * Fisker  _et al._ [2004] J. L. Fisker, F.-K. Thielemann, and M. Wiescher, Astrophys. J. 608, L61 (2004).
  * Huang  _et al._ [2021] W. Huang, M. Wang, F. Kondev, G. Audi, and S. Naimi, Chin. Phys. C 45, 030002 (2021).
  * van Wormer  _et al._ [1994] L. van Wormer, J. Görres, C. Iliadis, M. Wiescher, and F. K. Thielemann, Astrophys. J. 432, 326 (1994).
  * Randhawa  _et al._ [2020] J. S. Randhawa, Y. Ayyad, W. Mittig, Z. Meisel, T. Ahn, S. Aguilar, H. Alvarez-Pol, D. W. Bardayan, D. Bazin, S. Beceiro-Novo, D. Blankstein, L. Carpenter, M. Cortesi, D. Cortina-Gil, P. Gastis, M. Hall, S. Henderson, J. J. Kolata, T. Mijatovic, F. Ndayisabye, P. O’Malley, J. Pereira, A. Pierre, H. Robert, C. Santamaria, H. Schatz, J. Smith, N. Watwood, and J. C. Zamora, Phys. Rev. Lett. 125, 202701 (2020).
  * Wolf  _et al._ [2019] C. Wolf, C. Langer, F. Montes, J. Pereira, W.-J. Ong, T. Poxon-Pearson, S. Ahn, S. Ayoub, T. Baumann, D. Bazin, P. C. Bender, B. A. Brown, J. Browne, H. Crawford, R. H. Cyburt, E. Deleeuw, B. Elman, S. Fiebiger, A. Gade, P. Gastis, S. Lipschutz, B. Longfellow, Z. Meisel, F. M. Nunes, G. Perdikakis, R. Reifarth, W. A. Richter, H. Schatz, K. Schmidt, J. Schmitt, C. Sullivan, R. Titus, D. Weisshaar, P. J. Woods, J. C. Zamora, and R. G. T. Zegers, Phys. Rev. Lett. 122, 232701 (2019).
  * Morrissey  _et al._ [2003] D. Morrissey, B. Sherrill, M. Steiner, A. Stolz, and I. Wiedenhoever, Nucl. Instrum. Methods Phys. Res. B 204, 90 (2003).
  * Sumithrarachchi  _et al._ [2020] C. Sumithrarachchi, D. Morrissey, S. Schwarz, K. Lund, G. Bollen, R. Ringle, G. Savard, and A. Villari, Nucl. Instrum. Methods Phys. Res. B 463, 305 (2020).
  * Ringle  _et al._ [2013] R. Ringle, S. Schwarz, and G. Bollen, Int. J. Mass Spectrom. 349-350, 87 (2013).
  * Schwarz  _et al._ [2016] S. Schwarz, G. Bollen, R. Ringle, J. Savory, and P. Schury, Nucl. Instrum. Methods Phys. Res. A 816, 131 (2016).
  * Ringle  _et al._ [2009] R. Ringle, G. Bollen, A. Prinke, J. Savory, P. Schury, S. Schwarz, and T. Sun, Nucl. Instrum. Methods Phys. Res. A 604, 536 (2009).
  * Ringle  _et al._ [2007a] R. Ringle, G. Bollen, A. Prinke, J. Savory, P. Schury, S. Schwarz, and T. Sun, Int. J. Mass Spectrom. 263, 38 (2007a).
  * Blaum  _et al._ [2004] K. Blaum, D. Beck, G. Bollen, P. Delahaye, C. Guénaut, F. Herfurth, A. Kellerbauer, H.-J. Kluge, D. Lunney, S. Schwarz, L. Schweikhard, and C. Yazidjian, Europhys. Lett. (EPL) 67, 586 (2004).
  * Kwiatkowski  _et al._ [2015] A. Kwiatkowski, G. Bollen, M. Redshaw, R. Ringle, and S. Schwarz, Int. J. Mass Spectrom. 379, 9 (2015).
  * Bollen  _et al._ [1990] G. Bollen, R. B. Moore, G. Savard, and H. Stolzenberg, J. Appl. Phys. 68, 4355 (1990).
  * König  _et al._ [1995] M. König, G. Bollen, H.-J. Kluge, T. Otto, and J. Szerypo, Int. J. Mass Spectrom. Ion Processes 142, 95 (1995).
  * Kretzschmar [2007] M. Kretzschmar, Int. J. Mass Spectrom. 264, 122 (2007).
  * Birge [1932] R. T. Birge, Phys. Rev. 40, 207 (1932).
  * Gulyuz  _et al._ [2015] K. Gulyuz, J. Ariche, G. Bollen, S. Bustabad, M. Eibach, C. Izzo, S. J. Novario, M. Redshaw, R. Ringle, R. Sandler, S. Schwarz, and A. A. Valverde, Phys. Rev. C 91, 055501 (2015).
  * Brown and Gabrielse [1986] L. S. Brown and G. Gabrielse, Rev. Mod. Phys. 58, 233 (1986).
  * Ringle  _et al._ [2007b] R. Ringle, T. Sun, G. Bollen, D. Davies, M. Facina, J. Huikari, E. Kwan, D. J. Morrissey, A. Prinke, J. Savory, P. Schury, S. Schwarz, and C. S. Sumithrarachchi, Phys. Rev. C 75, 055503 (2007b).
  * Iliadis  _et al._ [2010] C. Iliadis, R. Longland, A. Champagne, A. Coc, and R. Fitzgerald, Nucl. Phys. A 841, 31 (2010).
  * Burleson  _et al._ [1980] G. R. Burleson, G. S. Blanpied, G. H. Daw, A. J. Viescas, C. L. Morris, H. A. Thiessen, S. J. Greene, W. J. Braithwaite, W. B. Cottingame, D. B. Holtkamp, I. B. Moore, and C. F. Moore, Phys. Rev. C 22, 1180 (1980).
  * Brown and Rae [2014] B. A. Brown and W. D. M. Rae, Nucl. Data Sheets 120, 115 (2014).
  * Henderson and Stroberg [2020] J. Henderson and S. R. Stroberg, Phys. Rev. C 102, 031303 (2020).
  * Fowler  _et al._ [1967] W. A. Fowler, G. R. Caughlan, and B. A. Zimmerman, Annu. Rev. Astron. Astrophys. 5, 525 (1967).
  * Longfellow  _et al._ [2020] B. Longfellow, A. Gade, J. A. Tostevin, E. C. Simpson, B. A. Brown, A. Magilligan, D. Bazin, P. C. Bender, M. Bowry, B. Elman, E. Lunderberg, D. Rhodes, M. Spieker, D. Weisshaar, and S. J. Williams, Phys. Rev. C 101, 031303 (2020).
  * Humblet  _et al._ [1987] J. Humblet, W. A. Fowler, and B. A. Zimmerman, Astron. Astrophys. 177, 317 (1987).
  * Longland  _et al._ [2010] R. Longland, C. Iliadis, A. Champagne, J. Newton, C. Ugalde, A. Coc, and R. Fitzgerald, Nucl. Phys. A 841, 1 (2010).
  * Richter  _et al._ [2020] W. A. Richter, B. A. Brown, R. Longland, C. Wrede, P. Denissenkov, C. Fry, F. Herwig, D. Kurtulgil, M. Pignatari, and R. Reifarth, Phys. Rev. C 102, 025801 (2020).
  * Banu  _et al._ [2012] A. Banu, F. Carstoiu, N. L. Achouri, W. N. Catford, M. Chartier, B. Fernández-Domínguez, M. Horoi, B. Laurent, N. A. Orr, S. Paschalis, N. Patterson, B. Pietras, B. T. Roeder, P. Roussel-Chomaz, J. S. Thomas, L. Trache, and R. E. Tribble, Phys. Rev. C 86, 015806 (2012).
  * Basunia [2015] M. S. Basunia, Nuclear Data Sheets 127, 69 (2015).
  * Merz and Meisel [2021] G. Merz and Z. Meisel, Mon. Not. R. Astron. Soc. 500, 2958 (2021).
  * Brown [2015] E. Brown,  _Stellar Astrophysics_ (Open Astrophysics Bookshelf, 2015).
  * Fisker  _et al._ [2008] J. L. Fisker, H. Schatz, and F.-K. Thielemann, Astrophys. J. Suppl. Ser. 174, 261 (2008).
  * [46] G. Lotay, Private communication.


††thanks:

]Received - -

# The Bardeen-Petterson Effect as an Observable Support for the Blanford-Payne
Process Black Hole Jet Production

Uicheol Jang Jang<EMAIL_ADDRESS>Department of Astronomy, Space Science and
Geology, Chungnam National Univ., Daejeon 34134  Hongsu Kim  Center for
Theoretical Astronomy , Korea Astronomy and Space Science Institute, Daejeon
34055

([)

###### Abstract

Relativistic jets are observed around accreting black holes, from stellar mass
to supermassive black holes. However, their origin has not been fully
understood. Although the Blanford-Payne process has been considered the most
reliable theory for jet production, it has no observational proof or support.
In this study, we suggest that the Bardeen-Petterson effect can provide
observational evidence for the Blanford-Payne process. We studied black hole
jet production using the Blanford-Payne process and the Bardeen-Peterson
effect. As a result, we were able to calculate the observable timescale of
black hole jet precession.

black hole jet, accretion disk

##  I Introduction

Accreting black holes are observed to produce relativistic jets. Although the
production of active galactic nuclei (AGN) jets has been known as the
Blanford-Payne processblandford1982hydromagnetic , observational evidence or
support has been rare. Therefore, we would like to suggest observable cases
that can support observational evidence of the Blanford-Payne process. The
Blanford-Payne process, as described, is a mechanism for producing AGN jets
that are driven by gas (corona) pressure blandford1982hydromagnetic . This
interpretation shows that the jet outflow is closely related to the accretion
disk. And the accretion disk around a spinning black hole is affected by frame
dragging, namely, the Einstein-Lense-Thirring effect lense1918einfluss . This
is well known as the Bardeen-Peterson effectbardeen1975lense . We suggest the
coupling of the Blanford-Payne process to the Bardeen-Peterson effect. The
concurrence of these two theories implies that the Blanford-Payne process can
be distinguished as a theory through the Bardeen-Peterson effect.

##  II Production of AGN Jet outflow.

###  II.1 Blandford-Payne process

As we remarked in the introduction, the Blanford-Payne process explains the
production of AGN jet outflows by employing magnetohydrodynamics (MHD) in
1982. According to the theory, a centrifugally driven outflow of matter from
the disk is possible if the poloidal component of the magnetic field makes an
angle of less than 60° with the disk surface. At large distances from the
disk, the toroidal component of the magnetic field becomes important and
collimates the outflow into a pair of antiparallel jets moving perpendicular
to the disk. In this way, magnetic stress can extract energy and angular
momentum from an accretion disk independently of the presence of
viscosity.blandford1982hydromagnetic
.https://doi.org/10.48550/arxiv.1508.05192 . Moreover General relativistic
magnetohydrodynamics(GRMHD) simulations have shifted the consensus towards the
understanding that the Blandford&Payne mechanism launches sub-relativistic
winds that collimate a relativistic jet powered by the BZ-
mechanism.koide_shibata_kudoh_1998 ,mizuno_2022 10.1007/978-94-011-4780-4_66
ART000959830

###  II.2 Timescale for Accretion jet flow.

As a simple estimate, the time scale of the whole outflow process, $\Delta
t_{BP}$, would be given by;

$\Delta t_{BP}\sim\frac{l}{c_{v}}$ (1)

Where we denote that $c_{v}$ is the Alfvén velocity, and $l$ is the length
scale over which the poloidal magnetic field extends, which are given as:

$c_{v}=(\frac{dp}{d\rho})^{\frac{1}{2}}=\frac{B}{\sqrt{\mu_{0}\rho}},$ (2)

$l=10^{3}\frac{GM}{c^{2}}$ (3)

Respectively, for the radio jet (or accretion disk) length scale of AGN, where
we denote that $\rho$ is the plasma density, $B$ is the magnetic flux, $p$ is
the gas pressure of plasma, $\mu_{0}$ is the vacuum permeability, $G$ is
Newton’s gravitational constant, $c$ is the speed of light, and $M$ is the
mass of the black hole.

##  III Precession of the accretion disk.

###  III.1 The Bardeen-Petterson effect

The deformation of the accretion disk has been demonstrated by General
relativistic magnetohydrodynamics(GRMHD)RN42 . the Bardeen-Petterson effect
presented the influence on tilted accretion disks around Kerr black holes as
astrophysical evidence for frame dragging (Lense-Thirring effect)
bardeen1975lense .This effect occurs within the Bardeen-Petterson radius,
$R_{BP}$, which is given, roughly, by;

$R_{BP}\gtrsim 10^{2}\frac{GM}{c^{2}}$ (4)

which reaches far outside the horizon.

###  III.2 Time scale for accretion disk precession

Since the total torque involved in the Bardeen-Petterson effect consists of
three parts: alignment, precession, and spin-down components, the precession
time scale and the alignment time scale would be of the same order;

$\Delta t_{LT}\sim\frac{2\pi}{\Omega_{GM}}$ (5)

Where $\Omega_{GM}$ is the precession angular velocity, which due to the
Lense-Thirring effect is given by Wilkins.wilkins1972bound

$\Omega_{GM}=\frac{2GJ}{c^{2}R^{3}_{BP}}\lesssim
a_{*}\frac{2G^{2}M^{2}}{c^{3}(100GM/c^{2})^{3}}=10^{-6}\frac{2c^{3}a_{*}}{GM};$
(6)

$J$ is the total angular momentum of the black hole, which is given by
$J=2GM^{2}a_{/}c$, and $a_{*}$ is the angular momentum ratio with $0\leq
a_{*}<1$. narayan2012observational .

##  IV Luminosity

Luminous AGN jets are some of the most energetic and powerful phenomena in the
universe, and they can produce high-energy radiation and particles that can be
detected across the electromagnetic spectrum, from radio waves to gamma rays.
AGN jets are classified into two main types based on their luminosity:
”quasar-like” jets, which are more luminous, and ”radio galaxy” jets, which
are less luminous. Quasar-like jets are typically associated with more massive
black holes and are more common in the early universe, while radio galaxy jets
are typically associated with less massive black holes and are more common in
the present-day universe. AGN jets are thought to be powered by the accretion
of matter onto a supermassive black hole at the center of a galaxy. As
material falls onto the black hole, it is heated to extremely high
temperatures and emits high-energy radiation and particles. Some of this
material is then accelerated to near the speed of light and ejected from the
vicinity of the black hole in the form of a jet. The particles in the jet are
primarily electrons, protons, and other charged particles, and they can be
accelerated to energies of up to $10^{20}eV$ or moreblandford1982hydromagnetic
. The jet is also accompanied by a magnetic field, which is thought to play a
key role in the acceleration and emission of particles from the jet.

Black hole jets are observed at different frequencies, depending on the
specific study and type of jet being observed. The Atacama Large Millimeter
Array (ALMA) is an interferometric array of radio telescopes that operate at
millimeter and submillimeter wavelengths, corresponding to frequencies that
are much higher than those observed by traditional radio telescopes. These
frequencies typically range from 30 to 950 GHz. This high-frequency range
allows ALMA to make very high-resolution observations of the jets, providing
detailed information about their structure and dynamics. For example, ALMA
used a 1.3mm (Band 6) wavelength to observe the jet of M87, as shown in (FIG
1) of the study.Goddi_2021  It is thought that jets from other black holes
could be observed using similar methods.

##  V Concurrence of the Blanford-Payne process with the Bardeen-Peterson
effect.

Close to the disk, the plasma flow is driven by gas pressure in a hot,
magnetically dominated corona blandford1982hydromagnetic . This produces the
jet outflow. The disk should be shifted periodically by frame dragging
bardeen1975lense . Therefore, when the Blanford-Payne process and the Bardeen-
Peterson effect work simultaneously, the AGN jet outflow is likely to rotate
along the time scale ($T_{LT}$).

We begin with a necessary condition for the process to occur. To make sure
that the events can be well resolved in a time sequence, we demand the
following condition in this context:

If $\Delta t_{BP}\lesssim\Delta t_{LT}$, the concurrence of the Blandford-
Payne process and the Bardeen-Petterson effect can be resolved (Figure 1).
This result can be applied to the supermassive black hole (SMBH) the M87,
which has been observed directly by the Event Horizon Telescope (EHT).
event2019first algaba2021broadband . SMBH M87 given $M\sim(6.5\pm 0.7)\times
10^{9}M_{\odot}$walsh2013m87  and $B=1\sim 34G$ akiyama2021first
akiyama2021first2 . The calculated time scales are $\Delta t_{BP}\sim 1.015$
years and $\Delta t_{LT}\sim 3986.7$ years, which is an impossible value for
observation. This is why no one has detected the jet precessions of M87. These
results imply that the observable properties can be calculated by selecting
smaller black holes. In the same manner, Sagittarius $A^{*}$ (which was
directly observed by EHT) has a mass of $M\sim 4\times 10^{6}M_{\odot}$
akiyama2022first , The time scales for this black hole are $\Delta t_{BP}\sim
35.6$ hours and $\Delta t_{LT}\sim 15.9$ years.

We acknowledge that the jet of Sagittarius A* is generally believed to have a
relatively low observable power compared to other known black hole jets. But
Regarding Sgr A*, recent observations have revealed that it undergoes
occasional flares in the submillimeter and infrared
bandsmurchikova_witzel_2021 , indicating the presence of a compact and
variable emission region close to the black hole. While no jet has been
directly observed from Sgr A*, it is still a viable candidate for exhibiting
jet precession, as the same physical processes that drive jet precession in
other systems could also operate in Sgr A*. Furthermore, recent theoretical
models have suggested that a Blandford-Znajek powered jet may be present in
Sgr A*, which could be relevant to our study.

In our paper, we discuss the timescale for jet precession in the context of
accreting black holes, regardless of whether or not they exhibit an observable
jet. Our analysis combines well-established theoretical frameworks to provide
a quantitative estimate of this timescale, which could help guide future
observational campaigns. We acknowledge that there are uncertainties in the
existence of Blandford-Payne winds and Bardeen-Petterson alignment, as well as
in the applicability of our study to specific objects. However, we believe
that our study provides a valuable contribution to the theoretical
understanding of jet precession in accreting systems, and we hope that it will
stimulate further discussions and investigations in this area.

If the mass of the black hole is smaller, the time scales will be more
reliable values to observe. Therefore, we suggest time series observations of
black hole jets. If the time series continuity of the concurrence of these two
events is observed, it will have a shape (FIG 2).

##  VI Discussion

We are discussing the relationship between the Blanford-Payne process and the
Bardeen-Peterson effect in the context of AGN jet outflow. We note that when
these two processes occur simultaneously, the jet outflow is expected to
rotate on a time scale referred to as $T_{LT}$. In order to observe this
rotation, We impose a necessary condition that $\Delta t_{BP}\lesssim\Delta
t_{LT}$.

We apply this condition to two supermassive black holes, M87 and Sagittarius
$A^{*}$, which have been observed by the Event Horizon Telescope. We find that
the calculated time scales for M87 make it impossible to observe the
precession of its jet, but the time scales for Sagittarius $A^{*}$ are more
feasible for observation. Based on this, the authors suggest that future
observations of black hole jets be made on smaller black holes, where the time
scales are more reliable.

Our overall aim is to provide robust support for the Blanford-Payne mechanism
and its effectiveness in real-world scenarios. We argue that the inclusion of
the Bardeen-Peterson effect is critical in their assessment of the viability
of the Blanford-Payne mechanism, as it enables them to accurately evaluate its
practical relevance.

In conclusion, We provides a theoretical analysis of the relationship between
the Blanford-Payne process and the Bardeen-Peterson effect, and its potential
impact on the observation of AGN jet outflows. We suggest that future
observations be made on smaller black holes to obtain reliable results, and
argue that the inclusion of the Bardeen-Peterson effect is necessary to assess
the viability of the Blanford-Payne mechanism.

## Acknowledgement

This research was supported by Basic Science Research Program through the
National Research Foundation of Korea(NRF) funded by the Ministry of Education
(NRF-2022R1A2C1092602).

## References

  * [1] Kazunori Akiyama, Antxon Alberdi, Walter Alef, Juan Carlos Algaba, Richard Anantua, Keiichi Asada, Rebecca Azulay, Uwe Bach, Anne-Kathrin Baczko, David Ball, et al.  First sagittarius a* event horizon telescope results. ii. eht and multiwavelength observations, data processing, and calibration.  The Astrophysical Journal Letters, 930(2):L13, 2022.
  * [2] Kazunori Akiyama, Juan Carlos Algaba, Antxon Alberdi, Walter Alef, Richard Anantua, Keiichi Asada, Rebecca Azulay, Anne-Kathrin Baczko, David Ball, Mislav Baloković, et al.  First m87 event horizon telescope results. vii. polarization of the ring.  The Astrophysical Journal Letters, 910(1):L12, 2021.
  * [3] Kazunori Akiyama, Juan Carlos Algaba, Antxon Alberdi, Walter Alef, Richard Anantua, Keiichi Asada, Rebecca Azulay, Anne-Kathrin Baczko, David Ball, Mislav Baloković, et al.  First m87 event horizon telescope results. viii. magnetic field structure near the event horizon.  The Astrophysical Journal Letters, 910(1):L13, 2021.
  * [4] JC Algaba, J Anczarski, K Asada, M Baloković, S Chandra, Y-Z Cui, AD Falcone, M Giroletti, C Goddi, K Hada, et al.  Broadband multi-wavelength properties of m87 during the 2017 event horizon telescope campaign.  The Astrophysical journal letters, 911(1):L11, 2021.
  * [5] James M Bardeen and Jacobus A Petterson.  The lense-thirring effect and accretion disks around kerr black holes.  The Astrophysical Journal, 195:L65, 1975.
  * [6] Roger D Blandford and DG Payne.  Hydromagnetic flows from accretion discs and the production of radio jets.  Monthly Notices of the Royal Astronomical Society, 199(4):883–903, 1982.
  * [7] Anderson Caproni, Zulema Abraham, and Herman J. Mosquera Cuesta.  Bardeen‐petterson effect and the disk structure of the seyfert galaxy ngc 1068.  The Astrophysical Journal, 638(1):120–124, 2006.
  * [8] Event Horizon Telescope Collaboration et al.  First m87 event horizon telescope results. i. the shadow of the supermassive black hole.  arXiv preprint arXiv:1906.11238, 2019.
  * [9] Ciriaco Goddi et al.  Polarimetric properties of event horizon telescope targets from alma.  The Astrophysical Journal Letters, 910(1):L14, mar 2021.
  * [10] Seiichi Kato, Takahiro Kudoh, and Kazunari Shibata.  Mhd simulation of active galactic nuclei jets.  In Shoken M. Miyama, Kohji Tomisaka, and Tomoyuki Hanawa, editors, Numerical Astrophysics, pages 213–214, Dordrecht, 1999. Springer Netherlands.
  * [11] Yoshiaki Kato, Kazunari Shibata, and Shin Mineshige.  Formation of magnetic-tower jets in magnetohydrodynamic accretion flows.  Journal of the Korean Physical Society, 45(3):1795–1798, 2004.
  * [12] Hongsu Kim and Yanghwan Kim.  Concurrence of the blandford-payne process and the bardeen-petterson effect: Theoretical prediction and its observational evidences, 2015.
  * [13] Shinji Koide, Kazunari Shibata, and Takahiro Kudoh.  General relativistic magnetohydrodynamic simulations of jets from black hole accretion disks: Two-component jets driven by nonsteady accretion of magnetized disks.  The Astrophysical Journal, 495(1):L63–L66, 1998.
  * [14] Josef Lense and Hans Thirring.  Über den einfluss der eigenrotation der zentralkörper auf die bewegung der planeten und monde nach der einsteinschen gravitationstheorie.  Physikalische Zeitschrift, 19:156, 1918.
  * [15] Yosuke Mizuno.  Grmhd simulations and modeling for jet formation and acceleration region in agns.  Universe, 8(2):85, 2022.
  * [16] Lena Murchikova and Gunther Witzel.  Second-scale submillimeter variability of sagittarius a* during flaring activity of 2019: On the origin of bright near-infrared flares.  The Astrophysical Journal Letters, 920(1):L7, 2021.
  * [17] Ramesh Narayan and Jeffrey E McClintock.  Observational evidence for a correlation between jet power and black hole spin.  Monthly Notices of the Royal Astronomical Society: Letters, 419(1):L69–L73, 2012.
  * [18] Jonelle L Walsh, Aaron J Barth, Luis C Ho, and Marc Sarzi.  The m87 black hole mass from gas-dynamical models of space telescope imaging spectrograph observations.  The Astrophysical Journal, 770(2):86, 2013.
  * [19] Daniel C Wilkins.  Bound geodesics in the kerr metric.  Physical Review D, 5(4):814, 1972.

Figure 1: polarized light of the M87 jet observed by ALMA. [9] Figure 2:
Schematic Light Curve of the coupled Blandford Payne process and the Bardeen
Petterson effect.


# Efficient Bayesian reduced rank regression using Langevin Monte Carlo
approach

The Tien Mai(1)

( (1) Oslo Centre for Biostatistics and Epidemiology, Department of
Biostatistics,
University of Oslo, Norway.
Email<EMAIL_ADDRESS>)

###### Abstract

The problem of Bayesian reduced rank regression is considered in this paper.
We propose, for the first time, to use Langevin Monte Carlo method in this
problem. A spectral scaled Student prior distrbution is used to exploit the
underlying low-rank structure of the coefficient matrix. We show that our
algorithms are significantly faster than the Gibbs sampler in high-dimensional
setting. Simulation results show that our proposed algorithms for Bayesian
reduced rank regression are comparable to the state-of-the-art method where
the rank is chosen by cross validation.

##  1 Introduction

Reduced rank regression [Anderson et al., 1951, Izenman, 1975, Velu and
Reinsel, 2013] is a widely used model in linear multivariate regression. In
this model, the low-rank constraint is imposed on the coefficient matrix to
promote estimation and prediction accuracy of multivariate regression. This
low-rank structure is building upon the belief that the response variables are
related to the predictors through only a few latent direction. This assumption
is also useful to extend the model to high dimension settings [Bunea et al.,
2011].

Variouse methods have been conducted for reduced rank regression based on
Bayesian approach. A first study can be found in [Geweke, 1996] in bayesian
econometrics which based on a low-rank matrix factorization approach. Since
then, various works have been studied, for example, [Kleibergen and Paap,
2002, Corander and Villani, 2004, Babacan et al., 2012, Schmidli, 2019,
Alquier, 2013] and developed to low-rank matrix estimation and completion [Lim
and Teh, 2007, Salakhutdinov and Mnih, 2008, Mai and Alquier, 2015]. Moreover,
Bayesian reduced rank regression has been successfully applied in other fields
such as genomics [Marttinen et al., 2014, Zhu et al., 2014]. More recently,
several works have expanded Bayesian methods for incoperating the sparsity
into the reduced rank regression [Goh et al., 2017, Chakraborty et al., 2020,
Yang et al., 2020].

It is noted however that most works in Bayesian reduced rank regression (BRRR)
usually employ conjugate priors as the conditional posterior distributions can
be explitcitely derived that allow to implement Gibbs sampling [Geweke, 1996,
Salakhutdinov and Mnih, 2008]. The details of these priors are reviewed and
discussed in [Alquier, 2013]. Nevertheless, these Gibbs sampling approaches
need to calculate a number of matrix inversions or singular value
decompositions at each iteration and thus it will be costly and can slow down
significantly the algorithm for large data. Different attempts, however based
on optimization, have been made to address this issue including variational
Bayesian methods and maximum a posteriori approach [Lim and Teh, 2007,
Marttinen et al., 2014, Yang et al., 2018].

In this paper, we consider a different way for choosing the prior distribution
on low-rank matrices rather than based on the traditional low-rank matrix
factorization. More specifically, a scaled Student prior is used in our
approach for which the rank of the coefficient matrix does not need to be
prespecified as in the matrix factorization approach. This prior has been
successfully used before in a context related to low-rank matrix completion
[Yang et al., 2018] and image denoising [Dalalyan et al., 2020a].

We develop, for the first time, a Langevin Monte Carlo (MC) approach in
Bayesian reduced rank regression. The Langevin MC method was introduced in
physics based on Langevin diffusions [Ermak, 1975] and it became popular in
statistics and machine learning following the paper [Roberts et al., 1996].
Recent advances in the study of Langevin Monte Carlo make it become more
popular in practice [Dalalyan, 2017, Durmus et al., 2017, Dalalyan et al.,
2020b] and a promissing approach for Bayesian statistics [Durmus et al.,
2019].

More specifically, in this work, we first present a naive (unadjusted)
Langevin MC algorithm and then a Metropolis–Hastings correction for Langevin
MC is proposed. Interestingly, our implementations does not require to perform
matrix inversion nor singular values decomposition and thus our algorithms can
deal with large data set efficiently. Numerical results from these two
algorithm are comparable to the frequentist approach for which the rank is
chosen using cross validation. More particularly, we further show that our
proposed Langevin MC algorithms are significantly faster than Gibbs sampler in
high-dimensional settings, see Section 4.

The paper is structured as follows. In Section 2 we present the reduced rank
regeression model, then the prior distribution is defined together with some
discussion regarding the low-rank factorization prior. In Section 3, the
implementations of the Langevin MC approach are given in details. Numerical
simulations and a real data application are presented in Section 4. Some
discussion and conclusion are given in Section 5 and 6 respectively.

##  2 Bayesian reduced rank regression

###  2.1 Model

We observe two matrices $X$ and $Y$ with

$Y=XB+\mathcal{E}$ (1)

where $Y\in\mathbb{R}^{n\times m}$ is the response matrix,
$X\in\mathbb{R}^{n\times p}$ is the predictor matrix, $B\in\mathbb{R}^{p\times
m}$ is the unknown regression coefficient matrix and $\mathcal{E}$ is an
$n\times m$ random noise matrix with $\mathbb{E}(\mathcal{E})=0$.

We assume that the entries $\mathcal{E}_{i,j}$ of $\mathcal{E}$ are i.i.d.
Gaussian $\mathcal{N}(0,\sigma^{2})$ and $\sigma^{2}$ is given. In this case,
note that the likelihood distribution is given by

$L(Y)\propto\exp\left\\{-\frac{1}{2\sigma^{2}}\|Y-XB\|^{2}_{F}\right\\}$

where $\|\cdot\|_{F}$ denotes the Frobenius norm, $\|M\|_{F}^{2}={\rm
Tr}(M^{T}M)$.

The objective is to estimate the parameter matrix $B$. In many applications,
it makes sense to assume that the matrix $B$ has low rank, i.e. ${\rm
rank}(B)\ll\min(p,m)$.

With a prior distribution $\pi(B)$, the posterior for the Bayesian reduced
rank regression is

$\displaystyle\mathcal{L}_{n}(B)\propto L(Y)\pi(B).$

Then, the Bayesian estimator is defined as

$\displaystyle\widehat{B}=\int B\mathcal{L}_{n}(B).$ (2)

Note that here we focus on estimating $B$ and thus $\sigma^{2}$ is fixed. On
can relax this assumption and place a prior distribution on $\sigma$.

###  2.2 A low-rank promoting prior

Borrowing motivation from a low-rank prior explored in recent works [Yang et
al., 2018, Dalalyan et al., 2020a], we consider the following prior,

$\displaystyle\pi(B)\propto\det(\lambda^{2}\mathbf{I}_{p}+BB^{\intercal})^{-(p+m+2)/2}$
(3)

where $\lambda>0$ is a tuning parameter and $\mathbf{I}_{p}$ is the $p\times
p$ identity matrix .

To illustrate that this prior has the potential to encourage the low-rankness
of $B$, one can check that

$\displaystyle\pi(B)\propto\prod_{j=1}^{m}(\lambda^{2}+s_{j}(B)^{2})^{-(p+m+2)/2},$

where $s_{j}(B)$ denotes the $j$the largest singular value of $B$. It is well
known that the log-sum function $\sum_{j=1}^{m}\log(\lambda^{2}+s_{j}(B)^{2})$
encourages a sparsity on $\\{s_{j}(B)\\}$, see [Candes et al., 2008, Yang et
al., 2018]. Thus the resulting matrix $B$ has a low-rank structrure,
approximately.

The following Lemma explains the reason why this prior is a spectral scaled
Student prior distribution.

######  Lemma 1.

_[Dalalyan et al., 2020a] _ If a random matrix $B$ has the density
distribution as in (3), then the random vectors $B_{i}$ are all drawn from the
$p$-variate scaled Student distribtuion $(\lambda/\sqrt{3})t_{3,p}$.

#### On low-rank factorization priors

The first idea about a low-rank prior was carried out in [Geweke, 1996]. That
is to express the matrix parameter $B$ as $B_{p\times m}=M_{p\times
k}N_{m\times k}^{\top}$ with $k\leq\min(p,m)$. The prior is difined on $M$ and
$N$ rather than on $B$ as

$\pi(M,N)\propto\exp\left\\{-\frac{\tau^{2}}{2}(\|M\|^{2}_{F}+\|N\|^{2}_{F})\right\\}$

for some $\tau>0$. This prior allows to obtain explicit forms for the marginal
posteriors that allows an implementation of the Gibbs algorithm to sample from
the posterior, see [Geweke, 1996]. However, the downside of this approach is
the problem of choosing $k$, the reduced rank, is not directly addressed.
Thus, one has to perform model selection for any possible $k$, as done in
[Kleibergen and Paap, 2002, Corander and Villani, 2004].

Recent approaches focus on fixing a large $k$, e.g. $k=\min(p,m)$, then
sparsity-promoting priors are placed on the columns of $M$ and $N$ such that
most columns are almost null. So that the resulting matrix $B=MN^{\top}$ is
approximately low-rank. This direction was first proposed in [Babacan et al.,
2012] in the context of matrix completion, but it still can be used in reduced
rank regression. See [Alquier, 2013] for the details and dicussions on low-
rank factorization priors.

With low-rank factorization priors, most authors simulate from the posterior
by using the Gibbs sampler as the conditional posterior distributions can be
explitcitely derived, e.g. in [Geweke, 1996, Salakhutdinov and Mnih, 2008].
However, these Gibbs sampling algorithms update the factor matrices in a row-
by-row fashion and invole a number of matrix inverse operations at each
iteration. This is expensive and slow down the algorithm for large data set.

##  3 Langevin Monte Carlo implementations

In this section, we propose to compute an approximation of the posterior with
the scaled multivariate Student prior by a suitable version of the Langevin
Monte Carlo algorithm.

###  3.1 Unadjusted Langevin Monte Carlo algorithm

Let us recall that the log posterior is of the following form

$\log\mathcal{L}_{n}(B)=-\frac{1}{2\sigma^{2}}\|Y-XB\|_{F}^{2}-\frac{p+m+2}{2}\log\det(\lambda^{2}\mathbf{I}_{p}+BB^{\intercal}),$

and consequently,

$\nabla\log\mathcal{L}_{n}(B)=-\frac{1}{\sigma^{2}}X^{\intercal}(Y-XB)-(p+m+2)(\lambda^{2}\mathbf{I}_{p}+BB^{\intercal})^{-1}B,$

We use the constant step-size unadjusted Langevin MC (denoted by LMC) [Durmus
et al., 2019]. It is defined by choosing an initial matrix $B_{0}$ and then by
using the recursion

$\displaystyle
B_{k+1}=B_{k}-h\nabla\log\mathcal{L}_{n}(B_{k})+\sqrt{2h}\,W_{k},\qquad
k=0,1,\ldots,$ (4)

where $h>0$ is the step-size and $W_{0},W_{1},\ldots$ are independent random
matrices with i.i.d. standard Gaussian entries. The detail of the algorithm is
given in the Algorithm 1.

Note that a direct application of the Langevin MC algorithm (4) needs to
calculate a $p\times p$ matrix inversion at each iteration. This can slow down
significantly the algorithm and might be expensive. However, one can easily
verify that the matrix
$\mathbf{M}=(\lambda^{2}\mathbf{I}_{p}+BB^{\intercal})^{-1}B$ is the solution
to the following convex optimization problem

$\min_{\mathbf{M}}\big{\\{}\|\mathbf{I}_{m}-B^{\top}\mathbf{M}\|_{F}^{2}+\lambda^{2}\|\mathbf{M}\|_{F}^{2}\big{\\}}.$

The solution of this optimization problem can be obtained by using the package
’glmnet’ [Friedman et al., 2010] (with family option ’mgaussian’). This does
not require matrix inversion nor other costly operation. However, it is noted
that in this case we are using the Langevin MC with approximate gradient
evaluation, theoretical assessment of this method can be found in [Dalalyan
and Karagulyan, 2019].

Algorithm 1 LMC for BRRR

1:Input: matrices $Y\in\mathbb{R}^{n\times m},X\in\mathbb{R}^{n\times p}$

2:Parameters: Positive real numbers $\lambda,h,T$

3:Onput: The matrix $\widehat{B}$

4:Initialize:
$B_{0}\leftarrow(X^{\top}X+0.1\mathbf{I}_{p})^{-1}X^{\top}Y;\widehat{B}=\mathbf{0}_{p\times
m}$

5:for $k\leftarrow 1$ to $T$ do

6: Simulate $B_{k}$ from (4);

7: $\widehat{B}\leftarrow\widehat{B}+B_{k}/T$

8:end for

######  Remark 1.

It seems that the Algorithm 1 looks like an iterative gradient descent for
minimizing the penalized Gaussian log-likelihood with a penalty. However,
Algorithm 1 computes the posterior mean and not a maximum a posteriori
estimator. More precisely, our algorithm includes a final step of averaging.

######  Remark 2.

For small values of $h$, the ouput $\widehat{B}$ is very close to the integral
(2) of interest. However, for some $h$ that may not small enough, the Markov
process is transient and thus the sum explodes [Roberts and Stramer, 2002]. To
address this problem, one have to take a smaller h and restart the algorithm
or a Metropolis–Hastings correction can be included in the algortihm. The
Metropolis–Hastings approach ensures the convergence to the desired
distribution, however, it greatly slows down the algorithm because of an
additional acception/rejection step at each iteration. The approach by taking
a smaller $h$ also slows down the algorithm but we keep some control on its
time of execution.

######  Remark 3.

Based our observations from numerical studies in Section 4, the initial matrix
$B_{0}$ can also effect the convergence of the algorithm. We suggest using
$B_{0}=(X^{\top}X+0.1\mathbf{I}_{p})^{-1}X^{\top}Y$ as a default alternative.

###  3.2 A Metropolis-adjusted Langevin algorithm

Here, we consider a Metropolis-Hasting correction to the Algorithm 1. This
approach guarantees the convergence to the posterior. More precisely, we
consider the update rule in (4) as a proposal for a new state,

$\displaystyle\widetilde{B}_{k+1}=B_{k}-h\nabla\log\mathcal{L}_{n}(B_{k})+\sqrt{2h}\,W_{k},\qquad
k=0,1,\ldots.$ (5)

Note that $\widetilde{B}_{k+1}$ is normally distributed with mean
$B_{k}-h\nabla\log\mathcal{L}_{n}(B_{k})$ and the covariance matrix equals to
$2h\mathbf{I}_{p}$. This proposal is accepted or rejected according to the
Metropolis-Hastings algorithm. That is the proposal is accepted with
probabiliy:

$\min\left\\{1,\frac{\mathcal{L}_{n}(\widetilde{B}_{k+1})q(B_{k}|\widetilde{B}_{k+1})}{\mathcal{L}_{n}(B_{k})q(\widetilde{B}_{k+1}|B_{k})}\right\\},$

where

$q(x^{\prime}|x)\propto\exp\left(-\frac{1}{4h}\|x^{\prime}-x+h\nabla\log\mathcal{L}_{n}(x)\|^{2}_{F}\right)$

is the transition probability density from $x$ to $x^{\prime}$. The detail of
the Metropolis-adjusted Langevin algorithm (denoted by MALA) for BRRR is given
Algorithm 2. Compared to random-walk Metropolis–Hastings, MALA has the
advantage that it usually proposes moves into regions of higher probability,
which are then more likely to be accepted.

######  Remark 4.

Following [Roberts and Rosenthal, 1998], the choice of the step-size $h$ is
tuned such that the acceptance rate is approximate $0.5$. See Section 4 for
some choices in special cases in our simulations.

Algorithm 2 MALA for BRRR

1:Input: matrices $Y\in\mathbb{R}^{n\times m},X\in\mathbb{R}^{n\times p}$

2:Parameters: Positive real numbers $\lambda,h,T$

3:Onput: The matrix $\widehat{B}$

4:Initialize:
$B_{0}\leftarrow(X^{\top}X+0.1\mathbf{I}_{p})^{-1}X^{\top}Y;\widehat{B}=\mathbf{0}_{p\times
m}$

5:for $k=1$ to $T$ do

6: Simulate $\widetilde{B}_{k}$ from (5)

7: Calculate
$\alpha=\min\left\\{1,\frac{\mathcal{L}_{n}(\widetilde{B}_{k})q(B_{k-1}|\widetilde{B}_{k})}{\mathcal{L}_{n}(B_{k})q(\widetilde{B}_{k}|B_{k-1})}\right\\}$

8: Sample $u\sim U[0,1]$

9: if $u\leq\alpha$ then

10: $B_{k}=\widetilde{B}_{k}$

11: else

12: $B_{k}=B_{k-1}$

13: end if

14: $\widehat{B}\leftarrow\widehat{B}+B_{k}/T$

15:end for

##  4 Numerical studies

###  4.1 Simulations setups and details

First, we perform some numerical studies on simulated data to access the
performance of our proposed algorithms. We consider the following model
setups:

  * •

Model I: A low-dimensional set up is studied with $n=100,p=12,m=8$ and the
true rank $r={\rm rank}(B)=3$. The design matrix $X$ is generated from
$\mathcal{N}(0,\Sigma)$ where the covariance matrix $\Sigma$ is with diagonal
entries 1 and off-diagonal entries $\rho_{X}\geq 0$. We consider $\rho_{X}=0$
and $\rho_{X}=0.5$, this creates a wide-range correlation in the predictors.
The true coefficient matrix is generated as $B=B_{1}B_{2}^{\top}$ where
$B_{1}\in\mathbb{R}^{p\times r},B_{2}\in\mathbb{R}^{m\times r}$ and all
entries in $B_{1}$ and $B_{2}$ are randomly sampled from $\mathcal{N}(0,1)$.

  * •

Model II: This model is similar to Model I, however, a high-dimensional set up
is considered with $n=100,p=150,m=90$.

  * •

Model III: An approximate low-rank set up is studied. This series of
simulation is similar to the Model II, except that the true coefficient is no
longer rank 3, but it can be well approximated by a rank 3 matrix:

$B=2\cdot B_{1}B_{2}^{\top}+E,$

where $E$ is matrix with entries sampled from $\mathcal{N}(0,1)$.

Under each setting, the entire data generation process described above is
replicated 100 times.

We compare our algorithms LMC, MALA to the naive reduced rank method (denoted
RRR, see [Velu and Reinsel, 2013]) where the rank is selected by 10-fold cross
validation. The RRR method is available from the R package
’rrpack’111https://cran.r-project.org/package=rrpack. We also compare LMC and
MALA to the Gibbs sampler from [Alquier, 2013], however, we are just able to
perform these comparisions in Model I as the Gibbs sampler is too slow for
large dimensions, see Figure 1. The R codes for the Gibbs sampler are kindly
provided by the author of the paper [Alquier, 2013].

The evaluations are done by using the estimation error (Est) and the
normalized mean square error (Nmse)

${\rm Est}:=\|B-\widehat{B}\|^{2}_{F}/(pm),\quad{\rm
Nmse}:=\|B-\widehat{B}\|^{2}_{F}/\|B\|^{2}_{F},$

that are calculated as the average of the mean squared errors from all 100
runs. We also evaluate the average over 100 runs of the prediction error
(Pred) as

${\rm Pred}:=\|Y_{test}-X_{test}\widehat{B}\|^{2}_{F}/(nm)$

where $X_{test}$ is a newly generated $n\times p$ test-sample matrix of
predictors and $Y_{test}$ is a newly generated $n\times m$ test-sample matrix
of responses. We also report the average of estimated rank, denote Rank, for
different methods over all the runs.

The choice of the step-size parameters is set as: for Model III, we take
$h=3/(\sqrt{m}np)$; with Model II $h=5/(mnp)$ and with Model I,
$h=2/(pm\sqrt{n})$. This choice is selected such that the acceptance rate of
MALA is approximate 0.5. We fixed $\lambda=3$ in all models. The LMC, MALA and
Gibbs sampler are run with $T=200$ iterations and we take the first $100$
steps as burn-in.

###  4.2 Simulation results

In low-dimensional setting where $pm<n$ as in Model I, Langevin MC algorithms
(LMC, MALA) are able to recover the true rank of the model, see Table 1. The
results of MALA are slightly better than LMC. The prediction errors of LMC and
MALA are comparable to RRR and Gibbs sampler. In terms of other errors (Est
and Nmse), LMC and MALA are twice worse than RRR and Gibbs sampler.

However, it is worth noting that the running time of our algorithms is
linearly with $p$ where $n$ and $m$ are fixed, while the Gibbs sampler is not.
More specifically, we conducted a comparision on the running time for these
four algortihms where the dimension $p$ is varied by $10,50,100,150$ with
fixed $n=100,m=90$. The results is given in Figure 1. It is clear that the
Gibbs sampler is several magnitude slower than our algorithms.

Results from high dimensional settings as in Model II and II reveal that our
algorithms perform quite similar the RRR method, see Table 2, in term of all
considered errors. Moreover, it is interesting that LMC and MALA are slightly
better than RRR method in Model III where coefficient matrix is approximately
low-rank. More specifically, MALA produces promissing results that are lightly
better than LMC as well as RRR.

Table 1: Simulation results on simulated data in Model I for different methods, with their standard error in parentheses. (Est: average of estimation error; Pred: average of prediction error; Rank: average of estimated rank). | $\rho_{X}=0.0$
---|---
Errors | LMC | MALA | RRR | Gibbs
$10^{2}\times$Est | 1.25 (0.21) | 1.26 (0.21) | 0.58 (0.14) | 0.59 (0.13)
Pred | 1.15 (0.06) | 1.15 (0.06) | 1.07 (0.05) | 1.07 (0.05)
$10^{3}\times$Nmse | 4.66 (2.21) | 4.80 (2.46) | 2.16 (1.11) | 2.18 (1.06)
Rank | 3 (0.0) | 3 (0.0) | 3 (0.0) | 3 (0.0)
| $\rho_{X}=0.5$
| LMC | MALA | RRR | Gibbs
$10^{2}\times$Est | 2.52 (0.41) | 2.39 (0.56) | 1.03 (0.24) | 1.12 (0.02)
Pred | 1.17 (0.07) | 1.16 (0.07) | 1.08 (0.06) | 1.08 (0.07)
$10^{2}\times$Nmse | 0.94 (0.41) | 0.96 (0.55) | 0.43 (0.31) | 0.46 (0.29)
Rank | 3 (0.0) | 3 (0.0) | 3 (0.0) | 3 (0.0)
Table 2: Simulation results on simulated data in Model II & III for different methods, with their standard error in parentheses. (Est: average of estimation error; Pred: average of prediction error; Rank: average of estimated rank). Models |  Errors | LMC | MALA | RRR
---|---|---|---|---
|  Est | 1.00 (0.13) | 1.00 (0.13) | 0.99 (0.13)
II, $\rho_{X}=0.0$ |  $10^{-2}\times$Pred | 1.51 (0.23) | 1.51 (0.23) | 1.49 (0.23)
|  Nmse | 0.34 (0.03) | 0.34 (0.03) | 0.33 (0.03)
|  Rank | 3 (0.0) | 3 (0.0) | 3 (0.0)
|  Est | 1.01 (0.14) | 1.01 (0.14) | 0.99 (0.14)
II, $\rho_{X}=0.5$ |  $10^{-1}\times$Pred | 7.65 (1.20) | 7.65 (1.20) | 7.45 (1.20)
|  Nmse | 0.34 (0.03) | 0.34 (0.03) | 0.33 (0.03)
|  Rank | 3 (0.0) | 3 (0.0) | 3 (0.0)
|  Est | 4.32 (0.64) | 4.32 (0.64) | 4.34 (0.64)
III, $\rho_{X}=0.0$ |  $10^{-2}\times$ Pred | 6.53 (1.05) | 6.53 (1.05) | 6.56 (1.05)
|  Nmse | 0.33 (0.03) | 0.33 (0.03) | 0.33 (0.03)
|  Rank | 3.04 (0.20) | 3.04 (0.20) | 3.04 (0.20)
|  Est | 4.41 (0.58) | 4.40 (0.57) | 4.43 (0.57)
III, $\rho_{X}=0.5$ |  $10^{-2}\times$ Pred | 3.31 (0.47) | 3.30 (0.47) | 3.32 (0.47)
|  Nmse | 0.34 (0.03) | 0.34 (0.03) | 0.34 (0.03)
|  Rank | 3.09 (0.35) | 3.08 (0.34) | 3.08 (0.34)
Figure 1: Plot to compare the running times for 10 iterations of LMC, MALA
Gibbs sampler and 10-fold cross validation RRR with fixed $n=100,m=90,r=2$ and
the dimension $p$ is varied.

###  4.3 Real data application

We apply our algorithms to a breast cancer dataset [Witten et al., 2009] to
acess its performance on real data set. This data consisting of gene
expression measurements and comparative genomic hybridization measurements for
$n=89$ samples. The dataset is available from the R packge ’PMA’ [Witten et
al., 2009]. This data were used before in [Chen et al., 2013] in the context
of reduced rank regression.

Following [Chen et al., 2013], we consider the gene expression profiles of a
chromosome as predictors and the copy-number variations of the same chromosome
as repsonse. The analysis is focused on chromosome $21$, for which $m=44$ and
$p=227$. The data are randomly divided into a training set of size
$n_{train}=79$ and a test set of size $n_{test}=10$. Model estimation is done
by using the training data. Then the predictive performance is calculated on
the test data by its mean squared prediction error
$\|Y_{test}-X_{test}\widehat{B}\|_{F}^{2}/(mn_{test})$, where
$(Y_{test},X_{test})$ denotes the test set. We repeat the random training/test
spliting process $100$ times and report the average mean squared prediction
error and the average rank estimate for each method. The results are given in
Table 3, we can see that MALA is better than LMC.

Table 3: Comparision of the model fits to the real data. The mean suqared prediction errors (MSPE) and the estimated ranks are reported, with their standard error in parentheses. | LMC | MALA | RRR
---|---|---|---
MSPE | 0.052 (.009) | 0.049 (.008) | 0.030 (.008)
Rank | 1.03 (.17) | 1.03 (.17) | 0.74 (.44)

##  5 Discussion

It is noted that our Langevin MC approaches for BRRR are using a different
prior on (approximate) low-rank matrix comparing with the matrix factorization
as in Gibbs sampler. There are several other ways to define such priors on a
whole matrix, see [Sundin, 2016]. For example, one could consider, with
$\lambda>0$,

$\pi(B)\propto\exp(\lambda{\rm Trace}(BB^{\top})^{1/2})$

where its log-prior is the nuclear norm that is also promoting the low-rank
structure on $B$. The application of Langevin MC method for such priors would
be interesting research directions in the future.

A vital part in the Langevin MC approach is choosing the step-size $h$. Here,
in this work, the choice of $h$ is picked such that the acceptance rate in
MALA is around $0.5$, motivating from [Roberts and Rosenthal, 1998]. We have
tried with a decreasing step-size as $h_{t}=h_{t-1}/t$, however this choice
does not improve the results at all compare to the choise defining through the
acceptance rate. It is noted that there are several other way for choosing
$h$, for example, $h$ is adaptively changed in each iteration as in [Marshall
and Roberts, 2012]. The study of such approach to BRRR is left for future
research.

Bayesian studies that incorporating sparsity into RRR model to account for
both rank selection and variable selection have been recently carried out in
[Goh et al., 2017, Yang et al., 2020]. However, these works are still based on
low-rank factorization priors and the implementation of the Gibbs sampler.
Thus, the application of Langevin MC to this problem would be another
interesting future work.

##  6 Conclusion

In the paper, we have proposed efficient Langevin MC approach for BRRR. The
performances of our algorithms are similar to the state-of-the-art method in
simulations. More importantly, we showed that the proposed algorithms are
significant faster than the Gibbs sampler. This is an interesting way that
makes BRRR become more applicable in large data set.

## Availability of data and materials

The R codes and data used in the numerical experiments are available at:
https://github.com/tienmt/BRRR .

## Acknowledgments

This research of T.T.M was supported by the European Research Council grant
no. 742158.

## References

  * Alquier, 2013 Alquier, P. (2013).  Bayesian methods for low-rank matrix estimation: short survey and theoretical study.  In International Conference on Algorithmic Learning Theory, pages 309–323. Springer.
  * Anderson et al., 1951 Anderson, T. W. et al. (1951).  Estimating linear restrictions on regression coefficients for multivariate normal distributions.  Annals of mathematical statistics, 22(3):327–351.
  * Babacan et al., 2012 Babacan, S. D., Luessi, M., Molina, R., and Katsaggelos, A. K. (2012).  Sparse bayesian methods for low-rank matrix estimation.  IEEE Transactions on Signal Processing, 60(8):3964–3977.
  * Bunea et al., 2011 Bunea, F., She, Y., Wegkamp, M. H., et al. (2011).  Optimal selection of reduced rank estimators of high-dimensional matrices.  The Annals of Statistics, 39(2):1282–1309.
  * Candes et al., 2008 Candes, E. J., Wakin, M. B., and Boyd, S. P. (2008).  Enhancing sparsity by reweighted $\ell_{1}$ minimization.  Journal of Fourier analysis and applications, 14(5-6):877–905.
  * Chakraborty et al., 2020 Chakraborty, A., Bhattacharya, A., and Mallick, B. K. (2020).  Bayesian sparse multiple regression for simultaneous rank reduction and variable selection.  Biometrika, 107(1):205–221.
  * Chen et al., 2013 Chen, K., Dong, H., and Chan, K.-S. (2013).  Reduced rank regression via adaptive nuclear norm penalization.  Biometrika, 100(4):901–920.
  * Corander and Villani, 2004 Corander, J. and Villani, M. (2004).  Bayesian assessment of dimensionality in reduced rank regression.  Statistica Neerlandica, 58(3):255–270.
  * Dalalyan, 2017 Dalalyan, A. S. (2017).  Theoretical guarantees for approximate sampling from smooth and log-concave densities.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 3(79):651–676.
  * Dalalyan et al., 2020a Dalalyan, A. S. et al. (2020a).  Exponential weights in multivariate regression and a low-rankness favoring prior.  In Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, volume 56, pages 1465–1483. Institut Henri Poincaré.
  * Dalalyan and Karagulyan, 2019 Dalalyan, A. S. and Karagulyan, A. (2019).  User-friendly guarantees for the langevin monte carlo with inaccurate gradient.  Stochastic Processes and their Applications, 129(12):5278–5311.
  * Dalalyan et al., 2020b Dalalyan, A. S., Riou-Durand, L., et al. (2020b).  On sampling from a log-concave density using kinetic langevin diffusions.  Bernoulli, 26(3):1956–1988.
  * Durmus et al., 2017 Durmus, A., Moulines, E., et al. (2017).  Nonasymptotic convergence analysis for the unadjusted langevin algorithm.  The Annals of Applied Probability, 27(3):1551–1587.
  * Durmus et al., 2019 Durmus, A., Moulines, E., et al. (2019).  High-dimensional bayesian inference via the unadjusted langevin algorithm.  Bernoulli, 25(4A):2854–2882.
  * Ermak, 1975 Ermak, D. L. (1975).  A computer simulation of charged particles in solution. i. technique and equilibrium properties.  The Journal of Chemical Physics, 62(10):4189–4196.
  * Friedman et al., 2010 Friedman, J., Hastie, T., and Tibshirani, R. (2010).  Regularization paths for generalized linear models via coordinate descent.  Journal of Statistical Software, 33(1):1–22.
  * Geweke, 1996 Geweke, J. (1996).  Bayesian reduced rank regression in econometrics.  Journal of econometrics, 75(1):121–146.
  * Goh et al., 2017 Goh, G., Dey, D. K., and Chen, K. (2017).  Bayesian sparse reduced rank multivariate regression.  Journal of multivariate analysis, 157:14–28.
  * Izenman, 1975 Izenman, A. J. (1975).  Reduced-rank regression for the multivariate linear model.  Journal of multivariate analysis, 5(2):248–264.
  * Kleibergen and Paap, 2002 Kleibergen, F. and Paap, R. (2002).  Priors, posteriors and bayes factors for a bayesian analysis of cointegration.  Journal of Econometrics, 111(2):223–249.
  * Lim and Teh, 2007 Lim, Y. J. and Teh, Y. W. (2007).  Variational bayesian approach to movie rating prediction.  In Proceedings of KDD cup and workshop, volume 7, pages 15–21. Citeseer.
  * Mai and Alquier, 2015 Mai, T. T. and Alquier, P. (2015).  A bayesian approach for noisy matrix completion: Optimal rate under general sampling distribution.  Electron. J. Statist., 9(1):823–841.
  * Marshall and Roberts, 2012 Marshall, T. and Roberts, G. (2012).  An adaptive approach to langevin mcmc.  Statistics and Computing, 22(5):1041–1057.
  * Marttinen et al., 2014 Marttinen, P., Pirinen, M., Sarin, A.-P., Gillberg, J., Kettunen, J., Surakka, I., Kangas, A. J., Soininen, P., O’Reilly, P., Kaakinen, M., et al. (2014).  Assessing multivariate gene-metabolome associations with rare variants using bayesian reduced rank regression.  Bioinformatics, 30(14):2026–2034.
  * Roberts and Rosenthal, 1998 Roberts, G. O. and Rosenthal, J. S. (1998).  Optimal scaling of discrete approximations to langevin diffusions.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(1):255–268.
  * Roberts and Stramer, 2002 Roberts, G. O. and Stramer, O. (2002).  Langevin diffusions and metropolis-hastings algorithms.  Methodology and computing in applied probability, 4(4):337–357.
  * Roberts et al., 1996 Roberts, G. O., Tweedie, R. L., et al. (1996).  Exponential convergence of langevin distributions and their discrete approximations.  Bernoulli, 2(4):341–363.
  * Salakhutdinov and Mnih, 2008 Salakhutdinov, R. and Mnih, A. (2008).  Bayesian probabilistic matrix factorization using markov chain monte carlo.  In Proceedings of the 25th international conference on Machine learning, pages 880–887.
  * Schmidli, 2019 Schmidli, H. (2019).  Bayesian reduced rank regression for classification.  In Applications in Statistical Computing, pages 19–30. Springer.
  * Sundin, 2016 Sundin, M. (2016).  Bayesian methods for sparse and low-rank matrix problems.  PhD thesis, KTH Royal Institute of Technology.
  * Velu and Reinsel, 2013 Velu, R. and Reinsel, G. C. (2013).  Multivariate reduced-rank regression: theory and applications, volume 136.  Springer Science & Business Media.
  * Witten et al., 2009 Witten, D. M., Tibshirani, R., and Hastie, T. (2009).  A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis.  Biostatistics, 10(3):515–534.
  * Yang et al., 2020 Yang, D., Goh, G., and Wang, H. (2020).  A fully bayesian approach to sparse reduced-rank multivariate regression.  Statistical Modelling, page 1471082X20948697.
  * Yang et al., 2018 Yang, L., Fang, J., Duan, H., Li, H., and Zeng, B. (2018).  Fast low-rank bayesian matrix completion with hierarchical gaussian prior models.  IEEE Transactions on Signal Processing, 66(11):2804–2817.
  * Zhu et al., 2014 Zhu, H., Khondker, Z., Lu, Z., and Ibrahim, J. G. (2014).  Bayesian generalized low rank regression models for neuroimaging phenotypes and genetic markers.  Journal of the American Statistical Association, 109(507):977–990.


# Towards a privacy-preserving distributed cloud service for preprocessing
very large medical images

Yuandou Wang1, Neel Kanwal2, Kjersti Engan2, Chunming Rong2, Zhiming Zhao1,3
1Multiscale Networked Systems, University of Amsterdam, the Netherlands
2Department of Electrical Engineering and Computer Science, University of
Stavanger, Norway
3LifeWatch ERIC Virtual Lab and Innovation Center (VLIC), Amsterdam, the
Netherlands
Email: {y.wang8<EMAIL_ADDRESS>{neel.kanwal, kjersti.engan,
<EMAIL_ADDRESS>

###### Abstract

Digitized histopathology glass slides, known as Whole Slide Images (WSIs), are
often several gigapixels large and contain sensitive metadata information,
which makes distributed processing unfeasible. Moreover, artifacts in WSIs may
result in unreliable predictions when directly applied by Deep Learning (DL)
algorithms. Therefore, preprocessing WSIs is beneficial, e.g., eliminating
privacy-sensitive information, splitting a gigapixel medical image into tiles,
and removing the diagnostically irrelevant areas. This work proposes a cloud
service to parallelize the preprocessing pipeline for large medical images.
The data and model parallelization will not only boost the end-to-end
processing efficiency for histological tasks but also secure the
reconstruction of WSI by randomly distributing tiles across processing nodes.
Furthermore, the initial steps of the pipeline will be integrated into the
Jupyter-based Virtual Research Environment (VRE) to enable image owners to
configure and automate the execution process based on resource allocation.

###### Index Terms:

Computational Pathology, Cloud Computing, Privacy-preserving, Image
Preprocessing, Virtual Research Environment, Infrastructure Planning

##  I Introduction

Deep Learning (DL) approaches have advanced and innovated automatic
diagnostics, such as quantifying the presence of cancerous cells in digitized
histopathology glass slides, called Whole Slide Images (WSIs). However,
running these diagnostic services over a large scale requires a significant
infrastructure capacity for storing and processing images using complex DL
models, e.g., cloud computing and High-Performance computing (HPC), are often
needed. The owners of the medical images, e.g., hospitals, often do not have
such an infrastructure and have to rely on collaborators with remote
infrastructure resources.

Establishing a DL-based pipeline for medical images on a remote infrastructure
is challenging; for instance, 1) WSIs often contain privacy-sensitive
information in their metadata and cannot be directly sent out to the public
cloud from the hospitals; 2) WSIs are often very large and require high
network bandwidth to upload; and 3) WSIs are split into tiles to process [1,
2] and require specialized hardware, e.g., GPUs, to run complex DL models.
Furthermore, it is often complicated to deploy an end-to-end pipeline and
create an efficient re-configurable workflow[3] on remote infrastructure.

Figure 1: An overview of the preprocessing pipeline for detecting artifacts in
whole slide images.

In this paper, we will tackle these challenges by proposing a cloud-based
service that will be integrated into a collaborative virtual research
environment based on the works [4, 5], and we will present a use case of a
medical image preprocessing application from digital pathology domain to
testify our methodology.

Figure 2: The methodology proposed in this work. Step ❶ - preparation of data
and workflow; step ❷ - resource allocation; step ❸ - deployment and execution;
step ❹ - load balance; and step ❺ - aggregating distributed data in a summary
view to visualize the final results.

##  II Case Study

Digital pathology overcomes the hurdles of traditional histopathology by
facilitating the diagnostic process using a WSI [1]. The preparation of
histological glass slides may result in the appearance of artifacts on the
obtained WSI due to improper handling of the tissue specimen during the tissue
processing stages. These histological artifacts are diagnostically irrelevant
and are usually ignored by pathologists in the diagnosis process [6, 2].
Therefore, it is vital to detect and remove them before applying diagnostic or
prognostic algorithms. Some frequently appearing artifacts are damaged tissue,
folded tissue, blur, air bubbles, and diagnostically irrelevant blood [6, 2].
Computational pathology (CPATH) researchers may run DL-based artifact
preprocessing algorithms over thousands of WSIs before applying diagnostic
algorithms, requiring powerful computational resources to process WSIs
efficiently. Fig. 1 presents an overview of such artifact preprocessing
pipeline, which is an ensemble of five DL models for blood, blur, damaged
tissue, folded tissue, and air bubbles detection tasks from WSI in a binary
fashion.

Traditionally, the artifact preprocessing pipeline runs over a single machine,
bringing the disadvantages of a single security breach or system failure.
Besides, handling gigapixel WSIs is time-consuming and resource-intensive,
which raises the demand for parallel distributed computing. Nevertheless, WSIs
processed on private clouds in research environments are de-identified or
pseudonymized under various regulations. This raises concerns about the
embedded privacy-sensitive metadata while making distributed processing over
public clouds.

##  III Methodology

We introduce a methodology to cope with the highlighted challenges, as shown
in Fig. 2. It consists of five main steps: preparation of data and workflow,
resource allocation, deployment and execution, load balance, and data
aggregation.

###  III-A Data and workflow preparation

Step ❶ aims to introduce parallelism and encryption available for the next
steps. To guarantee privacy-preserving requirements, we remove the _metadata_
from WSI before splitting the gigapixel image into many image tiles to
introduce data parallelism. Meanwhile, containerizing the computational tasks
as several reusable fine-grained services can improve scalability and security
since they are isolated from each other and from the host system. Besides, we
apply a matrix $A_{x}$ to record the distribution of the tiles over a grid
which can be encoded to hide its coordinates and divided into sub-matrices
$A_{e,1},A_{e,2},...,A_{e,K}$. Each sub-matrix is considered as an index of
distributed dataset $D_{e,k}$ for each service-based task.

###  III-B Resource allocation

Based on the prepared data and service-based tasks, step ❷ is to map available
resources (e.g., clusters at universities and commercial clouds) to various
tasks in a manner that optimizes their utilization and satisfies requirements
[7]. Related methods, such as IC-PCP [8] and machine learning-based approach
[9], can be improved for workflow scheduling. For bi-objective optimization,
such as reducing execution time and cost, there are trade-offs with time
performance and monetary cost over the cloud. On this basis, this paper looks
into workflow scheduling problems under the influence of privacy requirements
and the split data sets, so the research problem is more challenging.

###  III-C Deployment and execution

After a deployment plan is created at step ❷, the datasets and service-based
tasks will be assigned to planned infrastructures equipped with computation,
communication, and storage resources. The system should ensure that data
storage and task execution remain in place and continue to be effective even
among changes (such as downtime, errors, or attackers) to the system or
emerging threats, according to step ❸. Due to distributed processing, it
reduces the burden of a single machine and avoids a single security breach.

###  III-D Load Balance

Considering that computing nodes may unpredictably slow down or fail during
their execution, step ❹ aims to improve the performance, reliability, and load
balance of task-based applications [10]. This approach asymptotically achieves
near-ideal load balancing and computation cost in the presence of slow nodes
(stragglers), which could also be complementary to workflow scheduling.

###  III-E Data aggregation

Step ❺ takes the predicted distributed output, $b_{o,k}$ for each encoded
tile, from step ❸ for every service task, to reconstruct the encoded
distributions as $A_{d,k}$. Privacy preservation can be guaranteed by the
random-value perturbation technique [11]. This approach has solid theoretical
foundations and is easier to apply for the reconstruction of the encrypted
data than others (e.g., differential privacy and secure multiparty computation
[12]) especially considering a data matrix manipulation. Then by tracing back
coordinates, we can create a segmentation mask for detected artifacts. It
presents the results of the DL-based artifacts detection in a summary view,
incl. visualization, evaluation, and metrics.

##  IV System Model

We sketch out a privacy-preserving distributed processing pipeline for the
medical application, shown in Fig. 3. It is composed of three main steps - viz
_Splitting_ , _Computation_ , and _Aggregation_. Both Splitting (See in step
❶) and Aggregation (See in step ❺) will be executed on a _Trusted Server_.
Such distributed data processing application can be defined as a tuple:

$\mathcal{A}=(\mathcal{M},\varepsilon,\mathcal{D},\mathcal{R},\mathcal{I},req)$
(1)

where $\mathcal{M}$ denotes a set of lightweight interconnected
_microservices_. A _source microservice_ $m_{src}$ processing the data stream
produced by the source dataset $\mathcal{D}$. A _sink microservice_ $m_{snk}$
representing its final results $\mathcal{R}$. $\varepsilon$ indicates a set of
_data streams_ $d_{u,i}$ flowing from an upstream microservice $m_{u}$ to a
downstream microservice $m_{i}\in\mathcal{M}$. $\mathcal{I}$ denotes a set of
_cloud infrastructure_ and $req$ is a set of user requirements.

Figure 3: An overview of the privacy-preserving preprocessing pipeline for
whole slide image in a distributed processing manner.

Along with these lines, the research problem has turned the emphasis of
studying on privacy-preserving service orchestration – Or more precisely, how
to customize a virtual infrastructure and schedule the workflow execution
under privacy-preserving constraints while reducing its time and monetary
cost?

Privacy requirements: Reconstruction of a WSI from distributed resource nodes
can lead to finding similar medical images using content-based image retrieval
and extrapolating possible patient information from other sources. Therefore,
using this distributed scheme, their privacy will be preserved during the
process.

Bi-objective optimization: Reducing the execution time of the application over
the cloud can be crucial for many stakeholders as it can lead to significant
cost savings and improve overall processing time for WSI. We aim to reduce the
monetary cost $f_{1}$ and minimize the application’s maximum completion time
(i.e., makespan) $f_{2}$. Let us denote $ET(m_{src})$ and $ET(m_{snk})$ as the
execution time of splitting and aggregation services over the trusted server.
$m_{det}$ for the artifact detection microservices, which will be deployed to
the cloud where $ET(m_{det}(A_{e,k},\mathcal{D}_{e,k}),\mathcal{I}_{k})$
denotes its total execution time. Then, the bi-objective optimization problem
can be formulated as follows:

$\displaystyle\min f_{1}$
$\displaystyle=\sum_{k=1}^{K}ET(m_{det}(A_{e,k},\mathcal{D}_{e,k}),\mathcal{I}_{k})\times
p_{k}\times x_{k}$ (2) $\displaystyle\min f_{2}$
$\displaystyle=ET(m_{src}(\mathcal{D}))+makespan(m_{det})+ET(m_{snk}(\mathcal{R}))$
(3)

Here $K$ and $p_{k}$ represent the number of the split data sets and the unit
price of cloud infrastructure $\mathcal{I}_{k}$, subject to,

$\small\begin{gathered}makespan(m_{det})=\max\\{ET(m_{det}(A_{e,k},\mathcal{D}_{e,k}),\mathcal{I}_{k})\times
x_{k}\\}\\\ ET(m_{det}(A_{e,k},\mathcal{D}_{e,k}),\mathcal{I}_{k})>0\\\
\texttt{where}\hskip 9.00002ptx_{k}=\begin{cases}1,&\text{if }m_{det}\text{ is
mapping to }\mathcal{I}_{k},\\\ 0,&\text{otherwise}.\end{cases}\end{gathered}$

##  V Discussion and Future Work

This work-in-progress paper presents the methodology for privacy-preserving
task-based parallel applications for distributed cloud environments. Our
method enables domain-specific users to handle gigapixel medical images
efficiently, maintaining privacy among distributed nodes. In future work, we
will develop prototypes and demonstrate the benefits of our pipeline using
datasets from different hospitals and integrating the method with a Jupyter-
based virtual research environment.

## Acknowledgment

This work has been funded by the European Union project CLARIFY (860627),
$\text{ENVRI}^{\text{FAIR}}$ (824068), BlueCloud-2026 (101094227) and
LifeWatch ERIC.

## References

  * [1] N. Kanwal, F. Pérez-Bueno, A. Schmidt, R. Molina, and K. Engan, “The devil is in the details: Whole slide image acquisition and processing for artifacts detection, color variation, and data augmentation. a review.” _IEEE Access_ , 2022.
  * [2] N. Kanwal, T. Eftestøl, F. Khoraminia, T. C. M. Zuiverloon, and K. Engan, “Vision transformers for small histological datasets learned through knowledge distillation,” in _Advances in Knowledge Discovery and Data Mining_. Springer Nature Switzerland, 2023, pp. 167–179.
  * [3] Z. Zhao, A. Belloum, and M. Bubak, “Special section on workflow systems and applications in e-Science,” _Future Generation Computer Systems_ , vol. 25, no. 5, pp. 525–527, May 2009. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0167739X08001726
  * [4] Y. Wang, S. Koulouzis, R. Bianchi, N. Li, Y. Shi, J. Timmermans, W. D. Kissling, and Z. Zhao, “Scaling notebooks as re-configurable cloud workflows,” _Data Intelligence_ , vol. 4, no. 2, pp. 409–425, 2022.
  * [5] Z. Zhao, S. Koulouzis, R. Bianchi, S. Farshidi, Z. Shi, R. Xin, Y. Wang, N. Li, Y. Shi, J. Timmermans _et al._ , “Notebook-as-a-vre (naavre): From private notebooks to a collaborative cloud virtual research environment,” _Software: Practice and Experience_ , 2022.
  * [6] N. Kanwal, S. Fuster, F. Khoraminia, T. C. Zuiverloon, C. Rong, and K. Engan, “Quantifying the effect of color processing on blood and damaged tissue detection in whole slide images,” in _2022 IEEE 14th Image, Video, and Multidimensional Signal Processing Workshop (IVMSP)_. IEEE, 2022, pp. 1–5.
  * [7] L. Hurwicz, “The design of mechanisms for resource allocation,” _The American Economic Review_ , vol. 63, no. 2, pp. 1–30, 1973.
  * [8] A. Taal, J. Wang, C. de Laat, and Z. Zhao, “Profiling the scheduling decisions for handling critical paths in deadline-constrained cloud workflows,” _Future Generation Computer Systems_ , vol. 100, pp. 237–249, 2019.
  * [9] Y. Wang, H. Liu, W. Zheng, Y. Xia, Y. Li, P. Chen, K. Guo, and H. Xie, “Multi-objective workflow scheduling with deep-q-network-based multi-agent reinforcement learning,” _IEEE access_ , vol. 7, pp. 39 974–39 982, 2019\.
  * [10] E. Soljanin, “Technical perspective: Balancing at all loads,” _Communications of the ACM_ , vol. 65, no. 5, pp. 110–110, 2022.
  * [11] H. Kargupta, S. Datta, Q. Wang, and K. Sivakumar, “Random-data perturbation techniques and privacy-preserving data mining,” _Knowledge and Information Systems_ , vol. 7, pp. 387–414, 2005.
  * [12] Q. Li, J. S. Gundersen, R. Heusdens, and M. G. Christensen, “Privacy-preserving distributed processing: metrics, bounds and algorithms,” _IEEE Transactions on Information Forensics and Security_ , vol. 16, pp. 2090–2103, 2021.