# Optically-biased Rydberg microwave receiver enabled by hybrid nonlinear
interferometry

Sebastian Borówka  Centre for Quantum Optical Technologies, Centre of New
Technologies, University of Warsaw, Banacha 2c, 02-097 Warsaw, Poland  Faculty
of Physics, University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland  Mateusz
Mazelanik  Centre for Quantum Optical Technologies, Centre of New
Technologies, University of Warsaw, Banacha 2c, 02-097 Warsaw, Poland
Wojciech Wasilewski  Centre for Quantum Optical Technologies, Centre of New
Technologies, University of Warsaw, Banacha 2c, 02-097 Warsaw, Poland  Faculty
of Physics, University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland  Michał
Parniak<EMAIL_ADDRESS>Centre for Quantum Optical Technologies, Centre
of New Technologies, University of Warsaw, Banacha 2c, 02-097 Warsaw, Poland
Faculty of Physics, University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland

###### Abstract

The rapid development in Rydberg-atomic detection of microwave fields has
paved the way to a new type of intrinsically calibrated RF measurements. The
coupling of Rydberg vapors medium both to microwave and optical fields allows
harnessing the merits of all-optical detection, e.g. weak disruption of the
measured field and invulnerability to extremely strong fields, owing to the
lack of conventional RF antenna in the detector. The trouble with this
approach arises, as the greatest sensitivity is exhibited with the use of
additional microwave field acting as a local oscillator, and the measurement
can no longer be all-optical. Here we propose a different solution, optical-
bias detection, that allows truly all-optical operation, while retaining most
of the sensitivity. We tackle the issue of uncorrelated laser phase, emerging
in this type of detection, by facilitating separate measurement of the phase
spectrum in a nonlinear process and using the results in post-processing. We
report the sensitivity of $176\ \mathrm{nV/cm/\sqrt{Hz}}$ and reliable
operation up to $3.5\ \mathrm{mV/cm}$ of $13.9\ \mathrm{GHz}$ electric field,
while comparing with the state of the art (although not all-optical) method
realized in the same setup.

## Introduction

While the exceptional sensitivity of Rydberg atoms to microwave (MW) fields
has been known from atomic beam and cold-atom experiments, only recent
advances enabled sensing using hot-atom vapor cell systems. The original
approaches to hot-atom sensors were based on the splitting of an
electromagnetically induced transparency (EIT) line due to MW-induced Autler-
Townes (A-T) effect. This elegant and simple approach allows for a limited
sensitivity [1, 2, 3] and enables self-calibration based on atomic constants,
yet does not yield phase information. Since then, invaluable progress has been
made, extending this method to provide phase information and enhance sensing
[4, 5, 6, 7, 8, 9] or enable applications in communication protocols [10, 11,
12, 13, 14, 15, 16, 17, 18], as well as miniaturization of the receivers [19,
20].

The sensitivity metric has been revolutionized by treating the atoms as a MW
mixer with optical output [21, 22], effectively constituting a superheterodyne
receiver (superhet) [23]. Further developments allowed consecutive
enhancements of the method [24, 25, 26, 27, 28] and adapting it to specific
problems, such as measurements of angle of arrival [29] and polarization [30].

Atomic superhet reception provides an important link to standard radio
communication, where superheterodyne detection method is present in all modern
receivers. In the end, however, while the readout is optical, the Rydberg-atom
receiver still requires a strong MW local oscillator (LO) for operation. In
terms of practical applications, this solution is not ideal, as a MW antenna
needs to be a part of the receiver, rendering useless some of the advantages
of atomic sensors, such as weakly-disruptive (stealthy) measurement of
electric field. [31]

In this paper we present a different type of Rydberg-atom detection, based on
optically-biasing ("optical-bias") the mixing process, thus erasing the
requirement of a MW LO and making the receiver all-optical. Drawing
inspiration from our recent discoveries in MW-to-optical conversion in Rydberg
atoms [32], we extend the well-established two-photon Rydberg excitation
scheme by the addition of two optical fields in the near-infrared, realising
three-photon Rydberg excitation scheme via practically Doppler-free
$5\mathrm{D}_{5/2}$ energy level (in rubidium), already explored in some
earlier works [33, 34]. Effectively the optical fields play the same role as
the MW LO in superhet detection – the optical-bias induces beat modulation in
the probe field transmission, and the measurement of its amplitude serves as a
MW electrometer. Similar approaches have been presented with the use of two RF
fields [35, 36], where several all-optical solutions were proposed as well
[36].

The most important challenge arising in this approach is the collective laser
phase spectrum, which transfers into the modulation signal. In the superhet
detection this issue is mitigated by the addition of phase-stable MW LO, which
enables phase stabilization for the optical fields, so long as they are
reasonably near atomic two-photon resonance. Similar solutions have been
presented in cases of two RF transitions [35, 36], however, it is no longer
possible in all-optical scheme. We tackle the problem of laser phase by
measuring it in a separate process. To get the laser phase reference we employ
a difference frequency generation, using the same laser fields. With real-time
correlation performed on an field programmable gate array (FPGA) circuit
between measured signal and the reference, we are able to remove the
introduced phase fluctuations. This way we render the signal laser-noise-free
and achieve performance directly comparable to the superhet detection.

Figure 1:  A, Energy level structure utilized in the optical-bias detection.
Two photon (_probe_ – _coupling_) and three photon (_probe_ – _dressing_ –
_coupling_) Rydberg excitation paths are used to access both energy levels
connected by the MW transition. The $\sigma^{+}$ transitions ensure the
largest transition dipole moments. All of the optical fields are atomic
resonant, apart from the indicated detuning $\Delta_{5\mathrm{D}}=-1.8\
\mathrm{MHz}$. The beat modulation in probe transmission is also observed at
$\Delta=1.8\ \mathrm{MHz}$ due to the detuning of MW field. B, Experimental
setup of the optical-bias receiver. Three optical ($480\ \mathrm{nm}$, $776\
\mathrm{nm}$ and $1258\ \mathrm{nm}$) fields are divided into two paths. In
the first path they counterpropagate with respect to the $780\ \mathrm{nm}$
probe field, enabling partial Doppler effect cancellation. The fields,
propagating as Gaussian beams have circular polarizations and are focused to
waists of around $w_{0}=250\ \mathrm{\mu m}$, and combined with dichroic
mirrors and spectral filters inside ${}^{85}\mathrm{Rb}$ vapor cell. The
second path leads to frequency generation (DFG) setup for laser phase spectrum
detection. The $480\ \mathrm{nm}$ and $776\ \mathrm{nm}$ fields induce DFG at
$1258\ \mathrm{nm}$ shifted by the frequency of the detected MW field, as
dictated by the conservation of energy. Combined with the non-shifted $1258\
\mathrm{nm}$ field used in the detection setup, they induce beat modulation at
$13.9\ \mathrm{GHz}$ on a photodetector. Then downmixing with LO signal
enables the retrieval of laser-noise spectrum shifted to lower frequencies.
For measurements, the cell is placed inside a MW absorbing shield with a
helical MW antenna acting as a signal source. In the case of superhet
detection measured for comparison, the $776\ \mathrm{nm}$ and $1258\
\mathrm{nm}$ are switched off and MW LO is combined with the signal at the
antenna. C, Operating principle of a MW optical-bias receiver. An additional
measurement of the phase spectrum of optical fields allows the phase
compensation. D, Operating principle of a MW superhet receiver. The MW LO
allows the phase compensation of optical fields. Figure 2:  A, Spectrum of the
probe field modulation obtained in optical-bias detection of $720\ \mathrm{\mu
V/cm}$ MW field. The maximum of the signal is $25\ \mathrm{dB}$ above noise
level and the spectral width is FWHM. (The visible chime is a small superhet-
type signal resulting from cross-talk. It is outside the set detection
bandwidth.) B Respective spectrum obtained via DFG and beat modulation, as the
reference. The maximum of the signal is $41\ \mathrm{dB}$ above noise level.
C, Respective spectrum of the phase-correlated signal. The maximum of the
signal is $60\ \mathrm{dB}$ above noise level. The spectral width of the
signal is Fourier-limited at $10\ \mathrm{Hz}$. The artifact spurs are at
$-38\ \mathrm{dB}$ below the signal.

## Results

### Comparison between superhet and optical-bias MW detection setups

In the superhet detection, pictured schematically in the Fig. 1D, the optical
fields, typically probe and coupling beams in counterpropagation enable
precise detection of the power of the MW field. Like any fast power detector
this setup can be used as a mixer to detect a weak detuned MW signal using a
strong resonant MW LO. The total MW intensity beats at the frequency of the
detuning leading to a modulation in the probe field transmission. The
modulation can be then analyzed spectrally to constitute frequency-resolved
detection of weak MW electric fields.

On the contrary, in the optical-bias detection, pictured schematically in the
Fig. 1C, the weak MW signal is converted to weak optical signal, a sideband of
the probe signal, which is subjected to optical heterodyne detection. Thus,
the detection stage is conveniently decoupled from the volume containing MW
field.

Alternative picture involves analyzing MW polarization of the atoms (bias)
induced by the laser fields. In this picture, the beating occurs between said
polarization and incoming weak MW field. It causes modulation of atomic
coherences and modulation of the probe field transmission.

The rubidium energy levels employed in the optical-bias detection are shown in
the Fig. 1A. The standard two-photon Rydberg excitation at $780\ \mathrm{nm}$
and $480\ \mathrm{nm}$ is supplemented with two fields in near-infrared, $776\
\mathrm{nm}$ and $1258\ \mathrm{nm}$, both of which are convenient to work
with in terms of fiber-based solutions. All of the fields are resonant, apart
from the $\Delta_{5\mathrm{D}}$ detuning from the $5\mathrm{D}_{5/2}$ and the
beat note detuning, $\Delta$. We found the optimal working point when
$\Delta_{5\mathrm{D}}=-\Delta$, as indicated in the Fig. 1C, and set this
value to $1.8\ \mathrm{MHz}$.

### Phase reference in superhet and optical-bias detection

Both superhet and optical-bias detection are forms of mixing of the incoming
weak MW signal with either supplied MW LO or optical fields performing the
role of the LO in optical-bias detection. A mixer detects only the amplitude
quadrature of the incoming signal, i.e. only the component in phase with the
LO. In case of superhet detection the detection phase is fixed by the
reference MW LO. In case of optical-bias detection, the phase information is
provided with the polarization (bias) of atoms. This is a result of preparing
the atoms with optical fields and therefore every field involved in the
process add to the bias effect. There, a relative phase of involved fields
$\phi$ plays a crucial role. In our case, this phase can be expressed as

$\phi=\phi_{\mathrm{480}}-\phi_{\mathrm{1260}}-\phi_{\mathrm{776}},$ (1)

where in the subscripts we denote respective optical fields’ wavelengths (in
nm). Note that the contribution of $\phi_{\mathrm{780}}$ cancels out.

In general, each component of this phase is time-dependent and may be
separated into fast-varying (at the optical frequencies) and slow-varying (due
to the optical fields phase instabilities). In the interaction picture, via
considering rotating wave approximations in relation to each transition
between the detuned levels, pictured in the Fig. 1A, only the slow-varying
elements remain. These can be then combined into time-dependent slow-varying
phase, $\phi(t)$. In our approach, instead of minimizing this value (which may
be difficult in practice, as one would need to use extremely well stabilized
lasers), we decide to register it in a separate process and refer to it in
post-processing.

To register the $\phi(t)$ phase several solutions are possible, including the
use of a separate atomic vapor cell with MW LO. In our implementation, we
decided to take advantage of an optical nonlinear process. As shown in the
Fig. 1C, we realize a difference frequency generation (DFG) between $480\
\mathrm{nm}$ and $776\ \mathrm{nm}$ field and combine the DFG signal with
$1258\ \mathrm{nm}$ laser field to detect the beat modulation around $13.9\
\mathrm{GHz}$. The beat modulation is then downconverted to $4.6\
\mathrm{MHz}$ and processed electronically on FPGA circuit to yield $\phi(t)$
measurement in real time.

### Spectral characteristics

An example (detection of $720\ \mathrm{\mu V/cm}$ MW field) of uncompensated
signal in the optical-bias method is presented in the Fig. 2A in relation to
noise, which is shot noise of the probe optical field. The spectral FWHM (full
width at half maximum) of the signal is $62\ \mathrm{kHz}$. The uncompensated
spectral signal shape for monochromatic MW is dictated by collective phase
stability of the lasers, thus corresponding to the $\phi(t)$ phase reference
spectrum pictured in the Fig. 2B.

The phase-compensated signal of the example is presented in the Fig. 2C. The
measurement duration chosen in this case and in the onward analysis is $t=100\
\mathrm{ms}$, which corresponds to the resolution of $10\ \mathrm{Hz}$. Thus
the noise-compensated signal is Fourier-limited to this resolution. The
improvement in the signal in relation to noise, from $25\ \mathrm{dB}$ to $60\
\mathrm{dB}$, yielding $35\ \mathrm{dB}$, is consistent with the spectrum-
based estimation limit, which is $(62\ \mathrm{kHz})/(10\ \mathrm{Hz})=38\
\mathrm{dB}$.

The artifact spurs at the level of $-38\ \mathrm{dB}$ below the signal, are
the results of imperfect balancing of MW detection and noise detection setups
and can be eliminated with better alignment of the hybrid interferometer.

Note that in the Fig. 2C we present only a part of instantaneously acquired MW
spectrum. However in full, we are able to observe real-time measurement of $2\
\mathrm{MHz}$ MW bandwidth with $1\ \mathrm{Hz}$ resolution (for time base
$t=1\ \mathrm{s}$). To further improve resolution and bandwidth, better
electronic solutions need to be employed, and the hybrid interferometer has to
be stabilized, e.g. with shortening interferometric arms’ lengths of the
system or using fiber solutions.

### Dynamic range

With attenuating MW signal we measure the sensitivity of the optical-bias
detection method. The results are presented in the Fig. 3, where we facilitate
comparison with the superhet method. To achieve adequate comparison, we
normalize the signals to respective noise floors, which in both methods are
shot noise. In the optical-bias method we find saturation due to energy level
shifts at MW field $3.5\ \mathrm{mV/cm}$, while the sensitivity reaches $176\
\mathrm{nV/cm/\sqrt{Hz}}$. These results are compared with the superhet method
realized in the same setup, where we register the saturation at $1.84\
\mathrm{mV/cm}$ and the sensitivity of $87\ \mathrm{nV/cm/\sqrt{Hz}}$.
Notably, the optical-bias method saturates for larger MW fields, while having
lower sensitivity, though the dynamic range is almost equal at $76\
\mathrm{dB}$ for $t=100\ \mathrm{ms}$ measurement duration.

The behavior of atomic saturation is predicted by the numerical simulation
based on the theory model described in the Methods and Materials section. The
results comparing the experimental data with the simulated function are
presented in the Fig. 3. In this case the shape of the function is predicted
with measured parameters and the only free parameter is the constant scaling
the overall signal strength.

Figure 3: Comparison between results obtained in superhet measurement (blue
dots) and optical-bias (red dots). Both results are presented in relation to
their respective noise levels (black dashed line), which in both cases is shot
noise of the detected probe field transmission. The optical-bias measurement
method results in comparable, though smaller sensitivity and overall
efficiency. Notably however, it becomes saturated for larger MW fields than
the superhet method, thus retaining very similar dynamic range. The results
presented here are averaged over $n=8$ shots for each point to facilitate
better comparison. The yellow solid line represents a theoretical prediction.
Notably, the saturation point is predicted accurately and the experimental
results differ from the theoretical predictions only in the stronger MW field
regime, where the atomic response can no longer be considered instantaneous.

### Beat modulation transfer

We then perform an experimental analysis of the transfer of beat modulation in
probe field detuning domain, once again in the working point of $720\
\mathrm{\mu V/cm}$ MW signal field. We start with the superhet method, where
in the probe field absorption, Fig. 4A, we observe widened EIT feature (we
found the optimal working point of MW LO to be at only $1.6\ \mathrm{MHz}$
Rabi frequency for the beat modulation at $1.8\ \mathrm{MHz}$). Respectively,
the modulation transfer dependence is presented in the Fig. 4D, bearing
resemblance to the transfer of amplitude modulation of MW field [9].

The optical-bias method exhibits an EIT feature resembling off-resonant A-T
splitting due to MW field, although this effect is only due to optical fields,
Fig. 4B. The visible EIT outside the $\pm 10\ \mathrm{MHz}$ range is due to
interaction with the sublevels not taking part in the optical-bias detection
method, particularly the hyperfine splitting of the $5\mathrm{P}_{3/2}$ and
$5\mathrm{D}_{5/2}$ levels. The modulation transfer, Fig. 4E, is similar to
the superhet method, when the amplitude of the modulation is considered.
However, the phase dependence on probe detuning differs in the direction of
transitions.

The results obtained in the numerical simulation based on the theoretical
model are shown in the Fig. 4C,F. Both the probe field transmission and the
amplitude of modulation are reproduced in terms of shapes in the probe field
detuning domain. The difference in the direction of the transition of the
phase may be attributed to the interaction with other atomic sublevels visible
in the experimental probe transmission, Fig. 4B.

Figure 4: Comparison of EIT effects (upper row) and modulation transfer (lower
row), both amplitude $|S|$ (red lines) and phase $\arg(S)$ (blue lines)
between superhet (left column) and optical-bias (middle column) detection
results, and the theoretical prediction for optical-bias (right column) in the
domain of probe field detuning. The MW field in all cases is $720\ \mathrm{\mu
V/cm}$. Notably, despite similarities in the shape of modulation transfer
amplitude, the optical-bias method exhibits a different direction of the
transition of phase in the demodulated signal than the superhet method. The
opacity of blue lines representing phases $\arg(S)$ is determined by the
respective amplitude $|S|$ to signify measurement uncertainty.

## Discussion

We demonstrated a method of Rydberg atomic detection of MW fields that is both
all-optical and sensitive, as shown by the comparison with the standard
superhet method. This paves a way to achieving all of the merits of all-
optical measurement, i.e. weak disruption of the field and invulnerability to
strong interaction, while not compromising on the senstivity. Furthermore, as
the atomic medium interacts with a single MW frequency, it is possible to
precisely design a resonant MW cavity around it, which would allow for
surpassing the conventional MW receivers [37, 38].

In principle, the demonstrated optical-bias method can be extended to
different Rydberg transitions, thus enabling the tunable detection of wide
bandwidth of MW and mm-wave frequencies. However, the limitation of the
presented realization is the presence of the LO frequency needed to downmix
the phase reference signal, which may be difficult for high frequency. This
issue can be resolved e.g. with using another vapor cell with LO signal as a
reference instead of the DFG mixing demonstrated in this work. Alternatively,
one could use an electro-optic modulator with high order sideband modulation
to shift one of the laser frequencies in the phase reference setup.

On the other hand, the optical-bias detection method can be massively improved
with the use of phase-stable laser fields, e.g. with locking all of the lasers
to an optical frequency comb. In this case, if the desired sensitivity and
spectral resolution is achieved, the optical-bias would work without the need
for additional $\phi(t)$ phase referencing in a hybrid interferometer.

## Methods and Materials

### Theoretical model

To facilitate the comparison with the theoretical predictions, we prepare a
numerical simulation based on the following considerations. The atomic state
time evolution is described by Gorini–Kossakowski–Sudarshan–Lindblad (GKSL)
equation,

$\partial_{t}\hat{\rho}=\frac{1}{i\hbar}[\hat{H},\hat{\rho}]+\mathcal{L}[\hat{\rho}],$
(2)

where $\hat{H}$ is Hamiltonian and $\mathcal{L}[{\cdot}]$ is a superoperator
responsible for sources of decoherence, such as spontaneous emission. We
consider the steady state solution to the GKSL equation,
$\partial_{t}\hat{\rho}(t)=0$. Atomic numerical data, such as state lifetimes
and transition dipole moments, is found via Alkali Rydberg Calculator [39]. To
get full agreement with the experiment these operations have to be done for a
range of velocity classes present in a room-temperature atomic medium, as well
as for a range of laser field power, changing with the radial position in the
Gaussian beams.

### Details of optical-bias setup

The MW detection part of the optical-bias setup employed in the experiment is
pictured in the Fig. 1C. Four optical Gaussian beams are all focused with
Gaussian waists $w_{0}=250\ \mathrm{\mu m}$ inside a rubidium vapor cell. The
rubidium is in natural abundance isotope proportion, though only
${}^{85}\mathrm{Rb}$ is addressed in this work. The matched optimized Rabi
frequencies of the fields (at the beam centers) are respectively $5.5\
\mathrm{MHz}$, $7.5\ \mathrm{MHz}$, $6.2\ \mathrm{MHz}$ and $9.5\
\mathrm{MHz}$ for $780\ \mathrm{nm}$, $480\ \mathrm{nm}$, $776\ \mathrm{nm}$
and $1258\ \mathrm{nm}$ fields. All of the fields have matched circular
polarizations. The lasers are stabilized via either cavity transfer locks or
optical phase-locked loop to a common reference source (a narrowband fiber
laser at $1560\ \mathrm{nm}$, frequency-doubled to $780\ \mathrm{nm}$ and
referenced to Rb cell with modulation-transfer lock). This results in
collective spectral stability of around $62\ \mathrm{kHz}$, as seen in the
Fig. 2A,B.

Near the vapor cell acting as a detector, we placed a MW helical antenna
acting as a source of weak MW field (and additionally also as a source of MW
LO in the superhet measurements facilitated for comparisons). The antenna
emits circularly-polarized MW field colinearly with the optical fields. Both
the antenna and the vapor cell are placed inside a MW absorbing shield (made
from LeaderTech EA-LF500 material), providing protection from MW reflections
and on-air noise. The probe laser beam counterpropagates to all the other
fields and its power is registered on an avalanche photodiode (Thorlabs
APD430A). This detection is shot-noise limited with respect to the probe field
power.

### Details of laser phase detection setup

The optical beams $480\ \mathrm{nm}$ and $776\ \mathrm{nm}$ of respective
powers ${\sim}100\ \mathrm{mW}$ and ${\sim}15\ \mathrm{mW}$ are focused and
Rayleigh length-matched in a z-cut $30\ \mathrm{mm}$ MgO:PPLN crystal with
$5.17\ \mathrm{\mu m}$ poling period, at ${\sim}85\ {}^{\circ}\mathrm{C}$,
resulting in $1258\ \mathrm{nm}$ DFG signal at ${\sim}25\ \mathrm{\mu W}$. We
combine the DFG signal with ${\sim}1\ \mathrm{mW}$ of sampled $1258\
\mathrm{nm}$ laser field (which is also propagated through the crystal to
reduce the differences in optical paths). Both of $1258\ \mathrm{nm}$ fields
are then fiber coupled. The beat modulation between them is detected with a
fast photodiode (25G SFP28 module for $1310\ \mathrm{nm}$) and serves as an
electronic MW phase reference signal.

The phase reference from DFG photodiode is downmixed with additional MW LO to
a frequency of $4.6\ \mathrm{MHz}$ convenient for Red Pitaya STEMlab 125-14
ADC/FPGA board used as a measurement tool. The signal is digitized and
digitally demodulated, thus providing a complex signal encoding reference
phase $\phi(t)$. The reference phase is subtracted from the demodulated APD
signal yielding a phase corrected signal of optical-bias detection method.
Spectra of uncompensated signal, reference phase signal and compensated signal
are presented in the Fig. 2.

### Details of the MW setup and calibration

The MW frequency resonant with the $55^{2}\mathrm{D}_{5/2}\rightarrow
54^{2}\mathrm{F}_{7/2}$ transition is measured as $\omega_{0}=13912.4\
\mathrm{MHz}$ via A-T splitting. The MW signal is then set at
$\omega_{0}+\Delta$, where $\Delta=1.8\ \mathrm{MHz}$. This offset frequency
is chosen large enough to avoid any electronic problems at low frequencies and
fit withing detection bandwidth, resulting in a well resolved sideband on the
probe field. The MW signal power is attenuated with a programmable RF
attenuator (MiniCircuits RCDAT-18G-63). All MW fields are calibrated using A-T
splitting method in the resolvable regime.

In the case of superhet detection, facilitated for comparison, the MW LO is
added with a power splitter inserted in the antenna cable. The MW LO is
resonant with the atoms and at $800\ \mathrm{\mu V/cm}$ ($1.6\ \mathrm{MHz}$
Rabi frequency), which is measured experimentally as optimal in terms of SNR
for the set detection parameters, and then combined with the MW signal at the
antenna.

All the MW signals are generated with separate LMX2529 PLL frequency
synthesizers synchronized to the same clock reference. The signal from APD is
measured with the STEMlab 125-14 and digitally demodulated yielding a complex
amplitude and phase signal.

## Data availability

Data underlying the results presented in this paper are available in the Ref.
[40].

## Code availability

The codes used for the numerical simulation and the analysis of experimental
data are available from the authors upon request.

## Author contributions

S.B. built the optical, microwave and electronic processing setups. W.W.
developed the measurement system firmware and software, assisted by S.B. who
programmed the software. S.B. took the measurements and analyzed the data
assisted by other authors. S.B, M.M. and M.P. prepared the figures and the
manuscript. W.W. developed the theory and numerical simulation assisted by
S.B. who facilitated comparison with experiments. M.P. conceived the central
idea and led the project assisted by M.M. and W.W.

###### Acknowledgements.

We thank K. Banaszek for the generous support. The “Quantum Optical
Technologies” (MAB/2018/4) project is carried out within the International
Research Agendas programme of the Foundation for Polish Science co-financed by
the European Union under the European Regional Development Fund. This research
was funded in whole or in part by National Science Centre, Poland grant no.
2021/43/D/ST2/03114.

## References

  * Sedlacek  _et al._ [2012] J. A. Sedlacek, A. Schwettmann, H. Kübler, R. Löw, T. Pfau, and J. P. Shaffer, Nature Physics 8, 819 (2012).
  * Sedlacek  _et al._ [2013] J. A. Sedlacek, A. Schwettmann, H. Kübler, and J. P. Shaffer, Physical Review Letters 111, 063001 (2013).
  * Fan  _et al._ [2015] H. Fan, S. Kumar, J. Sedlacek, H. Kübler, S. Karimkashi, and J. P. Shaffer, Journal of Physics B: Atomic, Molecular and Optical Physics 48, 202001 (2015).
  * Kumar  _et al._ [2017] S. Kumar, H. Fan, H. Kübler, J. Sheng, and J. P. Shaffer, Scientific Reports 7, 42981 (2017).
  * Anderson  _et al._ [2018] D. A. Anderson, E. G. Paradis, and G. Raithel, Applied Physics Letters 113, 073501 (2018).
  * Meyer  _et al._ [2021] D. H. Meyer, P. D. Kunz, and K. C. Cox, Physical Review Applied 15, 014053 (2021).
  * Chopinaud and Pritchard [2021] A. Chopinaud and J. Pritchard, Physical Review Applied 16, 024008 (2021).
  * Liu  _et al._ [2021] X. Liu, F. Jia, H. Zhang, J. Mei, Y. Yu, W. Liang, J. Zhang, F. Xie, and Z. Zhong, AIP Advances 11, 085127 (2021).
  * Borówka  _et al._ [2022] S. Borówka, U. Pylypenko, M. Mazelanik, and M. Parniak, Applied Optics 61, 8806 (2022).
  * Deb and Kjærgaard [2018] A. B. Deb and N. Kjærgaard, Applied Physics Letters 112, 211106 (2018).
  * Meyer  _et al._ [2018] D. H. Meyer, K. C. Cox, F. K. Fatemi, and P. D. Kunz, Applied Physics Letters 112, 211108 (2018).
  * Cox  _et al._ [2018] K. C. Cox, D. H. Meyer, F. K. Fatemi, and P. D. Kunz, Physical Review Letters 121, 110502 (2018).
  * Song  _et al._ [2019] Z. Song, H. Liu, X. Liu, W. Zhang, H. Zou, J. Zhang, and J. Qu, Optics Express 27, 8848 (2019).
  * Holloway  _et al._ [2019] C. L. Holloway, M. T. Simons, J. A. Gordon, and D. Novotny, IEEE Antennas and Wireless Propagation Letters 18, 1853 (2019).
  * Jiao  _et al._ [2019] Y. Jiao, X. Han, J. Fan, G. Raithel, J. Zhao, and S. Jia, Applied Physics Express 12, 126002 (2019).
  * Anderson  _et al._ [2021] D. A. Anderson, R. E. Sapiro, and G. Raithel, IEEE Transactions on Antennas and Propagation 69, 2455 (2021).
  * Otto  _et al._ [2021] J. S. Otto, M. K. Hunter, N. Kjærgaard, and A. B. Deb, Journal of Applied Physics 129, 154503 (2021).
  * Elgee  _et al._ [2023] P. K. Elgee, J. C. Hill, K.-J. E. LeBlanc, G. D. Ko, P. D. Kunz, D. H. Meyer, and K. C. Cox, Applied Physics Letters 123, 084001 (2023).
  * Simons  _et al._ [2018] M. T. Simons, J. A. Gordon, and C. L. Holloway, Applied Optics 57, 6456 (2018).
  * Simons  _et al._ [2019a] M. T. Simons, A. H. Haddab, J. A. Gordon, D. Novotny, and C. L. Holloway, IEEE Access 7, 164975 (2019a).
  * Simons  _et al._ [2019b] M. T. Simons, A. H. Haddab, J. A. Gordon, and C. L. Holloway, Applied Physics Letters 114, 114101 (2019b).
  * Gordon  _et al._ [2019] J. A. Gordon, M. T. Simons, A. H. Haddab, and C. L. Holloway, AIP Advances 9, 045030 (2019).
  * Jing  _et al._ [2020] M. Jing, Y. Hu, J. Ma, H. Zhang, L. Zhang, L. Xiao, and S. Jia, Nature Physics 16, 911 (2020).
  * Jia  _et al._ [2021] F.-D. Jia, H.-Y. Zhang, X.-B. Liu, J. Mei, Y.-H. Yu, Z.-Q. Lin, H.-Y. Dong, Y. Liu, J. Zhang, F. Xie, and Z.-P. Zhong, Journal of Physics B: Atomic, Molecular and Optical Physics 54, 165501 (2021).
  * Prajapati  _et al._ [2021] N. Prajapati, A. K. Robinson, S. Berweger, M. T. Simons, A. B. Artusio-Glimpse, and C. L. Holloway, Applied Physics Letters 119, 214001 (2021).
  * Li  _et al._ [2022a] H. Li, J. Hu, J. Bai, M. Shi, Y. Jiao, J. Zhao, and S. Jia, Optics Express 30, 13522 (2022a).
  * Hu  _et al._ [2022] J. Hu, H. Li, R. Song, J. Bai, Y. Jiao, J. Zhao, and S. Jia, Applied Physics Letters 121, 014002 (2022).
  * Cai  _et al._ [2023] M. Cai, S. You, S. Zhang, Z. Xu, and H. Liu, Applied Physics Letters 122, 161103 (2023).
  * Robinson  _et al._ [2021] A. K. Robinson, N. Prajapati, D. Senic, M. T. Simons, and C. L. Holloway, Applied Physics Letters 118, 114001 (2021).
  * Wang  _et al._ [2023] Y. Wang, F. Jia, J. Hao, Y. Cui, F. Zhou, X. Liu, J. Mei, Y. Yu, Y. Liu, J. Zhang, F. Xie, and Z. Zhong, Optics Express 31, 10449 (2023).
  * Yuan  _et al._ [2023] J. Yuan, W. Yang, M. Jing, H. Zhang, Y. Jiao, W. Li, L. Zhang, L. Xiao, and S. Jia, Reports on Progress in Physics 86, 106001 (2023).
  * Borówka  _et al._ [2023] S. Borówka, U. Pylypenko, M. Mazelanik, and M. Parniak, Continuous wideband microwave-to-optical converter based on room-temperature rydberg atoms (2023).
  * Thaicharoen  _et al._ [2019] N. Thaicharoen, K. R. Moore, D. A. Anderson, R. C. Powel, E. Peterson, and G. Raithel, Physical Review A 100, 063427 (2019).
  * Li  _et al._ [2022b] W. Li, J. Du, M. Lam, and W. Li, Optics Letters 47, 4399 (2022b).
  * Anderson  _et al._ [2022] D. Anderson, R. Sapiro, L. Gonçalves, R. Cardman, and G. Raithel, Physical Review Applied 17, 044020 (2022).
  * Berweger  _et al._ [2023] S. Berweger, A. B. Artusio-Glimpse, A. P. Rotunno, N. Prajapati, J. D. Christesen, K. R. Moore, M. T. Simons, and C. L. Holloway, Physical Review Applied 20, 054009 (2023).
  * Santamaria-Botello  _et al._ [2022] G. Santamaria-Botello, S. Verploegh, E. Bottomley, and Z. Popovic, Comparison of noise temperature of rydberg-atom and electronic microwave receivers (2022).
  * Sandidge  _et al._ [2023] G. Sandidge, G. S. Botello, and Z. Popović, in  _2023 IEEE/MTT-S International Microwave Symposium - IMS 2023_ (IEEE, 2023).
  * Šibalić  _et al._ [2017] N. Šibalić, J. Pritchard, C. Adams, and K. Weatherill, Computer Physics Communications 220, 319–331 (2017).
  * Borówka  _et al._ [2024] S. Borówka, M. Mazelanik, W. Wasilewski, and M. Parniak, Replication data for: Optically-biased rydberg microwave receiver enabled by hybrid nonlinear interferometry, Harvard Dataverse (2024).


# On the Formation and Interaction of Multiple Supermassive Stars in
Cosmological Flows

Tyrone E. Woods National Research Council of Canada, Herzberg Astronomy &
Astrophysics Research Centre, 5071 West Saanich Road, Victoria, BC V9E 2E7,
Canada  Monash Centre for Astrophysics, School of Physics and Astronomy,
Monash University, VIC 3800, Australia  Samuel Patrick Institute of Cosmology
and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK  Daniel J.
Whalen Institute of Cosmology and Gravitation, University of Portsmouth,
Portsmouth PO1 3FX, UK  Alexander Heger Monash Centre for Astrophysics, School
of Physics and Astronomy, Monash University, VIC 3800, Australia  ARC Centre
of Excellence for Gravitational Wave Discovery (OzGrav), Melbourne, Australia
ARC Centre of Excellence for Astrophysics in Three Dimensions (ASTRO-3D),
Australia  Joint Institute for Nuclear Astrophysics, 1 Cyclotron Laboratory,
National Superconducting Cyclotron Laboratory, Michigan State University, East
Lansing, MI 48824-1321, USA

###### Abstract

Supermassive primordial stars with masses exceeding
$\mathord{\sim}10^{5}\,{\mathrm{M}_{\odot}}$ that form in atomically cooled
halos are the leading candidates for the origin of high-redshift quasars with
$z>6$. Recent numerical simulations, however, find that multiple accretion
disks can form within a halo, each of which can host a supermassive star.
Tidal interactions between the disks can gravitationally torque gas onto their
respective stars and alter their evolution. Later, when two satellite disks
collide, the two stars can come into close proximity. This may induce
additional mass exchange between them. We investigate the co-evolution of
supermassive stars in atomically-cooled halos driven by gravitational
interactions between their disks. We find a remarkable diversity of
evolutionary outcomes. The results depend on these interactions and how the
formation and collapse times of the stars in the two disks are correlated.
They range from co-evolution as main sequence stars to main sequence – black
hole pairs and black hole – black hole mergers. We examine the evolution of
these secondary supermassive stars in detail and discuss the prospects for
binary interactions on much smaller scales after the disks merge within their
host halos.

quasars: general — black hole physics — early universe — dark ages,
reionization, first stars — galaxies: formation — galaxies: high-redshift

##  1 Introduction

More than 200 quasars powered by supermassive black holes have now been
discovered at $z>6$ (Bosman, 2021), including nine at $z>7$ (Mortlock et al.,
2011; Wu et al., 2015; Bañados et al., 2018; Wang et al., 2021). A natural
explanation for the seeds of these quasars are direct collapse black holes
(DCBHs), which arise due to the formation of supermassive stars (SMSs) in
atomically cooled halos at $z\gtrsim 20$ (Bromm & Loeb, 2003; Lodato &
Natarajan, 2006; Wise et al., 2008; Regan & Haehnelt, 2009; Latif et al.,
2013). In this picture, a primordial halo grows to
$\mathord{\sim}10^{7}\,{\mathrm{M}_{\odot}}$ without forming stars because,
for example, it is immersed in strong Lyman-Werner (LW) UV fluxes (Dijkstra et
al., 2008; Agarwal et al., 2012; Latif et al., 2014) or highly supersonic
streaming motions of gas relative to dark matter (Tanaka & Li, 2014; Hirano et
al., 2017; Schauer et al., 2017), or it is dynamically heated by violent
mergers (Yoshida et al., 2003; Fernandez et al., 2014; Wise et al., 2019;
Regan et al., 2020). At this mass, virial temperatures reach
$\mathord{\sim}10^{4}\,\mathrm{K}$ that activate line cooling in hydrogen that
triggers rapid baryon collapse at initial rates of
$0.1-1\,{\mathrm{M}_{\odot}}\,\mathrm{yr}^{-1}$. Stellar evolution
calculations indicate that such flows, if they persist, build up
$10^{4}-10^{5}\,{\mathrm{M}_{\odot}}$ stars that later collapse to DCBHs
(Hirano et al. 2014; Umeda et al. 2016; Woods et al. 2017; Haemmerlé et al.
2018a – see Woods et al. 2019 and Maio et al. 2019 for recent reviews).

Cosmological simulations of the collapse of atomically-cooled halos are either
run at extremely high resolutions that can follow flows down to protostellar
radii but only for short times because of time step restrictions (e.g.,
Becerra et al., 2015, 2018; Ardaneh et al., 2018; Luo et al., 2018), or for
the longer times required to evolve the flows over many dynamical times but at
the cost of resolving fragmentation deep in the accretion disk of the star
(e.g., Chon et al., 2018; Regan & Downes, 2018; Suazo et al., 2019; Latif et
al., 2021). The first simulations to evolve atomically cooled flows for the
entire life of a SMS were Latif et al. (2020), which showed that some could
form in binaries or small clusters. Most recently, Patrick et al. (2020)
followed the collapse of a variety of atomically-cooled halos for
$2-4\,\mathrm{Myr}$ and found that multiple accretion disks could form in a
halo, each of which could harbor a SMS. These disks experienced multiple close
encounters that gravitationally torqued gas into their centers, triggering
brief bouts of massive accretion onto their respective stars. These stars thus
interacted indirectly with one another via encounters between their accretion
disks. Although the simulations suggested a variety of final outcomes for
these objects, they did not evolve the stars themselves.

Recently, Woods et al. (2021) modeled the growth of the stars of the primary
disks in Patrick et al. (2020) and found that they collapse at final masses of
$1.1\times 10^{5}\,{\mathrm{M}_{\odot}}-1.9\times
10^{5}\,{\mathrm{M}_{\odot}}$. Here, we investigate the co-evolution of the
stars in the most massive and stable of the companion disks in Patrick et al.
(2020) with the Kepler Lagrangian stellar evolution and hydrodynamics code.
Our Kepler models use the time-dependent accretion rates that were tallied at
the centers of the disks in the cosmological simulations. We follow the
evolution of the SMSs in the companion disks until the disks were subsumed
into the primary disks or torn apart by tidal interactions with them. In
Section 2 we describe our cosmological and stellar evolution models. The
evolution of the SMSs is examined in Section 3. From the interaction histories
of the disks and the formation and collapse times of the stars within them we
determine if they form SMS – SMS, SMS – DCBH or binary DCBH systems and
discuss possible consequences for their detection in Section 4. We discuss
prospects for binary interactions in these systems on much smaller scales at
later times in Section 5.

##  2 Numerical Methods

###  2.1 enzo

Accretion rates for the companion stars in our study were taken from
cosmological simulations with the enzo adaptive mesh refinement code by
Patrick et al. (2020). These simulations resolved accretion at the centers of
the disks on 0.01 pc scales with six-species primordial gas chemistry (H, He,
e-, H+, He+, and He2+) to approximate very high LW backgrounds that ensured
isothermal atomic cooling in the halos. The simulations included collisional
ionization and excitation cooling by H and He, recombination cooling by H and
He, bremsstrahlung cooling, and inverse Compton cooling by the cosmic
microwave background (Bryan et al., 2014). The accretion rates were tabulated
by computing mass fluxes through a 0.134 pc sphere at the center of each disk
at 10 kyr intervals.

The eight halos in Patrick et al. (2020) were chosen to span a range of
assembly histories and spin parameters and had masses of $\sim$ 1 - 9 $\times$
107 $\mathrm{M}_{\odot}$ at collapse at redshifts $z=$ 13.9 - 20.4. Resolution
at the smallest scales on which fragmentation is known to occur in the disks
had to be sacrificed in these simulations to follow their evolution for the
times required for the stars to form, evolve and collapse (e.g., Becerra et
al., 2015, 2018). Our accretion rates therefore exclude these fragments, but
they are expected to fall onto the protostar prior to the main sequence
(Inayoshi & Haiman, 2014). We show accretion rates only for the largest
companion disks that form in the halos in Figure 1, not smaller fragments that
could undergo partial collapse but are short lived because they are taken up
by the main disk or are tidally disrupted by it.

Figure 1: Accretion rates from Patrick et al. (2020) for the secondary stars.
Note that times here are measured from the formation of the secondary disk,
not the beginning of the simulation.

###  2.2 kepler

We follow the evolution of the stars with kepler (Weaver et al., 1978; Woosley
et al., 2002). kepler typically partitions each star into a few thousand
zones, over which it solves the energy and angular momentum equations:

$\displaystyle\frac{dv}{dt}$ $\displaystyle=$ $\displaystyle 4\pi
r^{2}\frac{\partial P}{\partial m_{r}}+\frac{4\pi}{r}\frac{\partial
Q}{\partial m_{r}}-\frac{G_{\rm{rel}}m_{r}}{r^{2}}$ (1)
$\displaystyle\frac{du}{dt}$ $\displaystyle=$ $\displaystyle-4\pi
P\frac{\partial}{\partial m_{r}}(vr^{2})+4\pi Q\frac{\partial}{\partial
m_{r}}\left(\frac{v}{r}\right)-\frac{\partial L}{\partial
m_{r}}+\epsilon~{}~{}~{}$ (2)

where $r$ is the radius, $v$ is the velocity, $P$ is the pressure, $u$ is the
internal energy per unit mass, $\epsilon$ is the local energy generation rate
per unit mass, and $L$ is the rate of energy flow through a shell of radius
$r$. We include a first-order post-Newtonian correction to the gravitational
constant, $G_{\rm{rel}}$:

$G_{\rm{rel}}=G\left(1+\frac{P}{\rho c^{2}}+\frac{4\pi
Pr^{3}}{m_{r}c^{2}}\right)\left(1-\frac{2Gm_{r}}{rc^{2}}\right)^{-1}$ (3)

In Equation 1, the factor $Q$ in the viscous term is

$Q=\frac{4}{3}\,\eta_{\nu}r^{4}\frac{\partial}{\partial
r}\left(\frac{v}{r}\right)\;,$ (4)

where $\eta_{\nu}$ is the dynamic viscosity as defined in Weaver et al.
(1978):

$\eta_{\nu}=\eta_{\mathrm{R}}+\frac{3}{4}l_{1}\rho
c_{\mathrm{s}}+\frac{3}{4}l_{2}^{2}\;\rm{max}(0,-\nabla\cdot\bf{v})\vspace{0.1in}$
(5)

which includes both the real viscosity, $\eta_{\mathrm{R}}$, and the
artificial viscosity, which can be modified arbitrarily to dampen acoustic
oscillations during quiescent phases of the evolution of a star. In this
study, we take the standard values of $l_{1}=0.1\Delta r$ and $l_{2}=2\Delta
r$, where $\Delta r$ is the width of a grid zone. At each time step, Kepler
initially attempts to make an arbitrarily large jump $\Delta t$ before
iterating to find the maximum time step permitted by preset restrictions on
the change in radius, temperature, density, luminosity, or velocity between
zones, allowing us to model long-lived, quiescent evolutionary phases and to
follow the emergence of shocks or the onset of collapse on short timescales.
The masses of the stars at collapse triggered by the the post-Newtonian
instability (Chandrasekhar, 1964) in Kepler are consistent with analytic
predictions (Haemmerlé, 2021).

To close these equations, nuclear-burning and energy generation are coupled to
the hydrodynamics and solved with an adaptive network (Woosley et al., 2004),
convection is treated in a time-dependent manner with heat transport following
the Ledoux criterion (Weaver et al., 1978), and a Helmholtz-like equation of
state is used that incorporates electron-positron pair production,
relativistic and non-relativistic degenerate and non-degenerate electrons, and
radiation (Timmes & Swesty, 2000). We neglect mass loss because wind and
pulsational mass losses are expected to be negligible relative to the
accretion rates in these stars (Vink et al., 2001; Baraffe et al., 2001;
Hosokawa et al., 2013). For simplicity, we neglect rotation although it could
affect the stellar structure (Haemmerlé et al., 2018b).

As in previous studies (Hosokawa et al., 2013; Woods et al., 2017; Haemmerlé
et al., 2018a; Woods et al., 2021), we initialize all models as 10
$\mathrm{M}_{\odot}$, $n=$ 3 polytropes with central densities
$\rho_{\mathrm{c}}=$ 10-3 g cm-3 and temperatures $T_{\mathrm{c}}=$ $1.2\times
10^{6}\,$K, i.e., at the onset of deuterium burning. We treat accretion onto
the star and the associated advection of entropy as described in Woosley et
al. (2004) and Woods et al. (2017). The star is evolved until the onset of
collapse, the end of the enzo simulation, or if the evolution time corresponds
to a moment in the enzo simulation when the companion disk interacted strongly
with or merged with the primary disk.

##  3 Evolution of Multiple Supermassive Stars in Primordial Halos

ct0cmHalo 01 |  ct0cmHalo 02
---|---
ct0cmHalo 08 |  ct0cmHalo 10
ct0cmHalo 12 |  ct0cmHalo 16
ct0cmHalo 19 |  ct0cmHalo 20
Figure 2: Kippenhahn diagrams of the internal structures of the stars in the satellite disks over time. The blue colors indicate energy generation rates from nuclear-burning, green lines denote convective regions, red lines mark semi-convective regions, and light blue lines indicate radiative/convectively neutral regions. Halo | $M_{1}$ | $t_{sim}$ | Evol Stage (1) | $X_{c,1}$ | $M_{2}$ | $t_{2}$ | Evol Stage (2) | $X_{c,2}$ | Close Pair?
---|---|---|---|---|---|---|---|---|---
| k$\mathrm{M}_{\odot}$ | Myr |  |  | k$\mathrm{M}_{\odot}$ | Myr |  |  |
01 | 134 | 0.947 | MS | 0.38 | 165 | 0.571 | MS | 0.40 | yes
02 | 143 | 1.46 | BH | 0.34 | 62 | 0.378 | MS | 0.58 | yes
08 |  $\geq$186 | 1.014 | BH | — | 99 | 0.503 | MS | 0.56 | yes
10 |  $\geq$132 | 2.006 | BH | — | 53 | 0.639 | MS | 0.55 | no
12 |  $\gtrsim$178 | 1.395 | BH | — | 79 | 0.380 | MS | 0.57 | yes
16 |  $\geq$109 | 3.054 | BH | — | 186 | 1.27 | MS (near collapse?) | 0.39 | no
19 | 46 | 1.439 | MS | 0.28 | 133 | 0.426 | MS | 0.54 | yes
20 |  $\geq$178 | 1.770 | BH | — | 154 | 1.103 | MS (near collapse?) | 0.34 | no
Table 1: Final masses, ages, evolutionary stages (ES), and central hydrogen
fractions for the primary (1) and secondary (2) SMSs in each halo, followed by
whether or not they form a close pair by the end of the simulation. Note that
the SMS in the primary disk in Halo 02 collapses to a black hole during the
merger with the secondary disk.

We find that all the secondary disks in the halos in Patrick et al. (2020)
form long-lived nuclear-burning stars. Like the stars in the primary disks
(Woods et al., 2021), we find a striking variety of internal structures for
these secondary SMSs as shown in Figure 2, in marked contrast to stars that
evolve under constant accretion rates (e.g., Hosokawa et al., 2013; Woods et
al., 2017; Haemmerlé et al., 2018a). The stars all initially exhibit a deep
radiative envelope corresponding to the surge in accretion associated with the
formation of their natal disks, which proceeds on much shorter times than the
star’s thermal timescale and leads to the buildup of a steep entropy gradient
(e.g., Begelman, 2010; Hosokawa et al., 2013; Woods et al., 2017; Haemmerlé et
al., 2018a). In several cases, such as Halos 01 and 16, however, almost all of
mass of the star eventually lies within its convective core because there is a
long ($\gtrsim$ 100 kyr) quiescent phase ($\gg$ Kelvin-Helmholtz time-scale)
in accretion in which the star can thermally relax and its structure can
approach that of an $n=3$ polytrope (Chandrasekhar, 1964; Woods et al., 2020).
In other models, we see the formation of both transient and long-lived
convective cells in the otherwise deep, high-entropy radiative envelopes of
some rapidly-accreting stars, similar to those in some constant-accretion rate
models (e.g., Woods et al., 2017; Haemmerlé et al., 2018a).

In Halo 01, a secondary disk forms at 400 kyr by breaking off from one of the
spiral arms of the primary disk. Throughout its evolution, this companion disk
remains in an elliptical orbit around the first with typical separations of
0.5 - 1 pc. The sharp spike in accretion $\sim$300 kyr after formation is due
to a merger with a clump in the halo. Accretion continues until the second
disk merges with the first 650 kyr after formation. At this time, the SMS in
the second disk is a somewhat evolved main sequence star with a central
hydrogen fraction of $\sim$0.4. The SMS in the first disk is only slightly
more evolved, with a central hydrogen fraction of 0.38. Although the
subsequent evolution of these stars is beyond the scope of this work, Woods et
al. (2021) assumed that the collision of the disks would trigger a surge in
accretion onto the first star but the presence of the second SMS may lead to
other outcomes, as discussed in $\S$ 4.1.

Halo 02 is an example of a particularly turbulent, chaotic system in which
three clumps form and interact with the primary and companion disks (Patrick
et al., 2020). Interactions between these clumps and the second disk lead to
two large bursts of accretion at 130 kyr and 250 kyr that create a SMS that
grows to $\mathord{\sim}60$ k$\mathrm{M}_{\odot}$ by 275 kyr and has a deep
radiative envelope. At this point, there is a brief quiescent phase in
accretion due to a close encounter between the two disks.The two disks then
merge 378 kyr after the formation of the second disk. The first and second
SMSs are 62 k$\mathrm{M}_{\odot}$ and 143 k$\mathrm{M}_{\odot}$, respectively,
and have central hydrogen fractions of 0.58 and 0.34. Notably, the collision
of the disks produces a large spike in accretion onto the first SMS that
quickly brings it up to the post-Newtonian instability and causes it to
collapse (Woods et al., 2021).

Only one satellite disk forms in Halo 08, which lasts from 685 kyr to about 1
Myr after the formation of the first disk. Accretion onto the star is highly
variable because of the eccentric orbit of the satellite disk, and the
fluctuations correlate with the smallest and largest distances between the two
disks, which are $\sim$ 0.3 pc and 2.0 pc, respectively. This accretion
history produces the distinctive step-like structure in the Kippenhahn diagram
of the star in Figure 2. About 1 Myr after the formation of the primary disk,
the star in the second disk has an age of 503 kyr, a mass of 99
k$\mathrm{M}_{\odot}$, and a central hydrogen fraction of 0.56. The SMS in the
primary disk, however, has by this time collapsed to a black hole, having
encountered the post-Newtonian instability at an age of 950 kyr when it
reached a mass of 186 k$\mathrm{M}_{\odot}$ and a central hydrogen fraction of
0.38.

The first disk in Halo 10 is stable and does not begin to fragment for 1.14
Myr. The second disk forms $\sim$ 1.2 pc from the center of the first in an
initially highly elliptical orbit. Two major episodes of accretion drive the
growth of the second SMS that in each case build up a deep, high-entropy
envelope. Between these episodes there is a long (nearly 100 kyr) quiescent
phase during which the star thermally relaxes without becoming entirely
convective. This star reaches a mass of 53 k$\mathrm{M}_{\odot}$ and a central
hydrogen fraction of 0.55 by the end of the simulation. At this time the
secondary disk is still in an elliptical orbit around the first with
separations that vary from 0.5 – 1 pc. The SMS in the first disk, however, has
collapsed because of the post-Newtonian Chandrasekhar instability late on the
main sequence at an age of 1.95 Myr, a final mass of 132
k$\mathrm{M}_{\odot}$, and a central hydrogen fraction of 0.06.

The primary disk in Halo 12 is highly turbulent and frequently fragments, but
most of the clumps soon migrate to the center of the disk. The longest-lived
of the clumps forms at $\sim$ 1 Myr and produces its own disk that survives
for 380 kyr. Rapid accretion in this disk creates a SMS with a deep radiative
envelope that grows to a mass 79 k$\mathrm{M}_{\odot}$ with a central hydrogen
fraction of 0.57 by the time the two disks merge. About 200 kyr before the
destruction of the second disk, the SMS in the primary disk collapsed via the
Chandrasekhar instability while still on the main sequence at a final mass of
178 k$\mathrm{M}_{\odot}$ and central hydrogen fraction of 0.32.

Like Halo 10, Halo 16 has a relatively stable disk that does not fragment for
$\sim$ 1.5 Myr. A key difference, however, is that the SMS in the primary disk
collapses to a BH at about the time the second disk forms because it
encounters the post-Newtonian instability at the very end of the main sequence
at a mass of 109 k$\mathrm{M}_{\odot}$, when its core hydrogen is exhausted.
It is unclear from the Enzo simulation how X-rays from this black hole would
affect the evolution of the companion disk but we discuss some possible
outcomes in $\S$4.2. As the primary disk continues to fragment, 3-body
interactions fling the companion disk into a highly elliptical orbit $\sim$
300 kyr after its formation. The initial burst of accretion is followed by a
long quiescent phase ($\sim$ 0.5 Myr) during which the second SMS, which is
now $\sim$ 40 k$\mathrm{M}_{\odot}$, becomes almost completely thermally
relaxed and almost fully convective (Woods et al., 2020). The rapid accretion
beginning at $\sim$ 750 kyr builds up a massive radiative envelope on top of
this convective core and the mass of the star grows to 186
k$\mathrm{M}_{\odot}$ by the end of the simulation. At this point the SMS is
still on the main sequence in a disk on a long orbit around the first disk,
but it appears to be on the verge of collapse. With a mass nearing the upper
limit for SMSs with similar accretion rates (e.g., Woods et al., 2017;
Haemmerlé et al., 2018a; Woods et al., 2021), the star is unlikely to survive
for much longer so the system may soon produce a DCBH binary.

As the primary disk in Halo 19 begins to fragment, two clumps merge and form a
stable companion disk after $\sim$ 800 kyr. The initial burst of accretion due
to the formation of this disk is particularly strong and builds up a star with
a massive convective envelope in $\sim$ 500 kyr. Accretion in the disk is so
rapid that the evolution of the second SMS overtakes that of the first. By the
time the disks merge, the first SMS has reached a mass of 46
k$\mathrm{M}_{\odot}$ and a central hydrogen fraction of 0.28 but the second
is at a mass of 133 k$\mathrm{M}_{\odot}$ and a central hydrogen fraction of
0.54. Both stars are still on the main sequence.

The primary disk in Halo 20 begins to fragment soon after formation and
quickly forms a single, massive companion disk. Accretion rates in the disk
are fairly high, and are rejuvenated at one point when a clump collides with
the disk. These high rates quickly build up a SMS with a deep, high-entropy
envelope. The companion disk later exchanges mass with the primary disk over a
number of orbits that mostly halts accretion onto the second star and allows
it to thermally relax. In the meantime, the SMS in the primary disk collapses
at 1.48 Myr to a BH via the post-Newtonian/Chandrasekhar instability late in
the main sequence at a core hydrogen fraction of 0.14 and a mass of 189
k$\mathrm{M}_{\odot}$. At 1.77 Myr, the end of the simulation, the companion
disk is in a relatively long elliptical orbit ($\sim$ 1.8 pc separation) that
is growing in radius. At this point, the companion is still on the main
sequence but, as in Halo 16, appears to be about to collapse, having reached a
mass of 154 k$\mathrm{M}_{\odot}$.

##  4 Discussion

In most of the cases above, the accretion history and evolution of the second
SMS ends with the merger of its host disk with the primary disk. Although the
subsequent evolution of the stars or DCBHs in the disks cannot be determined
here, a number of outcomes can be inferred for them for follow-up with future
simulations, as we outline below.

###  4.1 MS-MS Pairs

The primary and secondary disks in Halo 01, Halo 19, and nominally Halo 02
merge while their central stars are still on the main sequence (MS). With a
maximum physical resolution of 0.014 pc ($\sim$ 3 kAU), the Enzo simulations
cannot determine if they later interact. Drag forces in the inner disk could
bring about a swift merger between the two SMSs or the stars could carve out a
gap in the circumbinary disk. Just before the merger of the two disks in Halo
02 one of the SMSs appears to be on the verge of collapse via the post-
Newtonian instability. Examples such as this make it clear that a more careful
treatment of stellar structure during the resulting enormous surge in
accretion is essential (Menon & Heger, 2017) and could be informed by
smoothed-particle hydrodynamics calculations (e.g., Glebbeek et al., 2013).
Such mergers may also produce unique observational transients that signal the
formation of SMSs in the early Universe. If instead these supermassive stellar
pairs evolve into relatively long-lived binaries, they could later become BH –
MS or BH – BH pairs.

###  4.2 MS-BH Pairs

In Halos 08 and 12 the primary and companion disks merge after the star in the
primary disk has collapsed to a DCBH while the SMS in the other disk is still
on the main sequence. Given that the masses of both objects are comparable in
both halos, later interactions between the two are unlikely to produce tidal
disruption events. The tidal disruption radius

$R_{\rm t}=R_{*}\left(\frac{M_{\rm BH}}{M_{*}}\right)^{\frac{1}{3}}$ (6)

in such cases is of the order of the stellar radius and would therefore only
happen during a collision. However, fragmentation on scales smaller than those
resolved here may lead to the formation of much less massive stars that could
then be torn apart by the DCBH and produce a powerful radio or NIR transient
(Kashiyama & Inayoshi, 2016; Regős et al., 2021).

If these objects form a long-lived binary, they will continue to evolve
subject to additional torques from the circumbinary disk. Depending on the
orbital evolution, the companion SMS may overflow its Roche lobe while still
on the main sequence or as it evolves, in which case its deep radiative
envelope could sustain long-lived mass exchange and, perhaps, the formation of
a “supermassive” X-ray binary (e.g., Soberman et al., 1997). Eventually, the
collapse of the second SMS due to exhaustion of core nuclear fuel or the
Chandrasekhar instability will produce a DCBH binary as discussed below.

###  4.3 BH-BH Pairs

Some of our models may produce long-period DCBH binaries. At the end of the
simulation in Halo 10 the second SMS was still on the main sequence but its
host disk was in a highly elliptical orbit with the first disk, whose star had
collapsed to a DCBH. This long-period orbit suggests that the SMS will
collapse to a DCBH before interacting closely with the other BH. Although the
secondary SMSs in Halos 16 and 20 are also still on the main sequence at the
end of the simulation they are both on the verge of collapse via the post-
Newtonian instability unless accretion stops. Both of these halos could thus
produce long-period ($\sim 1$ pc separation) binary DCBHs. The timescale for a
merger between the two BHs would depend on their dynamics in the halo and drag
forces due to gas flows therein, and the strength of their gravitational wave
(GW) signal depends on their masses when they finally do merge. Here, we
estimate the GW emission from the merger of the DCBH binary in Halo 20 (Robson
et al., 2019)111LISA Sensitivity Calculator available for download from:
github.com/eXtremeGravityInstitute/LISA_Sensitivity.git, considering a sky and
polarization averaged transfer function and assuming that the orbit has
circularized by prior GW emission. The expected signal is shown in Figure 3
for the redshift at which the halo collapses, $z=17.7$. This merger will be
detectable by LISA, as would mergers of the other BH – BH pairs in our halos
at their respective redshifts.

##  5 Conclusions

A growing number of recent numerical simulations indicate that the formation
of multiple very massive stars (e.g., Wise et al., 2019; Regan et al., 2020)
or SMSs (e.g., Latif et al., 2020; Patrick et al., 2020) was common in
primordial, atomic-cooling halos, and that these could be the origin of the
first quasars in the Universe. Here, we have shown that secondary disks in
these halos form SMSs whose accretion histories can be coupled to those of the
stars in the primary disks via tidal interactions. Most satellite disks merge
with the primary disk in just 1 – 2 Myr, which could yield supermassive
stellar mergers, supermassive X-ray binaries, massive black hole seed mergers,
and SMS binaries. Such interactions could profoundly affect the evolution of
the stars, the SMBH seeds they produce, and their observational signatures.
Simulations that follow the dynamics of these binaries and mass exchange
between them with stellar evolution calculations will be the aim of future
studies.

Figure 3: Sky-averaged LISA sensitivity curve (orange) and the observable
characteristic strain (blue) from the merger of the DCBHs in Halo 20.

T.E.W. acknowledges support from the National Research Council Canada’s
Plaskett Fellowship. S.P. was supported by STFC grant ST/N504245/1 and D.J.W.
was supported by the Ida Pfeiffer Professorship at the Institute of
Astrophysics at the University of Vienna. A.H. was supported by the Australian
Research Council (ARC) Centre of Excellence (CoE) for Gravitational Wave
Discovery (OzGrav) through project number CE170100004, by the ARC CoE for All
Sky Astrophysics in 3 Dimensions (ASTRO 3D) through project number
CE170100013, and by the National Science Foundation under Grant No.
PHY-1430152 (JINA Center for the Evolution of the Elements, JINA-CEE).

## References

  * Agarwal et al. (2012) Agarwal, B., Khochfar, S., Johnson, J. L., et al. 2012, MNRAS, 425, 2854, doi: 10.1111/j.1365-2966.2012.21651.x
  * Ardaneh et al. (2018) Ardaneh, K., Luo, Y., Shlosman, I., et al. 2018, MNRAS, 479, 2277, doi: 10.1093/mnras/sty1657
  * Bañados et al. (2018) Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473, doi: 10.1038/nature25180
  * Baraffe et al. (2001) Baraffe, I., Heger, A., & Woosley, S. E. 2001, ApJ, 550, 890, doi: 10.1086/319808
  * Becerra et al. (2015) Becerra, F., Greif, T. H., Springel, V., & Hernquist, L. E. 2015, MNRAS, 446, 2380, doi: 10.1093/mnras/stu2284
  * Becerra et al. (2018) Becerra, F., Marinacci, F., Bromm, V., & Hernquist, L. E. 2018, MNRAS, 480, 5029, doi: 10.1093/mnras/sty2210
  * Begelman (2010) Begelman, M. C. 2010, MNRAS, 402, 673, doi: 10.1111/j.1365-2966.2009.15916.x
  * Bosman (2021) Bosman, S. E. I. 2021, All z¿5.7 quasars currently known, 1.11, Zenodo, doi: 10.5281/zenodo.4610490
  * Bromm & Loeb (2003) Bromm, V., & Loeb, A. 2003, ApJ, 596, 34, doi: 10.1086/377529
  * Bryan et al. (2014) Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19, doi: 10.1088/0067-0049/211/2/19
  * Chandrasekhar (1964) Chandrasekhar, S. 1964, ApJ, 140, 417, doi: 10.1086/147938
  * Chon et al. (2018) Chon, S., Hosokawa, T., & Yoshida, N. 2018, MNRAS, 475, 4104, doi: 10.1093/mnras/sty086
  * Dijkstra et al. (2008) Dijkstra, M., Haiman, Z., Mesinger, A., & Wyithe, J. S. B. 2008, MNRAS, 391, 1961, doi: 10.1111/j.1365-2966.2008.14031.x
  * Fernandez et al. (2014) Fernandez, R., Bryan, G. L., Haiman, Z., & Li, M. 2014, MNRAS, 439, 3798, doi: 10.1093/mnras/stu230
  * Glebbeek et al. (2013) Glebbeek, E., Gaburov, E., Portegies Zwart, S., & Pols, O. R. 2013, MNRAS, 434, 3497, doi: 10.1093/mnras/stt1268
  * Haemmerlé (2021) Haemmerlé, L. 2021, A&A, 647, A83, doi: 10.1051/0004-6361/202039686
  * Haemmerlé et al. (2018a) Haemmerlé, L., Woods, T. E., Klessen, R. S., Heger, A., & Whalen, D. J. 2018a, MNRAS, 474, 2757, doi: 10.1093/mnras/stx2919
  * Haemmerlé et al. (2018b) —. 2018b, ApJ, 853, L3, doi: 10.3847/2041-8213/aaa462
  * Hirano et al. (2017) Hirano, S., Hosokawa, T., Yoshida, N., & Kuiper, R. 2017, Science, 357, 1375, doi: 10.1126/science.aai9119
  * Hirano et al. (2014) Hirano, S., Hosokawa, T., Yoshida, N., et al. 2014, ApJ, 781, 60, doi: 10.1088/0004-637X/781/2/60
  * Hosokawa et al. (2013) Hosokawa, T., Yorke, H. W., Inayoshi, K., Omukai, K., & Yoshida, N. 2013, ApJ, 778, 178, doi: 10.1088/0004-637X/778/2/178
  * Inayoshi & Haiman (2014) Inayoshi, K., & Haiman, Z. 2014, MNRAS, 445, 1549, doi: 10.1093/mnras/stu1870
  * Kashiyama & Inayoshi (2016) Kashiyama, K., & Inayoshi, K. 2016, ApJ, 826, 80, doi: 10.3847/0004-637X/826/1/80
  * Latif et al. (2014) Latif, M. A., Bovino, S., Van Borm, C., et al. 2014, MNRAS, 443, 1979, doi: 10.1093/mnras/stu1230
  * Latif et al. (2021) Latif, M. A., Khochfar, S., Schleicher, D., & Whalen, D. J. 2021, MNRAS, doi: 10.1093/mnras/stab2708
  * Latif et al. (2020) Latif, M. A., Khochfar, S., & Whalen, D. 2020, ApJ, 892, L4, doi: 10.3847/2041-8213/ab7c61
  * Latif et al. (2013) Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013, MNRAS, 430, 588, doi: 10.1093/mnras/sts659
  * Lodato & Natarajan (2006) Lodato, G., & Natarajan, P. 2006, MNRAS, 371, 1813, doi: 10.1111/j.1365-2966.2006.10801.x
  * Luo et al. (2018) Luo, Y., Ardaneh, K., Shlosman, I., et al. 2018, MNRAS, 476, 3523, doi: 10.1093/mnras/sty362
  * Maio et al. (2019) Maio, U., Borgani, S., Ciardi, B., & Petkova, M. 2019, PASA, 36, e020, doi: 10.1017/pasa.2019.10
  * Menon & Heger (2017) Menon, A., & Heger, A. 2017, MNRAS, 469, 4649.  https://arxiv.org/abs/1703.04918
  * Mortlock et al. (2011) Mortlock, D. J., Warren, S. J., Venemans, B. P., et al. 2011, Nature, 474, 616, doi: 10.1038/nature10159
  * Patrick et al. (2020) Patrick, S. J., Whalen, D. J., Elford, J. S., & Latif, M. A. 2020, arXiv e-prints, arXiv:2012.11612.  https://arxiv.org/abs/2012.11612
  * Regan & Downes (2018) Regan, J. A., & Downes, T. P. 2018, MNRAS, 478, 5037, doi: 10.1093/mnras/sty1289
  * Regan & Haehnelt (2009) Regan, J. A., & Haehnelt, M. G. 2009, MNRAS, 396, 343, doi: 10.1111/j.1365-2966.2009.14579.x
  * Regan et al. (2020) Regan, J. A., Wise, J. H., Woods, T. E., et al. 2020, The Open Journal of Astrophysics, 3, 15, doi: 10.21105/astro.2008.08090
  * Regős et al. (2021) Regős, E., Vinkó, J., & Stermeczky, Z. V. 2021, ApJ, 909, 64, doi: 10.3847/1538-4357/abdf51
  * Robson et al. (2019) Robson, T., Cornish, N. J., & Liu, C. 2019, Classical and Quantum Gravity, 36, 105011, doi: 10.1088/1361-6382/ab1101
  * Schauer et al. (2017) Schauer, A. T. P., Agarwal, B., Glover, S. C. O., et al. 2017, MNRAS, 467, 2288, doi: 10.1093/mnras/stx264
  * Soberman et al. (1997) Soberman, G. E., Phinney, E. S., & van den Heuvel, E. P. J. 1997, A&A, 327, 620.  https://arxiv.org/abs/astro-ph/9703016
  * Suazo et al. (2019) Suazo, M., Prieto, J., Escala, A., & Schleicher, D. R. G. 2019, ApJ, 885, 127, doi: 10.3847/1538-4357/ab45eb
  * Tanaka & Li (2014) Tanaka, T. L., & Li, M. 2014, MNRAS, 439, 1092, doi: 10.1093/mnras/stu042
  * Timmes & Swesty (2000) Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501, doi: 10.1086/313304
  * Umeda et al. (2016) Umeda, H., Hosokawa, T., Omukai, K., & Yoshida, N. 2016, ApJ, 830, L34, doi: 10.3847/2041-8205/830/2/L34
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574, doi: 10.1051/0004-6361:20010127
  * Wang et al. (2021) Wang, F., Yang, J., Fan, X., et al. 2021, ApJ, 907, L1, doi: 10.3847/2041-8213/abd8c6
  * Weaver et al. (1978) Weaver, T. A., Zimmerman, G. B., & Woosley, S. E. 1978, ApJ, 225, 1021, doi: 10.1086/156569
  * Wise et al. (2019) Wise, J. H., Regan, J. A., O’Shea, B. W., et al. 2019, Nature, 566, 85, doi: 10.1038/s41586-019-0873-4
  * Wise et al. (2008) Wise, J. H., Turk, M. J., & Abel, T. 2008, ApJ, 682, 745, doi: 10.1086/588209
  * Woods et al. (2020) Woods, T. E., Heger, A., & Haemmerlé, L. 2020, MNRAS, 494, 2236, doi: 10.1093/mnras/staa763
  * Woods et al. (2017) Woods, T. E., Heger, A., Whalen, D. J., Haemmerlé, L., & Klessen, R. S. 2017, ApJ, 842, L6, doi: 10.3847/2041-8213/aa7412
  * Woods et al. (2021) Woods, T. E., Patrick, S., Elford, J. S., Whalen, D. J., & Heger, A. 2021, ApJ, 915, 110, doi: 10.3847/1538-4357/abfaf9
  * Woods et al. (2019) Woods, T. E., Agarwal, B., Bromm, V., et al. 2019, Publications of the Astronomical Society of Australia, 36, e027, doi: 10.1017/pasa.2019.14
  * Woosley et al. (2002) Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Reviews of Modern Physics, 74, 1015, doi: 10.1103/RevModPhys.74.1015
  * Woosley et al. (2004) Woosley, S. E., Heger, A., Cumming, A., et al. 2004, ApJS, 151, 75, doi: 10.1086/381533
  * Wu et al. (2015) Wu, X.-B., Wang, F., Fan, X., et al. 2015, Nature, 518, 512, doi: 10.1038/nature14241
  * Yoshida et al. (2003) Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003, ApJ, 592, 645, doi: 10.1086/375810