# Single-Photon Storage in a Ground-State Vapor Cell Quantum Memory

Gianni Buser<EMAIL_ADDRESS>Departement Physik, Universität Basel,
Klingelbergstr. 82, 4056 Basel, Switzerland.  Roberto Mottola  Departement
Physik, Universität Basel, Klingelbergstr. 82, 4056 Basel, Switzerland.  Björn
Cotting  Departement Physik, Universität Basel, Klingelbergstr. 82, 4056
Basel, Switzerland.  Janik Wolters  Departement Physik, Universität Basel,
Klingelbergstr. 82, 4056 Basel, Switzerland.  Deutsches Zentrum für Luft- und
Raumfahrt e.V. (DLR), Institute of Optical Sensor Systems, Rutherfordstr. 2,
12489 Berlin, Germany.  Technische Universität Berlin, Institut für Optik und
Atomare Physik, Hardenbergstr. 36, 10623 Berlin, Germany.  Philipp Treutlein
Departement Physik, Universität Basel, Klingelbergstr. 82, 4056 Basel,
Switzerland.

###### Abstract

Interfaced single-photon sources and quantum memories for photons together
form a foundational component of quantum technology. Achieving compatibility
between heterogeneous, state-of-the-art devices is a long-standing challenge.
We built and successfully interfaced a heralded single-photon source based on
cavity-enhanced spontaneous parametric down-conversion in ppKTP and a matched
memory based on electromagnetically induced transparency in warm ^87Rb vapor.
The bandwidth of the photons emitted by the source is
$370\text{\,}\mathrm{MHz}$, placing its speed in the technologically relevant
regime while remaining well within the acceptance bandwidth of the memory.
Simultaneously, the experimental complexity is kept low, with all components
operating at or above room temperature. Read-out noise of the memory is
considerably reduced by exploiting polarization selection rules in the
hyperfine structure of spin-polarized atoms. For the first time, we
demonstrate single-photon storage and retrieval in a ground-state vapor cell
memory, with $g_{c,\,\text{ret}}^{(2)}=$0.177(23)$$ demonstrating the single-
photon character of the retrieved light. Our platform of single-photon source
and atomic memory is attractive for future experiments on room-temperature
quantum networks operating at high bandwidth.

##  I Introduction

Quantum memories combined with single photons from high quality sources are
versatile and indispensable building blocks across the fields of quantum
communication and information. They are at the heart of each node and
interconnect in visions of a quantum internet [1, 2], and central to the
standard paradigm of quantum repeaters [3, 4, 5]. They can be used to
synchronize probabilistic gate operations and sources [6, 7], and can even
improve the indistinguishability of photons emitted by quantum dots through
filtering [8]. Further prospective applications include linear optical quantum
computing, metrology, and photon detection [9, 10]. These applications put
different requirements on a memory, calling for a quantitative assessment of
memory performance with numerous figures of merit [11], including fidelity,
efficiency, storage time, bandwidth, and various compatibility parameters.

Memories implemented in the ground state of room-temperature atomic vapors
perform well in terms of fidelity, efficiency, storage times, and bandwidth
[12, 13, 14, 15]. Moreover, they are compatible with high quality single-
photon sources based on spontaneous parametric downconversion [15, 16] or
semiconductor quantum dots [17, 18, 19, 20]. Together with their technological
simplicity, this renders atomic vapor cells a promising memory system for
quantum networks, potentially even ones deployed in space [21, 22]. However, a
long-standing problem of warm atomic vapor memories is read-out noise arising
from four-wave mixing [23, 24] and collisional fluorescence [25, 26], which
degrades the quality of the retrieved photons. Consequently, such memories are
commonly tested with laser pulses attenuated to the single-photon level,
circumventing the stricter requirements on memory noise imposed by real single
photon sources with imperfect efficiencies. Readout noise can be suppressed in
cold atom systems [27, 28, 29, 30], or with excited-state storage schemes [31,
32], but these approaches come at the price of much higher experimental
complexity, or fundamentally limited storage times, respectively.
Demonstrating storage and retrieval of single photons in ground-state vapor
cell memories, as characterized by non-classical photon number statistics of
the retrieved light, has so far remained elusive [15, 33, 12, 34].

Here we report the storage and retrieval of single photons in a ground-state
atomic vapor cell quantum memory. Our memory scheme suppresses readout noise
by exploiting polarization selection rules in the atomic hyperfine structure
and by operating at a bandwidth much higher than the excited state’s radiative
decay rate. We interface the atomic memory with a single photon source based
on cavity-enhanced spontaneous parametric downconversion (SPDC), which we
built for this purpose with improved operation and performance characteristics
compared to our earlier work [16]. Single photons from this source are stored
in the atomic memory and retrieved with decidedly non-classical photon number
statistics, opening up many further possibilities for quantum networking
experiments at high bandwidth in a room-temperature system.

Figure 1:  (a) Basic scheme of the experiment. (b) Energy levels of the photon
source. By spontaneous parametric downconversion in ppKTP,
$404\text{\,}\mathrm{n}\mathrm{m}$ pump photons are converted to
$795\text{\,}\mathrm{n}\mathrm{m}$ signal and
$822\text{\,}\mathrm{n}\mathrm{m}$ idler photons. (c) Energy level scheme for
the atomic memory. Optical pumping and selection rules are used to isolate a
four-level lambda system on the ^87Rb D1 line. The $\sigma^{-}$ polarized
signal is detuned by $\Delta$ from the $F=2\rightarrow F^{\prime}=1$
transition and the $\sigma^{+}$ polarized control is equally detuned from
$F=1\rightarrow F^{\prime}=1$. Photons can thus be stored as a spin wave
excitation between the initially prepared ground state $\ket{g}$ and the
storage state $\ket{s}$ via the excited states $\ket{e_{1}}$ and
$\ket{e_{2}}$.

##  II Memory Scheme

An overview of the single-photon source and quantum memory setup is shown in
Fig. 1. The memory operates on the ^87Rb D1 line at $795\text{\,}\mathrm{nm}$
in a hot atomic vapor. We initially prepare the atoms in the stretched Zeeman
ground state $\ket{g}=\ket{F=2,m_{F}=2}$ by optical pumping. This allows us to
exploit polarization selection rules to isolate a four-level lambda system
formed by the two ground states $\ket{g}$ and $\ket{s}=\ket{F=1,m_{F}=0}$ and
the excited states $\ket{e_{1}}=\ket{F^{\prime}=1,m_{F}^{\prime}=1}$ and
$\ket{e_{2}}=\ket{F^{\prime}=2,m_{F}^{\prime}=1}$ as shown in Fig. 1 (c). In
the storage process, a circularly polarized ($\sigma^{-}$) signal, that is the
single photon to be stored, is reversibly mapped to an atomic ground-state
superposition between $\ket{g}$ and $\ket{s}$ by the ($\sigma^{+}$) control
laser. Applying the control laser pulse again after the storage time recreates
the photon in the signal mode [35].

Our scheme overcomes two significant limitations of lambda-scheme atomic
memories that do not control the Zeeman state, namely their susceptibility to
four-wave mixing noise [23, 24, 15] and the presence of parasitic single-
photon transitions [33]. The former arises from off-resonant coupling of the
strong control laser to the initially prepared atomic state. The latter occur
when the signal is absorbed on a transition that would require a selection-
rule forbidden mapping to the storage state by the control. In hyperfine
lambda-schemes using $\pi$-polarized light this can occur with atoms initially
in the $m_{F}=0$ state [36]. Both of these problems are addressed
simultaneously by controlling the Zeeman state of the atoms and exploiting
polarization selection rules [37].

Since the two storage pathways involving $\ket{e_{1}}$ and $\ket{e_{2}}$
interfere destructively [38], this approach only works if the detuning of
signal and control light is lower or comparable to the excited state
splitting, and not in between the states, so that one of the transitions
dominates. This interference leads to an effective reduction in the optical
depth of the ensemble, but not to absorption without storage. Therefore
detunings within a $\mathrm{GHz}$ range red from $F^{\prime}=1$ or blue from
$F^{\prime}=2$ can be considered and optimized for signal-to-noise ratio (SNR)
and efficiency. The final working point used,
$\Delta=-2\pi\times$700\text{\,}\mathrm{MHz}$$, is the result of empirical
optimization, and the exact value is less crucial than minimizing the two-
photon detuning between the signal and control.

We operate our memory in the technologically relevant regime of large
bandwidths, typically several hundred MHz, much larger than the excited-state
radiative decay rate of the ^87Rb D1 line of
$2\pi\times$5.75\text{\,}\mathrm{MHz}$$. This allows us to also significantly
suppress noise due to collisional fluorescence [25, 26] by time-gating the
signal. Overall, our memory scheme thus eliminates several main limitations of
previous attempts to store single photons in the ground-state of atomic
vapors.

Figure 2:  Experimental setup: ECDL external cavity diode laser, PC Pockels
cell, ppKTP monolithic periodically poled potassium titanyl phosphate cavity,
LP optical longpass, PBS polarizing beamsplitter, IF interference filter, SPAD
single photon avalanche diode, DL $60\text{\,}\mathrm{m}$ fiber delay line,
DDG digital delay generator, AWG arbitrary waveform generator, SOA
semiconductor optical amplifier (fiber connections omitted), DM dichroic
mirror, EOM electro-optic modulator (fiber connections omitted), TA tapered
amplifier, CP calcite polarizing prism, $\lambda/2$ half-wave plate,
$\lambda/4$ quarter-wave plate, $\upmu$MS 4-layer mu-metal magnetic shield,
HBT Hanbury Brown and Twiss configured single photon detectors, M fiber
connection for monitoring the control. The labels S, P, and C represent the
fiber connections of the signal, pump, and control to the memory respectively.

##  III Memory Setup

The implementation of the memory is further detailed in Fig. 2. The atomic
vapor cell at the heart of the memory is a commercial $75\text{\,}\mathrm{mm}$
long quartz cylinder with $19\text{\,}\mathrm{mm}$ outer diameter and wedged
windows, which contains enriched ^87Rb ($<$1\text{\,}\mathrm{\char 37\relax}$$
^85Rb specified), $5\text{\,}\mathrm{Torr}$ of N2 buffer gas, and paraffin
coating on the walls. It is housed inside a 4-layer magnetic shield (Twinleaf
MS-1L). A simple heater maintains an atomic temperature of
$50(1)\text{\,}\mathrm{\SIUnitSymbolCelsius}$, yielding an optical depth on
the signal transition of about 25. The atomic state preparation is performed
with $1\text{\,}\mathrm{mm}$ $e^{-2}$ diameter, circularly polarized laser
beams, pumping on the ^87Rb D1 hyperfine transition
$\ket{F=2,m_{F}}\rightarrow\ket{F^{\prime}=2,m_{F}^{\prime}=m_{F}+1}$ and
repumping on the ^87Rb D2 hyperfine transition
$\ket{F=1,m_{F}}\rightarrow\ket{F^{\prime},m_{F}^{\prime}=m_{F}+1}$, with
around $20\text{\,}\mathrm{mW}$ and $10\text{\,}\mathrm{mW}$ CW power,
respectively. The effectiveness of the state preparation was characterized at
an atomic temperature of $70\text{\,}\mathrm{\SIUnitSymbolCelsius}$ with pump-
probe measurements to exclude concerns of radiation trapping at high density
[39], and is estimated to be $>$98\text{\,}\mathrm{\char 37\relax}$$
initially, decaying exponentially with
$\tau=$5(1)\text{\,}\mathrm{\SIUnitSymbolMicro s}$$.

The laser setup for pumping and control is schematically shown in Fig. 2. The
pumping beams are generated by CW external cavity diode lasers seeding fiber-
integrated semiconductor optical amplifiers (SOAs) for fast switching via the
SOA current in response to a herald trigger. They are combined on a dichroic
mirror, then fiber coupled and overlapped with the control beam under a small
angle of $2.95(15)\text{\,}\mathrm{mrad}$. The Gaussian control pulses are
generated on demand with a fiber-integrated electro-optical amplitude
modulator (Jenoptik AM785), controlled by a fast arbitrary pulse generator
(PicoQuant PPG512). The FWHM of these pulses measured before the vapor cell is
$3.77(4)\text{\,}\mathrm{ns}$. These pulses are amplified with a tapered
amplifier (TA), then spectrally filtered to remove the background of amplified
spontaneous emission of the TA with two narrowband interference filters (IFs,
$0.5\text{\,}\mathrm{nm}$ FWHM specified by manufacturer at
$795\text{\,}\mathrm{nm}$) and two monolithic etalons
($550(10)\text{\,}\mathrm{MHz}$ bandwidth, $25.5(4)\text{\,}\mathrm{GHz}$ free
spectral range), and finally fiber coupled to bring them to the memory with a
maximum possible peak power on the atoms around $680(40)\text{\,}\mathrm{mW}$.
In the experiments detailed below the power is adjusted to yield a peak Rabi
frequency of $\Omega=2\pi\times$400(30)\text{\,}\mathrm{MHz}$$ on the
$\ket{s}\rightarrow\ket{e_{1}}$ transition.

Control and signal are initially in orthogonal linear polarizations, and are
thus combined on a single calcite prism. Then a quarter-wave plate prepares
the required circular polarizations. The signal (control) is focused to a
$e^{-2}$ diameter of $480(6)\text{\,}\mathrm{\SIUnitSymbolMicro m}$
($520(6)\text{\,}\mathrm{\SIUnitSymbolMicro m}$) in the center of the vapor
cell by the fiber outcoupling lenses. After the cell further waveplates
linearize and align the polarizations, then a second calcite prism separates
the signal from the control with a polarization extinction ratio of
$>$80\text{\,}\mathrm{dB}$$ (characterization limited by detection) after
which it is fiber coupled to a spectral filtering stage. These spectral
filters consist of three monolithic etalons ($550(10)\text{\,}\mathrm{MHz}$
bandwidth, $25.5(4)\text{\,}\mathrm{GHz}$ free spectral range) in series,
temperature tuned to the signal frequency. For CW light at the control
frequency each of these etalons delivers $-26\text{\,}\mathrm{dB}$ reduction
in intensity. Including the polarization filtering, this results in a total CW
control suppression in the signal channel by more than
$160\text{\,}\mathrm{dB}$. Finally, another fiber coupling after the etalons
sends the signal to the detection system consisting of two single photon
avalanche diodes (SPADs, Excelitas SPCM-AQRH-16) arranged in Hanbury Brown and
Twiss configuration. Despite the meticulous control filtration, the total
transmission of a strong CW probe at the signal frequency from the fiber input
of the memory to this output, leaving the atoms warm and unpumped, is
$T=$30(3)\text{\,}\mathrm{\char 37\relax}$$.

##  IV SPDC Source

The single photon source is based on cavity-enhanced SPDC, which provides
heralded single photons with high quality, high efficiency, and tunable
bandwidth [40, 41, 42]. The source is an evolution of the one described in
[16, 43] with improved reliability and performance. It is carefully tailored
to emit signal photons compatible with the atomic memory, at a wavelength of
$795\text{\,}\mathrm{nm}$ fine-tuned to the ^87Rb D1 line and a bandwidth of
$\sim$ $370\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$ matching the acceptance
bandwidth of the memory.

The heart of the source is a $5\text{\,}\mathrm{m}\mathrm{m}$ long
periodically poled potassium titanyl phosphate (ppKTP) crystal, polished and
coated to form a doubly resonant hemispherical monolithic cavity, see Fig. 2.
A periodic poling of $10.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ is chosen so
that the quasi-phase-matching conditions for type-II SPDC are met for signal
photons wavelength-matched with the ^87Rb D1 line whereas the bandwidth-
matching is given by the cavity linewidth. In the non-degenerate process
$404\text{\,}\mathrm{n}\mathrm{m}$ pump photons are down-converted to
$795\text{\,}\mathrm{n}\mathrm{m}$ ($822\text{\,}\mathrm{n}\mathrm{m}$) signal
(idler) photons, illustrated in Fig. 1 (b). We pump the crystal with
$4.5\text{\,}\mathrm{mW}$ in a double pass configuration to reach heralding
rates of
$1.5\text{\times}{10}^{5}\text{\,}\mathrm{counts}\text{\,}{\mathrm{s}}^{-1}$
on average. The pump frequency is stabilized through a sideband-offset lock on
a passively stable reference cavity, allowing for a tunable but locked laser.
The orthogonally polarized signal and idler photons are split on a polarizing
beam splitter cube, see Fig. 2, and coupled each into a polarization
maintaining fiber. The herald is filtered with a temperature-stabilized
monolithic etalon ($1150(20)\text{\,}\mathrm{MHz}$ bandwidth,
$51(1)\text{\,}\mathrm{GHz}$ free spectral range) and a narrow-band IF
($0.57(5)\text{\,}\mathrm{nm}$ FWHM, measured at $822\text{\,}\mathrm{nm}$)
before being detected with a SPAD.

We observe that a constant noise floor of uncorrelated photons is emitted by
the photon source. In order to suppress this background during photon
retrieval a fast Pockels cell acts as switch for the down-conversion process.
By rotating the polarization of the pump beam by $None$ the down-conversion
process is highly suppressed since the phase-matching conditions are not met
anymore. It takes about $140\text{\,}\mathrm{ns}$ upon the detection of an
idler photon for this switch to turn off the source. As the pump light remains
incident on the crystal no thermal drifts are induced by the switching.
Furthermore, optical isolators are placed both in the signal and idler arms.
They are used to prevent crosstalk from the lasers preparing the atomic state
of the memory and to suppress external cavity modes in the herald path, both
of which originate due to spurious back-reflections on the crystal’s plane
surface. Even with these additional optical elements the heralding efficiency,
that is the probability of having a signal photon exiting the optical fiber to
the memory upon the detection of an idler photon, is as high as
$\eta_{h}=40(4)\%$ for a coincidence window of $6.48\text{\,}\mathrm{ns}$.
Such a high $\eta_{h}$ is crucial as even small amounts of memory read-out
noise can accumulate when no photon is present most of the time [44]. We
measure the conditional second-order autocorrelation of the signal photons at
the memory input to be
$g^{(2)}_{\text{input}}=$4.21(2)\text{\times}{10}^{-2}$$, confirming their
high quality.

##  V Interfacing source and memory

Memories for heralded single photons need to react to the detection of the
herald, in our case the idler photon from the SPDC source. This places
stringent limits on the reaction time of all switching electronics and optics
of the memory setup. The time between the detection of an idler photon and the
arrival of the signal photon in the vapor cell is about
$270\text{\,}\mathrm{ns}$, as a $60\text{\,}\mathrm{m}$ optical fiber is used
to route the signal photons from the source to the memory setup. In a storage
and retrieval experiment the memory is first initialized by the pumping and
repumping lasers for a minimum time of $2\text{\,}\mathrm{\SIUnitSymbolMicro
s}$ to ensure that the desired atomic polarization is achieved. During this
pumping stage the detection of idler photons is rejected on a hardware level.
After this minimum duty cycle of the state preparation, the detection of an
idler photon triggers a digital delay generator (DDG, Highland Technology
T564, $21\text{\,}\mathrm{ns}$ insertion delay, $<$35\text{\,}\mathrm{ps}$$
RMS timing jitter) to switch off both the optical pumping and repumping beams
as well as the pumping of the source until after the photon is retrieved. All
this switching is prioritized to minimize noise, as the extinction ratio of
the switches improves over time, and the cables are kept as short as possible.
This reactive configuration ensures that the memory remains ready to accept a
photon at any time after initialization. The trigger is also relayed to a
second DDG for less time critical tasks, including triggering the generation
of control pulses and time stamping the idler photon detection with a time-to-
digital converter (qutools quTAU, timing resolution $81\text{\,}\mathrm{ps}$).
The typical rate of these experiments is
$1.5\text{\times}{10}^{5}\text{\,}{\mathrm{s}}^{-1}$, set by the chosen
heralding rate.

##  VI Single-photon storage results

Figure 3:  (a) Photon arrival-time histograms for storage experiments showing
the initial leakage of photons through the memory and the retrieval after
$160\text{\,}\mathrm{ns}$ of storage. The data are histogrammed in
$162\text{\,}\mathrm{ps}$ bins. As no intermediate features are visible on a
linear scale the time-axis is broken to show only these peaks. The shaded area
marks the $6.48\text{\,}\mathrm{ns}$ ($80\times$81\text{\,}\mathrm{ps}$$) wide
region of interest for the retrieval for which all figures of merit in the
text are specified. Almost no noise counts are visible on a linear scale. (b)
Memory lifetime as measured by a drop in efficiency and modeled exponentially
as $\eta_{e2e}(\tau_{s})=\eta_{e2e,0}\exp(-\tau_{s}/\tau)$. The shaded region
marks the $1\sigma$ confidence interval of the fit. Each point represents a
storage and retrieval experiment as depicted in (a), where each curve is
integrated for $5\text{\,}\mathrm{min}$. These data were collected a day after
the main results without optimizing the filters anew, and have thus been
scaled to account for a slightly lower setup transmission of about
$26\text{\,}\mathrm{\char 37\relax}$. Figure 4:  Photon arrival-time
histograms for storage experiments on a log-scale. The data are histogrammed
in $1\text{\,}\mathrm{ns}$ ($12\times$81\text{\,}\mathrm{ps}$$) bins, and
$\tau_{p}=0$ represents the time at which an input photon arrives at the
detector if it leaks through the setup without being stored or lost. The
$y$-axis is normalized to the peak of the no read-in curve in the retrieval
window (8325 counts) to yield a proxy for the worst-case signal-to-noise
ratio, and is plotted logarithmically so that small features, labeled A-F, can
be identified. These illustrate the experimental sequence and are described
further in the main text. Colored bars above the plot signify the time windows
in terms of detection time wherein the state preparation beams and the photon
source are effectively turned off.

Photon arrival-time histograms of the photons detected within
$20\text{\,}\mathrm{min}$ of integration time are recorded for three scenarios
and shown in Figs. 3 (a) and 4. In the first scenario repeated attempts of
storing a heralded photon for $160\text{\,}\mathrm{ns}$ are performed, which
we present as a test case for detailed analysis. In a read-out window of
$6.48\text{\,}\mathrm{ns}$ (shaded region) and for a total of
$N_{\text{herald}}=$159\,752\,941$$ storage attempts
$N_{\text{ret}}=$454\,030$$ photons are retrieved.

The other two scenarios represent distinct measurements of the noise
performance. The source blocked scenario represents the typical readout noise
estimation performed for memories, measuring the amount of noise produced when
the input is physically blocked. This serves as a comparative measure to other
memories, characterizing it in isolation. When the source is blocked the
number of memory induced noise counts for the same number of experiments as
above is $N_{\text{noise,}\,\text{mem}}=$29\,075$$. The readout noise floor is
then
$\mu_{\text{mem}}=N_{\text{noise,}\,\text{mem}}/N_{\text{herald}}=$1.82(18)\text{\times}{10}^{-4}$$.
This however does not capture noise stemming either from the source itself in
the form of uncorrelated photons at the signal frequency, nor from spurious
back reflections of light originating from the memory setup in the source.

The total noise is estimated in a measurement omitting the read-in control
pulse, labeled as the no read-in scenario. This induces a small systematic
error as the atomic response to the read-out pulse is also influenced by the
read-in pulse. The no read-in curve therefore delivers a slight overestimation
of the total amount of noise present in the read-out window. The number of
noise counts detected in this curve’s read-out window is
$N_{\text{noise,}\,\text{tot}}=$38\,634$$ yielding a total noise floor of
$\mu_{\text{tot}}=$2.42(24)\text{\times}{10}^{-4}$$. With this we calculate
the end-to-end efficiency defined as

$\eta_{\text{e2e}}=\frac{N_{\text{ret}}-N_{\text{noise,}\,\text{tot}}}{N_{\text{herald}}\eta_{h}\eta_{\text{det}}},$
(1)

to be $\eta_{e2e}=$1.1(2)\text{\,}\mathrm{\char 37\relax}$$. Here
$\eta_{\text{det}}$ is the quantum efficiency of the single photon detectors
specified as $60(6)\text{\,}\mathrm{\char 37\relax}$, which is the dominant
source of uncertainty. The signal-to-noise ratio of the combined experiment is
$\text{SNR}=\frac{N_{\text{ret}}-N_{\text{noise,tot}}}{N_{\text{noise,tot}}}=$10.8(15)$$
which bodes well for the quality of the retrieved photons. Indeed, the
conditioned autocorrelation of the retrieved photons is
$g_{c,\,\text{ret}}^{(2)}=$0.177(23)$$ confirming that the memory emission is
dominated by single photons. This is the first experimental confirmation of
the viability of retrieving single photons from ground-state memories in hot
atomic vapor.

A measurement of the exponential $1/e$ lifetime of the memory is shown in Fig.
3 (b). The fit yields $\tau=$680(50)\text{\,}\mathrm{ns}$$ lifetime and an
initial efficiency of $\eta_{e2e,0}=$1.4(4)\text{\,}\mathrm{\char 37\relax}$$.
This lifetime is close to what is expected for a memory limited by atomic
motion out of the interaction region, but may be slightly shortened by the
lifetime of the initial ground state preparation which has not been
exhaustively optimized for longevity. Accounting for the technical losses in
the setup by dividing out the measured transmission and extrapolating to zero
storage yields a total internal efficiency of
$\eta_{\text{int}}=\eta_{e2e,0}/T=$4.7(14)\text{\,}\mathrm{\char 37\relax}$$.
The conditioned autocorrelation of the retrieved photons eventually increases
for long storage times. For $\tau_{s}=$280\text{\,}\mathrm{ns}$$ it is still
low at $g_{c,\,\text{ret}}^{(2)}=$0.171(28)$$, but by
$\tau_{s}=$700\text{\,}\mathrm{ns}$$ it increases to
$g_{c,\,\text{ret}}^{(2)}=$0.503(93)$$. The lifetime could therefore also be
seen as roughly the time for which $g_{c,\,\text{ret}}^{(2)}\leq 0.5$.

The photon number statistics of the retrieved light depend on the photon
number statistics of the signal pulse, as well as those of the noise and
potentially the process at its origin too. An exact and physically simple
model for the case of an incoherent admixture of noise is derived in [15]. It
treats the noise as being generated independently from the storage and
retrieval processes, the additional light being added to the read-out by
straightforward superposition. It reads

$g_{c,\,\text{ret,}\,\text{theo}}^{(2)}=\frac{\left(N_{\text{ret}}-N_{\text{noise}}\right)^{2}g^{(2)}_{\text{input}}+2N_{\text{noise}}\left(N_{\text{ret}}-N_{\text{noise}}\right)+N_{\text{noise}}^{2}g^{(2)}_{\text{noise}}}{N_{\text{ret}}^{2}}.$
(2)

Here, $g^{(2)}_{\text{input}}$ has been established during the
characterization of the source; the statistics of the noise
$g_{c,\,\text{noise}}^{(2)}$, however, are not measured directly, as
insufficient noise counts accumulate within a reasonable integration time to
evaluate them meaningfully. Sources of noise for which we would expect
$g_{c,\,\text{noise}}^{(2)}=1$ are limited to leaked control laser light.
Known possible thermal noise sources include uncorrelated SPDC photons from
the source, as well as collisional fluorescence and four-wave mixing induced
by the control. We distinguish these scenarios on the memory side by scanning
the final filter etalon with the SPDC source blocked. We observe collisional
fluorescence peaks on resonance with the natural oscillators
$F^{\prime}\rightarrow F=2$ [25] which seem to constitute the entirety of the
noise at the signal frequency. Scanning the etalon to the control frequency,
we also observe the control laser and collisional fluorescence on
$F^{\prime}\rightarrow F=1$, but these are well suppressed when the final
etalon is set to the signal frequency. No further peaks are visible after the
filters, in particular confirming that four-wave mixing is well suppressed by
the memory scheme. We therefore expect $g_{c,\,\text{noise}}^{(2)}=2$ as non-
thermal noise sources are well excluded. Moreover, if the noise can be modeled
incoherently this is a conservative assumption as it represents the worst
case. This approach yields an expected value of
$g_{c,\,\text{ret,}\,\text{theo}}^{(2)}=$0.204(29)$$, which is in excellent
agreement with our observed statistics. In the limits
$g^{(2)}_{\text{input}}\rightarrow 0$, $g_{c,\,\text{noise}}^{(2)}\rightarrow
2$ the exact model reduces to $g_{c,\,\text{ret,}\,\text{theo}}^{(2)}\approx
2/(\text{SNR}+1)$, illustrating the link between the statistical measure and
the memory performance directly visible in the data. Note that were the memory
limited by a coherent noise source such as four-wave mixing gain, incoherent
models would drastically underestimate the $g^{(2)}$ of the retrieved light
[15].

Figure 4 makes small features in the data visible yielding additional insights
into the main features. The peaks at $\tau_{p}=0$ and
$\tau_{p}=$161\text{\,}\mathrm{ns}$$, also displayed in Fig. 3 (a), are now
revealed to be followed by exponential tails. These counts correspond to
fluorescence from the atomic excited state. Due to the high bandwidth, the
majority of this noise is avoided thanks to its temporal separation from the
signal, highlighting an innate advantage of high-speed memories. Further
labeled features correspond to steps in the experimental sequence and
technical effects. In the region A (B) the herald photon has already triggered
the electronics, and we see exponentially decaying fluorescence from the atoms
as the pumping (repumping) beam is switched off. In the region C, after the
photon is stored, the discrepancy between the storage and no read-in curves
corresponds to unintentional read-out of the stored photon. The end-to-end
(internal) efficiency of this unintentional read-out is
$0.38(5)\text{\,}\mathrm{\char 37\relax}$ $($1.26(17)\text{\,}\mathrm{\char
37\relax}$)$. It is caused by ringing in the EOM switching the control and its
limited extinction ratio. At longer storage times the end-to-end efficiency of
this read-out saturates at about $1\text{\,}\mathrm{\char 37\relax}$. The
feature present in both curves right beneath the label C at around
$\tau_{p}=$30\text{\,}\mathrm{ns}$$ is due to afterpulsing of the SPADs. The
steep feature labeled D corresponds to the Pockels cell switching off the SPDC
source. The remaining discrepancy between the no-read and source blocked
curves estimates the secondary effect of the first control pulse on the atoms.
The rising feature labeled E corresponds to both pumping lasers being switched
back on after the retrieval is complete. Finally, the feature labeled F marks
the source switching back on.

##  VII Discussion and Simulation

To analyze limitations of our experiment and moreover guide future development
we simulate the storage and retrieval process numerically. We consider an
atomic four-level system following the description in [45], also taking into
account the transverse profile of both light fields, similar to [46]. We
include both collisional broadening induced by the buffer gas and
inhomogeneous Doppler-broadening of the vapor. The latter is taken into
account by introducing different velocity classes of the atoms, as shown in
[47]. We focus on the scenario of forward retrieval. The storage and the
retrieval processes are much faster than the mean free time between collisions
of about $20\text{\,}\mathrm{ns}$ in the vapor cell used for this experiment.
This allows us to assume that the individual atoms do not change their
velocities during these processes. However, during the storage time of
$160\text{\,}\mathrm{ns}$ we assume that the atoms rethermalize fully. The
numerical simulation is implemented by solving the temporal derivatives
appearing in the equations of motion with a partially implicit second-order
Runge-Kutta method. For the spatial derivatives spectral collocation is used,
allowing us to replace them with Chebyshev differentiation matrices as
described in [46].

We use the simulation to evaluate the total efficiency as a function of the
peak control Rabi frequency. Herein the temporal width of the Gaussian control
pulse is fixed to $3.77\text{\,}\mathrm{ns}$, and the time alignment between
signal and control pulse is optimized for each simulation point. Seeking
plausible routes to improvements, other experimentally controllable parameters
are varied. We find that the two-photon detuning has a strong impact on the
attainable performance, as can be seen in Fig. 5 (a). Minimization of the two-
photon detuning intuitively corresponds to bringing the memory interaction
into resonance, note however, that following the convention of [48] we define
this parameter as the difference between the control and signal detunings on
their respective single-photon transitions,
$\Delta_{\text{tp}}=\Delta_{c}-\Delta_{s}$. In the regime of our experiment
where $\Omega\approx\Delta<\Delta_{\text{hfs}}$, the presence of the control
pulse induces a significant, time-dependent level-shift, which by this
definition appears as a shift of the optimal detuning from
$\Delta_{\text{tp}}=0$. Furthermore, different control beam waists are
considered. The simulation confirms that by choosing a larger control beam
waist the total efficiency can be improved. Fig. 5 (b) shows the performance
for different control waist for optimal two-photon detuning. With a control
beam significantly larger than the signal, more atoms experience the optimal
control Rabi frequency, leading to a more efficient process.

Figure 5:  (a) Simulation of the internal memory efficiency using the known
experimental parameters, as a function of the peak control Rabi frequency for
different two-photon detunings. (b) Simulation of the internal memory
efficiency as above, but for different control beam waists and setting
$\Delta_{\text{tp}}/2\pi=$-130\text{\,}\mathrm{MHz}$$, where
$w_{c}=$260(3)\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ is the experimental
waist of the control in the vapor cell center.

For the parameters of the experiment shown in Figs. 3 (a) and 4 the measured
value of the internal efficiency falls short of that predicted by our model.
This leads us to identify effects not captured by the simulation. One is the
unintentional read-out of the signal during the storage time by the leaking
control light which is present due to a combination of ringing in the EOM as
well as its finite extinction ratio (about $25\text{\,}\mathrm{dB}$). A second
issue is the quality of the temporal alignment of the control pulse to the
signal, mainly limited by the $>$350\text{\,}\mathrm{ps}$$ jitter of the
single photon detector measuring the idler photon. Based on tests with laser
pulses where this jitter source can be eliminated, we expect this amount of
variance in the temporal alignment to decrease the efficiency by a factor of
$<0.90(5)$. These issues, alongside the difficulty of stating an experimental
value of $\Delta_{\text{tp}}$ with any certainty, makes modeling the exact
experimental situation difficult. Experimentally, the signal frequency is set
once, whereupon the control frequency is optimized on the storage efficiency
at fixed $\Omega$ to compensate for drifts, without the possibility of a
direct measurement of the two-photon detuning. Our simulations will guide
future improvements by identifying the regimes were the best performance is
expected, accounting for the interdependence of the many optimization
parameters. This analysis highlights three main areas where technical
improvements can still be implemented: lasers, optical switches, and
detectors. A more powerful control laser would enable us to define a larger
and more homogeneous interaction region. Optical switches without ringing and
with a better on/off ratio would protect the spinwave during storage. Faster
single photon detectors would improve the time alignment of the control pulse
to the single signal photon. In particular, the currently remaining jitter on
a control pulse, removing the jitter of the idler detection, is measured to be
only $126(3)\text{\,}\mathrm{ps}$. The state-of-the-art in timing resolution
of single photon detection is more than an order of magnitude smaller.

Upon addressing technical matters, a further improvement of the efficiency is
possible by increasing the optical depth. Technically this is straightforward
via the atomic temperature, but is bound to affect other figures of merit as
well, in particular the readout noise. Due to their interdependence,
parameters such as the detuning from the excited states and Rabi frequency
would again require optimization, and accurate performance predictions
concerning noise, fidelity, and lifetime predicate dedicated models not
currently implemented. Additionally, some gains could be expected from shaping
the temporal profile of the control pulses. This is fraught with ambitious
demands on the time resolution of the pulse generation. A more straightforward
change, still requiring some modifications to the electronic systems though,
would be to optimize the amplitude of the read-in and read-out pulses
independently. Although the optimal processes are well known to be linked by
time-reversal symmetry both theoretically [49] and experimentally [50], the
matter is less cut-and-dried in the experimentally implemented high-bandwidth
case with realizable Gaussian control pulses, where a full pulse shape
optimization is technically prohibitive and a practical distinction between
programmed and effective pulses is warranted.

##  VIII Summary and Outlook

In summary, we have for the first time demonstrated the interfacing of a
heralded single photon source and a technologically simple atomic quantum
memory with ground state storage successfully maintaining the single photon
nature of the input upon retrieval. Through simulations we have laid a road
map for future improvements that will simultaneously realize state-of-the-art
efficiency. Moreover, the fidelity of the memory will be determined more
directly with Hong-Ou-Mandel experiments [51]. Modern vapor cell fabrication
techniques can produce cells with diameters $<$1\text{\,}\mathrm{mm}$$ [52].
Even with the currently available control power, which determines the mode
size yielding the optimal control intensity, this would match the cell size to
the optical mode. Under such conditions the atoms would remain in the
interaction region despite their motion during storage. With corresponding
improvements to the longevity of the Zeeman state preparation, which is known
to be short but not currently limiting, a memory lifetime on the order of
milliseconds is well within reach judging by similar systems [14]. Moreover,
Zeeman-state preserving anti-relaxation coatings on cells have shown far
longer relaxation times [53]. Possible applications of our source-memory
system thus include photon synchronization and use in high-bandwidth quantum
networks for quantum information processing and simulation. The multi-photon
rate enhancement provided by memories is given by the time-bandwidth product
$B$, scaled by the memory efficiency [6]. Due to the high bandwidth operation
of our memory we achieve $B=$250(20)$$ despite the non-optimized storage time,
which would be sufficient to demonstrate enhanced synchronized rates directly.

###### Acknowledgements.

We thank C. Müller, T. Kroh, A. Ahlrichs, S. Ramelow, and O. Benson for
initiating us into the mysteries of SPDC sources. We thank R. Warburton, L.
Zhai, G. Nguyen, and C. Spinnler for fruitful discussions and acknowledge
financial support from NCCR QSIT.

## References

  * Kimble [2008] H. J. Kimble, The quantum internet, Nature 453, 1023 (2008).
  * Wehner  _et al._ [2018] S. Wehner, D. Elkouss, and R. Hanson, Quantum internet: A vision for the road ahead, Science 362, eaam9288 (2018).
  * Briegel  _et al._ [1998] H.-J. Briegel, W. Dür, J. I. Cirac, and P. Zoller, Quantum repeaters: The role of imperfect local operations in quantum communication, Physical Review Letters 81, 5932 (1998).
  * Sangouard  _et al._ [2011] N. Sangouard, C. Simon, H. de Riedmatten, and N. Gisin, Quantum repeaters based on atomic ensembles and linear optics, Reviews of Modern Physics 83, 33 (2011).
  * Munro  _et al._ [2015] W. J. Munro, K. Azuma, K. Tamaki, and K. Nemoto, Inside quantum repeaters, IEEE Journal of Selected Topics in Quantum Electronics 21, 78 (2015).
  * Nunn  _et al._ [2013] J. Nunn, N. K. Langford, W. S. Kolthammer, T. F. M. Champion, M. R. Sprague, P. S. Michelberger, X.-M. Jin, D. G. England, and I. A. Walmsley, Enhancing multiphoton rates with quantum memories, Physical Review Letters 110, 133601 (2013).
  * Kaneda  _et al._ [2017] F. Kaneda, F. Xu, J. Chapman, and P. G. Kwiat, Quantum-memory-assisted multi-photon generation for efficient quantum information processing, Optica 4, 1034 (2017).
  * Gao  _et al._ [2019] S. Gao, O. Lazo-Arjona, B. Brecht, K. T. Kaczmarek, S. E. Thomas, J. Nunn, P. M. Ledingham, D. J. Saunders, and I. A. Walmsley, Optimal coherent filtering for single noisy photons, Physical Review Letters 123, 213604 (2019).
  * O’Brien [2007] J. L. O’Brien, Optical quantum computing, Science 318, 1567 (2007).
  * Bussières  _et al._ [2013] F. Bussières, N. Sangouard, M. Afzelius, H. de Riedmatten, C. Simon, and W. Tittel, Prospective applications of optical quantum memories, Journal of Modern Optics 60, 1519 (2013).
  * Simon  _et al._ [2010] C. Simon, M. Afzelius, J. Appel, A. B. de la Giroday, S. J. Dewhurst, N. Gisin, C. Y. Hu, F. Jelezko, S. Kröll, J. H. Müller, J. Nunn, E. S. Polzik, J. G. Rarity, H. D. Riedmatten, W. Rosenfeld, A. J. Shields, N. Sköld, R. M. Stevenson, R. Thew, I. A. Walmsley, M. C. Weber, H. Weinfurter, J. Wrachtrup, and R. J. Young, Quantum memories, The European Physical Journal D 58, 1 (2010).
  * Namazi  _et al._ [2017] M. Namazi, C. Kupchak, B. Jordaan, R. Shahrokhshahi, and E. Figueroa, Ultralow-noise room-temperature quantum memory for polarization qubits, Physical Review Applied 8, 034023 (2017).
  * Borregaard  _et al._ [2016] J. Borregaard, M. Zugenmaier, J. M. Petersen, H. Shen, G. Vasilakis, K. Jensen, E. S. Polzik, and A. S. Sørensen, Scalable photonic network architecture based on motional averaging in room temperature gas, Nature Communications 7, 11356 (2016).
  * Katz and Firstenberg [2018] O. Katz and O. Firstenberg, Light storage for one second in room-temperature alkali vapor, Nature Communications 9, 2074 (2018).
  * Michelberger  _et al._ [2015] P. S. Michelberger, T. F. M. Champion, M. R. Sprague, K. T. Kaczmarek, M. Barbieri, X. M. Jin, D. G. England, W. S. Kolthammer, D. J. Saunders, J. Nunn, and I. A. Walmsley, Interfacing GHz-bandwidth heralded single photons with a warm vapour Raman memory, New Journal of Physics 17, 043006 (2015).
  * Mottola  _et al._ [2020] R. Mottola, G. Buser, C. Müller, T. Kroh, A. Ahlrichs, S. Ramelow, O. Benson, P. Treutlein, and J. Wolters, An efficient, tunable, and robust source of narrow-band photon pairs at the 87Rb D1 line, Optics Express 28, 3159 (2020).
  * Akopian  _et al._ [2011] N. Akopian, L. Wang, A. Rastelli, O. G. Schmidt, and V. Zwiller, Hybrid semiconductor-atomic interface: Slowing down single photons from a quantum dot, Nature Photonics 5, 230 (2011).
  * Ulrich  _et al._ [2014] S. M. Ulrich, S. Weiler, M. Oster, M. Jetter, A. Urvoy, R. Löw, and P. Michler, Spectroscopy of the D1 transition of cesium by dressed-state resonance fluorescence from a single (In,Ga)As/GaAs quantum dot, Physical Review B 90, 125310 (2014).
  * Jahn  _et al._ [2015] J.-P. Jahn, M. Munsch, L. Béguin, A. V. Kuhlmann, M. Renggli, Y. Huo, F. Ding, R. Trotta, M. Reindl, O. G. Schmidt, A. Rastelli, P. Treutlein, and R. J. Warburton, An artificial Rb atom in a semiconductor with lifetime-limited linewidth, Physical Review B 92, 245439 (2015).
  * Zhai  _et al._ [2020] L. Zhai, M. C. Löbl, J.-P. Jahn, Y. Huo, P. Treutlein, O. G. Schmidt, A. Rastelli, and R. J. Warburton, Large-range frequency tuning of a narrow-linewidth quantum emitter, Applied Physics Letters 117, 083106 (2020).
  * Gündoğan  _et al._ [2021] M. Gündoğan, J. S. Sidhu, V. Henderson, L. Mazzarella, J. Wolters, D. K. L. Oi, and M. Krutzik, Proposal for space-borne quantum memories for global quantum networking, npj Quantum Information 7, 10.1038/s41534-021-00460-9 (2021).
  * Wallnöfer  _et al._ [2021] J. Wallnöfer, F. Hahn, M. Gündoğan, J. S. Sidhu, F. Krüger, N. Walk, J. Eisert, and J. Wolters, Simulating quantum repeater strategies for multiple satellites, (2021), arXiv:2110.15806 [quant-ph] .
  * Phillips  _et al._ [2011] N. B. Phillips, A. V. Gorshkov, and I. Novikova, Light storage in an optically thick atomic ensemble under conditions of electromagnetically induced transparency and four-wave mixing, Physical Review A 83, 063823 (2011).
  * Lauk  _et al._ [2013] N. Lauk, C. O’Brien, and M. Fleischhauer, Fidelity of photon propagation in electromagnetically induced transparency in the presence of four-wave mixing, Physical Review A 88, 013823 (2013).
  * Rousseau  _et al._ [1975] D. L. Rousseau, G. D. Patterson, and P. F. Williams, Resonance Raman scattering and collision-induced redistribution scattering in I2, Physical Review Letters 34, 1306 (1975).
  * Manz  _et al._ [2007] S. Manz, T. Fernholz, J. Schmiedmayer, and J.-W. Pan, Collisional decoherence during writing and reading quantum states, Physical Review A 75, 040101(R) (2007).
  * Chanelière  _et al._ [2005] T. Chanelière, D. N. Matsukevich, S. D. Jenkins, S.-Y. Lan, T. A. B. Kennedy, and A. Kuzmich, Storage and retrieval of single photons transmitted between remote quantum memories, Nature 438, 833 (2005).
  * Choi  _et al._ [2008] K. S. Choi, H. Deng, J. Laurat, and H. J. Kimble, Mapping photonic entanglement into and out of a quantum memory, Nature 452, 67 (2008).
  * Zhou  _et al._ [2012] S. Zhou, S. Zhang, C. Liu, J. F. Chen, J. Wen, M. M. T. Loy, G. K. L. Wong, and S. Du, Optimal storage and retrieval of single-photon waveforms, Optics Express 20, 24124 (2012).
  * Wang  _et al._ [2019] Y. Wang, J. Li, S. Zhang, K. Su, Y. Zhou, K. Liao, S. Du, H. Yan, and S.-L. Zhu, Efficient quantum memory for single-photon polarization qubits, Nature Photonics 13, 346 (2019).
  * Kaczmarek  _et al._ [2018] K. T. Kaczmarek, P. M. Ledingham, B. Brecht, S. E. Thomas, G. S. Thekkadath, O. Lazo-Arjona, J. H. D. Munns, E. Poem, A. Feizpour, D. J. Saunders, J. Nunn, and I. A. Walmsley, High-speed noise-free optical quantum memory, Physical Review A 97, 042316 (2018).
  * Finkelstein  _et al._ [2018] R. Finkelstein, E. Poem, O. Michel, O. Lahad, and O. Firstenberg, Fast, noise-free memory for photon synchronization at room temperature, Science Advances 4, eaap8598 (2018).
  * Wolters  _et al._ [2017] J. Wolters, G. Buser, A. Horsley, L. Béguin, A. Jöckel, J.-P. Jahn, R. J. Warburton, and P. Treutlein, Simple atomic quantum memory suitable for semiconductor quantum dot single photons, Physical Review Letters 119, 060502 (2017).
  * Thomas  _et al._ [2019] S. E. Thomas, T. M. Hird, J. H. D. Munns, B. Brecht, D. J. Saunders, J. Nunn, I. A. Walmsley, and P. M. Ledingham, Raman quantum memory with built-in suppression of four-wave-mixing noise, Physical Review A 100, 033801 (2019).
  * Fleischhauer and Lukin [2002] M. Fleischhauer and M. D. Lukin, Quantum memory for photons: Dark-state polaritons, Physical Review A 65, 022314 (2002).
  * Yan  _et al._ [2001] M. Yan, E. G. Rickey, and Y. Zhu, Electromagnetically induced transparency in cold rubidium atoms, Journal of the Optical Society of America B 18, 1057 (2001).
  * Walther  _et al._ [2007] P. Walther, M. D. Eisaman, A. André, F. Massou, M. Fleischhauer, A. S. Zibrov, and M. D. Lukin, Generation of narrow-bandwidth single photons using electromagnetically induced transparency in atomic ensembles, International Journal of Quantum Information 05, 51 (2007).
  * Vurgaftman and Bashkansky [2013] I. Vurgaftman and M. Bashkansky, Suppressing four-wave mixing in warm-atomic-vapor quantum memory, Physical Review A 87, 063836 (2013).
  * Rosenberry  _et al._ [2007] M. A. Rosenberry, J. P. Reyes, D. Tupa, and T. J. Gay, Radiation trapping in rubidium optical pumping at low buffer-gas pressures, Physical Review A 75, 023401 (2007).
  * Slattery  _et al._ [2015] O. Slattery, L. Ma, P. Kuo, and X. Tang, Narrow-linewidth source of greatly non-degenerate photon pairs for quantum repeaters from a short singly resonant cavity, Applied Physics B 121, 413 (2015).
  * Ahlrichs and Benson [2016] A. Ahlrichs and O. Benson, Bright source of indistinguishable photons based on cavity-enhanced parametric down-conversion utilizing the cluster effect, Applied Physics Letters 108, 021111 (2016).
  * Tsai and Chen [2018] P.-J. Tsai and Y.-C. Chen, Ultrabright, narrow-band photon-pair source for atomic quantum memories, Quantum Science and Technology 3, 034005 (2018).
  * Cotting [2021] B. Cotting,  _Spontaneous parametric down-conversion heralded single-photon source for quantum memory applications_ , Master’s thesis, University of Basel, École polytechnique fédérale de Lausanne (2021).
  * Jobez  _et al._ [2015] P. Jobez, C. Laplane, N. Timoney, N. Gisin, A. Ferrier, P. Goldner, and M. Afzelius, Coherent spin control at the quantum level in an ensemble-based optical memory, Physical Review Letters 114, 230502 (2015).
  * Rakher  _et al._ [2013] M. T. Rakher, R. J. Warburton, and P. Treutlein, Prospects for storage and retrieval of a quantum-dot single photon in an ultracold 87Rb ensemble, Physical Review A 88, 053834 (2013).
  * Nunn [2008] J. Nunn,  _Quantum memory in atomic ensembles_ , Ph.D. thesis, St. John’s College, University of Oxford (2008).
  * Gorshkov  _et al._ [2007a] A. V. Gorshkov, A. André, M. D. Lukin, and A. S. Sørensen, Photon storage in $\Lambda$-type optically dense atomic media. III. Effects of inhomogeneous broadening, Physical Review A 76, 033806 (2007a).
  * Gorshkov  _et al._ [2007b] A. V. Gorshkov, A. André, M. Fleischhauer, A. S. Sørensen, and M. D. Lukin, Universal approach to optimal photon storage in atomic media, Physical Review Letters 98, 123601 (2007b).
  * Gorshkov  _et al._ [2008] A. V. Gorshkov, T. Calarco, M. D. Lukin, and A. S. Sørensen, Photon storage in $\Lambda$-type optically dense atomic media. IV. optimal control using gradient ascent, Phys. Rev. A 77, 043806 (2008).
  * Phillips  _et al._ [2008] N. B. Phillips, A. V. Gorshkov, and I. Novikova, Optimal light storage in atomic vapor, Physical Review A 78, 023801 (2008).
  * Hong  _et al._ [1987] C. K. Hong, Z. Y. Ou, and L. Mandel, Measurement of subpicosecond time intervals between two photons by interference, Physical Review Letters 59, 2044 (1987).
  * Møller  _et al._ [2017] C. B. Møller, R. A. Thomas, G. Vasilakis, E. Zeuthen, Y. Tsaturyan, M. Balabas, K. Jensen, A. Schliesser, K. Hammerer, and E. S. Polzik, Quantum back-action-evading measurement of motion in a negative mass reference frame, Nature 547, 191 (2017).
  * Corsini  _et al._ [2013] E. P. Corsini, T. Karaulanov, M. Balabas, and D. Budker, Hyperfine frequency shift and Zeeman relaxation in alkali-metal-vapor cells with antirelaxation alkene coating, Physical Review A 87, 022901 (2013).


# Analytical Solutions for Radiation-Driven Winds in Massive Stars II: The
$\delta$-slow Regime

I. Araya1, A. Christen2, M. Curé3, L. S. Cidale4,5, R. O. J. Venero4,5, C.
Arcos3, A. C. Gormaz-Matamala3,7, M. Haucke4,5, P. Escárate 8 , and H.
Clavería9
1 Centro de Investigación DAiTA Lab, Facultad de Estudios Interdisciplinarios,
Universidad Mayor, Chile
2 Instituto de Estadística, Facultad de Ciencias, Universidad de Valparaíso,
Av. Gran Bretaña 1111, Valparaíso, Chile
3 Instituto de Física y Astronomía, Facultad de Ciencias, Universidad de
Valparaíso, Av. Gran Bretaña 1111, Valparaíso, Chile
4 Departamento de Espectroscopía, Facultad de Ciencias Astronómicas y
Geofísicas, Universidad Nacional de La Plata (UNLP), Argentina
5 Instituto de Astrofísica La Plata, CCT La Plata, CONICET-UNLP, La Plata,
Argentina
7 Departamento de Ciencias, Facultad de Artes Liberales, Universidad Adolfo
Ibáñez, Av. Padre Hurtado 750, Viña del Mar, Chile
8 Instituto de Electricidad y Electrónica, Facultad de Ciencias de la
Ingeniería, Universidad Austral de Chile, General Lagos 2086, Valdivia, Chile
9 Instituto de Estadística, Pontificia Universidad Católica de Valparaíso,
Chile E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Accurate mass-loss rates and terminal velocities from massive stars winds are
essential to obtain synthetic spectra from radiative transfer calculations and
to determine the evolutionary path of massive stars. From a theoretical point
of view, analytical expressions for the wind parameters and velocity profile
would have many advantages over numerical calculations that solve the complex
non-linear set of hydrodynamic equations. In a previous work, we obtained an
analytical description for the fast wind regime. Now, we propose an
approximate expression for the line-force in terms of new parameters and
obtain a velocity profile closed-form solution (in terms of the Lambert $W$
function) for the $\delta$-slow regime. Using this analytical velocity
profile, we were able to obtain the mass-loss rates based on the m-CAK theory.
Moreover, we established a relation between this new set of line-force
parameters with the known stellar and m-CAK line-force parameters. To this
purpose, we calculated a grid of numerical hydrodynamical models and performed
a multivariate multiple regression. The numerical and our descriptions lead to
good agreement between their values.

###### keywords:

hydrodynamics — methods: analytical — stars: early-type — stars: mass-loss —
stars: winds, outflows

††pubyear: 2021††pagerange: Analytical Solutions for Radiation-Driven Winds in
Massive Stars II: The $\delta$-slow Regime–A

##  1 Introduction

The knowledge of stellar wind properties of massive stars is fundamental for
understanding stellar evolution processes, different evolutionary scenarios
and enrichment of star’s nearby environments.

Accurate wind parameters (mass-loss rate and terminal velocity) are crucial
for the study of the wind properties of massive stars. Insights into the
physics of stellar winds are attained by studying the effects of wind
parameters on the emergent line spectrum and by comparing the latter with
observations. From a theoretical point of view, this implies to solve highly
non-linear equations in which the radiation field and hydrodynamics are
strongly coupled.

Winds of massive stars are driven by the transfer of momentum from the
radiation field to the plasma by scattering processes in the spectral lines
(Lucy & Solomon, 1970). Currently, these winds are best described by the m-CAK
theory (Castor et al., 1975; Friend & Abbott, 1986; Pauldrach et al., 1986).

Generally, there are many approximations that reduce considerably the
complexity of the computation of the hydrodynamic and the NLTE radiative
transfer solutions. One example is the extensive use of a simple analytical
approximation for the velocity field, the so-called $\beta$-law, first
proposed by Lamers & Rogerson (1978). A value of $\beta\simeq 0.8-1.2$,
generally agrees very well with the m-CAK numerical hydrodynamic solution
(Lamers & Cassinelli, 1999). This value of $\beta$ is determined empirically
by fitting the observed line profile with a synthetic one. This approximation
has been proved to be very effective and efficient to describe the winds of O-
and early B-type supergiants. However, in the case of late B- and A-type
supergiants there is a clear tendency towards higher values of $\beta$, even
with values larger than $3$, leading to inconsistencies with respect to the
hydrodynamic theory (Stahl et al., 1991; Verdugo et al., 1999; Crowther et
al., 2006; Lefever et al., 2007; Markova & Puls, 2008; Searle et al., 2008;
Haucke et al., 2018). Therefore, accurate analytical approximations of the
m-CAK hydrodynamic equations are indispensable to have a self-consistent
coupling between the hydrodynamics and multidimensional radiative transfer
problems in moving media.

For the case of the fast regime (standard m-CAK solution), this issue was
addressed by Villata (1992), Müller & Vink (2008) and Araya et al. (2014). The
aim of this work is to extend the procedure of Araya et al. (2014) to the
$\delta$-slow regime. The $\delta$-slow solution 111Previously, Curé (2004)
found another type of slow solution for rapidly rotating stars, called
$\Omega$-slow solution., found by Curé et al. (2011), is based on the m-CAK
theory, that describes the wind velocity profile when the ionization-related
line-force parameter $\delta$ takes higher values than the ones provided by
the standard m-CAK solution (see, e.g., Lamers & Cassinelli, 1999, and
references therein). High values of $\delta$, even larger than $1/3$, which
corresponds to a wind with neutral hydrogen as a trace element (Puls et al.,
2000), are expected in strong ionization gradients (see also Kudritzki, 2002).
The $\delta$-slow solution is characterized by low terminal speeds
($v_{\infty}$) and might explain the obtained values for late-B and A-type
supergiants. This solution also seems to fit quite well the observed anomalous
correlation between the terminal and escape velocities found in A supergiants,
as well as their corresponding wind momentum-luminosity relationship (Curé et
al., 2011).

With the purpose to have an approximate solution from the hydrodynamic, Araya
et al. (2014) developed an expression in terms of the stellar and m-CAK line-
force parameters ($\alpha$, $k$, and $\delta$) and applied it to the fast
regime. This expression, based on the works of Müller & Vink (2008) and
Villata (1992), describes the line acceleration as function of the radial
distance, allowing to solve analytically the hydrodynamic stationary equation
of motion. The use of expressions for both radiation force and velocity
profile as a function of the line-force parameters can provide a clear view
into how the line-driven mechanism is related with the hydrodynamics.

On the other hand, it is important to obtain a simple representation of the
radiation force and the derived slow solutions under such different ionization
conditions. Therefore, a significant contribution of this work consists in
offering a quick way to generate an analytical expression to estimate mass-
loss rates for these alternative wind regimes. There are currently no
parametric expressions that can be used for this purpose without the need to
fully solve the hydrodynamic equations.

This work is organized as follows: Section 2 presents briefly the hydrodynamic
equations for line-driven winds and the dimensionless form of the equation of
motion. In Section 3, the basic concepts developed by Müller & Vink (2008) are
recapitulated including their line acceleration term as function of the radial
distance. Then, this line acceleration term is modified with the purpose to
obtain a better agreement with the $\delta$-slow solution. In Section 4, a
recipe to obtain the line acceleration parameters (required by the line
acceleration term) is developed, based on a grid of hydrodynamic models and a
multivariate multiple regression. Then, an analytical expression for the
$\delta$-slow solution is developed and compared with the numerical models
described in Section 5. In Section 6, we give our conclusions. In addition, a
recipe to derive the mass-loss rate based on our expression is provided in
Appendix A.

##  2 The Standard Hydrodynamical Wind Model

The CAK theory for line-driven winds was originally developed by Castor et al.
(1975). This theory describes, for a point source, a stationary, one-
dimensional, non-rotating, isothermal, outflowing wind with spherical
symmetry. Adopting these assumptions, and neglecting the effects of viscosity,
heat conduction and magnetic fields, the equations of mass conservation and
radial momentum state:

$4\,\pi\,r^{2}\,\rho\,v=\dot{M},$ (1)

and

$v\,\frac{dv}{dr}=-\frac{1}{\rho}\frac{dp}{dr}-\frac{G\,M_{*}(1-\Gamma_{\rm{E}})}{r^{2}}+g^{\rm{line}}.$
(2)

Here $v$ is the fluid radial velocity, $dv/dr=v^{\prime}$ is the velocity
gradient and $g^{\rm{line}}$ is the line acceleration. All other variables
have their standard meaning (for a detailed derivation and definitions of
variables, constants and functions, see Curé, 2004).

The so called m-CAK theory, which include the effects of rotation and a disk-
like source, was developed by Friend & Abbott (1986) and Pauldrach et al.
(1986), based on a general expression from Abbott (1982) for the line force:

$g^{\rm{line}}=\frac{C}{r^{2}}\,f_{\rm{FD}}(r,v,v^{\prime})\,\left(r^{2}\,v\,v^{\prime}\right)^{\alpha}\left(\frac{n_{E11}}{W(r)}\right)^{\delta},$
(3)

where the coefficient $C$ (eigenvalue) depends on the mass-loss rate $\dot{M}$
and the line-force parameter $k$ (see Eq. 25). $W(r)$ is the dilution factor,
$n_{E11}$ is the electron number density $n_{E}$ in units of
$10^{-11}\,\rm{cm^{-3}}$, and $f_{\rm{FD}}$ is the finite disk correction
factor. The m-CAK line-force parameters are: $\alpha$ , $k$ and $\delta$.

The momentum equation (Eq. 2) can be expressed in a dimensionless form (see
e.g., Müller & Vink, 2008; Araya et al., 2014) as:

$\hat{v}\,\frac{d\hat{v}}{d\hat{r}}=-\frac{\hat{v}_{\rm{crit}}^{2}}{\hat{r}^{2}}+\hat{g}^{\rm{line}}-\frac{1}{\rho}\frac{d\rho}{d\hat{r}},$
(4)

with $\hat{r}=r/R_{*}$, $\hat{v}=v/a$ and
$\hat{v}_{\rm{crit}}=v_{\rm{esc}}/a\sqrt{2}$. Here $R_{*}$ is the stellar
radius, $a$ is the isothermal sound speed, $\hat{v}_{\rm{crit}}$ is the
dimensionless rotational break-up velocity and $v_{\rm{esc}}$ is the escape
velocity. The dimensionless line acceleration reads:

$\hat{g}^{\rm{line}}=\frac{R_{*}}{a^{2}}\,g^{\rm{line}}.$ (5)

Using Eq. 1 together with the equation of state for an ideal gas
($p=a^{2}\rho$), the dimensionless equation of motion is:

$\left(\hat{v}-\frac{1}{\hat{v}}\right)\frac{d\hat{v}}{d\hat{r}}=-\frac{\hat{v}_{\rm{crit}}^{2}}{\hat{r}^{2}}+\frac{2}{\hat{r}}+\hat{g}^{\rm{line}}.$
(6)

In general, the calculation of the line acceleration involves the coupling of
hydrodynamics with the radiative transport in NLTE. A very successful approach
is to calculate the line acceleration using the Sobolev approximation. The
pioneering work of Castor et al. (1975) laid the foundations of CAK theory and
later improvements (m-CAK). A further description was done by Feldmeier (1998)
who extended the CAK approach using a second order Sobolev approximation, i.e,
$g^{\rm{line}}=g^{\rm{line}}(r,v,v^{\prime},v^{\prime\prime})$.
However, in this work, to obtain an analytical expression of the $\delta$-slow
solution, we will use a radial dependence for the line acceleration following
the methodology used by Araya et al. (2014), i.e.,
$g^{\rm{line}}=g^{\rm{line}}(r)$. This approach allows to obtain an analytical
expression for the velocity field in terms of the Lambert $W$ function (see
Section 3).

##  3 Line Acceleration

In this section we review the basic concepts developed by Müller & Vink (2008,
hereafter MV08) to derive, later on, a general analytical expression for the
velocity profile in the frame of the $\delta$-slow radiation-driven wind
regime for massive stars. We demonstrate that this expression enables to
integrate the equation of motion (Eq. 6) leading to an analytical expression
for the $\delta$-slow wind velocity profile.

###  3.1 The Fast Regime Approximation

In the framework of m-CAK stellar wind theory, MV08 present a mathematical
expression for the line acceleration via a parameterized description that
depends only on the radial coordinate. Using Monte Carlo multi-line radiative
transfer calculations (de Koter et al., 1997; Vink et al., 1999) and a
velocity profile from a $\beta$-law, these authors computed the line
acceleration. Then, the numerical line acceleration, which collect all the
physically motivated mathematical properties for the radiative line
acceleration term, is expressed by by the following function:

$\hat{g}^{\rm{line}}_{\mathrm{MV08}}(\hat{r})=\frac{\hat{g}_{0}}{\hat{r}^{1+\delta_{1}}}\left(1-\frac{\hat{r_{0}}}{\hat{r}^{\delta_{1}}}\right)^{\gamma},$
(7)

where $\hat{g_{0}}$, $\delta_{1}$, $\hat{r_{0}}$, and $\gamma$ are the MV08
line acceleration parameters. It is important to note that these parameters,
lack of any physical meaning, and besides, are not directly related to $k$,
$\alpha$ and $\delta$ parameters from m-CAK theory.

Replacing Eq. 7 in Eq. 6, the dimensionless equation of motion are derived and
a fully analytical velocity profile is obtained (see MV08 for details about
the methodology used to obtain this solution) by means of the Lambert
W-function (Corless et al., 1993; Corless et al., 1996; Cranmer, 2004).

The line acceleration expression given by MV08 (Eq. 7) results in a good
approximation for the m-CAK line force for $\delta\leq 0.2$, but this
expression fails for $\delta$-slow solutions, when $\delta\gtrsim 0.25$.
Overall, this approximation gives a poor agreement with respect to the
numerical $\delta$-slow solution (from m-CAK theory). The numerical solutions
are obtained from the stationary hydrodynamic code Hydwind (Curé, 2004).

Araya et al. (2014) developed an analytical solution for the velocity of the
fast wind regime in terms of the stellar and m-CAK line-force parameters
combining the methodology from MV08 and the line acceleration proposed by
Villata (1992). Unfortunately, this expression also fails when the line force
parameter $\delta$ is higher than about $0.3$, because in this case a term
from the proposed line acceleration expression turns complex. From a
mathematical point of view, high values of $\delta$ would require high values
of $\alpha$ in order to obtain an expression with real values, but such kind
of $\alpha$ values would be totally unphysical.

###  3.2 The New $\delta$-slow Regime Approximation

In view of the unsatisfactory results obtained when applying the approximate
description of the wind velocity for the $\delta$-slow case, we decided to
modify the functional form of the line acceleration given by MV08 in order to
obtain a better description of the $\delta$-slow wind. Thus, our proposed line
acceleration is the following:

$\hat{g}^{\rm{line}}_{\rm{new}}(\hat{r})=\frac{\hat{g}_{0}}{\hat{r}^{1+\delta_{1}}}\left(1-\frac{1}{\hat{r}^{\delta_{2}}}\right)^{\gamma},$
(8)

where $\hat{g}_{0}$, $\delta_{1}$, $\delta_{2}$, and $\gamma$ are the new set
of line acceleration parameters.

The new expression follows the same mathematical properties as MV08’s but the
inclusion of the $\delta_{2}$ parameter yields to a better agreement with the
numerical line acceleration from the m-CAK model.

Based on this new definition for the radiation force, the new dimensionless
equation of motion reads:

$\left(\hat{v}-\frac{1}{\hat{v}}\right)\frac{d\hat{v}}{d\hat{r}}=-\frac{\hat{v}_{{\rm
crit}}^{2}}{\hat{r}^{2}}+\frac{2}{\hat{r}}+\frac{\hat{g}_{0}}{\hat{r}^{1+\delta_{1}}}\left(1-\frac{1}{\hat{r}^{\delta_{2}}}\right)^{\gamma}.$
(9)

The same methodology developed by MV08 is employed to solve the new equation
of motion and the solution is given through the Lambert $W$ function,

$\hat{v}(\hat{r})=\sqrt{-W_{j}(x(\hat{r}))},$ (10)

with

$\displaystyle x(\hat{r})$ $\displaystyle=$
$\displaystyle-\left(\frac{\hat{r}_{\rm
c}}{\hat{r}}\right)^{4}\,\exp\left[-2\,\hat{v}^{2}_{\rm{crit}}\left(\frac{1}{\hat{r}}-\frac{1}{\hat{r}_{\rm
c}}\right)\right.$ (11)
$\displaystyle\left.-2\left(I_{\hat{g}^{\rm{line}}}(\hat{r})-I_{\hat{g}^{\rm{line}}}(\hat{r}_{\rm
c})\right)-1\right],$

where

$\displaystyle I_{\hat{g}^{\rm{line}}}$ $\displaystyle\equiv$
$\displaystyle\int\hat{g}^{\rm{line}}(\hat{r})d\hat{r}$ (12) $\displaystyle=$
$\displaystyle-\frac{g_{0}\,\hat{r}^{-\delta_{1}}\,{{}_{2}F_{1}}\left[-\gamma,\frac{\delta_{1}}{\delta_{2}},1+\frac{\delta_{1}}{\delta_{2}},\hat{r}^{-\delta_{2}}\right]}{\delta_{1}},$

being ${{}_{2}F_{1}}$ the Gauss hypergeometric function. Note that the
constant of integration vanishes due to the subtraction between the integrals
at $\hat{r}$ and $\hat{r}_{\rm c}$. The critical (or sonic) point,
$\hat{r}_{\rm c}$, is obtained numerically making the RHS of Eq. 9 equal zero.

Finally, taking into account the numerical solution from Hydwind as reference,
a good agreement is obtained with our expression (Eq. 10) for the velocity
profile.

##  4 Line Acceleration Parameters

In Araya et al. (2014) a relationship between the MV08 line-force parameters
($\hat{g_{0}}$, $\delta_{1}$, $\hat{r_{0}}$, and $\gamma$) and the stellar and
m-CAK line-force parameters was given. This relationship is an easy-to-use and
versatile method to compute the velocity profile analytically, because both
stellar and m-CAK line force parameters are already available for a wide range
of spectral types (see, Abbott, 1982; Pauldrach et al., 1986; Lamers &
Cassinelli, 1999; Noebauer & Sim, 2015; Gormaz-Matamala et al., 2019; Lattimer
& Cranmer, 2021).

To derive a similar relationship, now for the $\delta$-slow regime, we created
a grid of m-CAK hydrodynamic models and develop that relationship applying a
multivariate multiple regression (MMR Rencher & Christensen, 2012; Mardia et
al., 1980).

###  4.1 Grid of Hydrodynamic Models

We built a Hydwind grid of stellar models for $\delta$-slow solutions. The
grid points were selected to cover the region of the $T_{\rm{eff}}$–$\log\,g$
diagram where the B- and A-type supergiants are located.

For each given pair of stellar parameters ($T_{\rm{eff}}$, $\log\,g$), the
stellar radius was calculated from $M_{\rm{bol}}$ by means of the flux-
weighted gravity-luminosity relationship (Kudritzki et al., 2003; Kudritzki et
al., 2008), but in addition we added 20 values for stellar radius (from 5
$R_{\sun}$ to 100 $R_{\sun}$ in steps of 5 $R_{\sun}$). The surface gravities
comprise the range of $\log g=2.7$ down to about $90\%$ of the Eddington
limit, in steps of 0.15 dex. We considered 22 effective temperature grid
points, ranging from $9\,000$ K to $19\,500$ K, in steps of 500 K. These
$T_{\rm{eff}}$ and $\log\,g$ ranges were adopted to describe mainly the wind
of intermediate and late B supergiants.

Figure 1: Location of the grid models in the $T_{\rm{eff}}$–$\log\,g$ plane.
Blue and green dots represent the non-converged and converged solutions,
respectively. Red solid lines represent the evolutionary tracks for stars of
$7M_{\sun}$ to $60M_{\sun}$ without rotation (Ekström et al., 2008), while the
black lines correspond to the zero age main-sequence (ZAMS) and the terminal
age main-sequence (TAMS).

The m-CAK line-force parameters used for each set of $(T_{\rm{eff}}$,
$\log\,g$) values are given in Table 1. We considered only high values of
$\delta$ in order to obtain $\delta$-slow solutions.

Then, a huge combinations of parameters were executed in Hydwind, considering
the standard boundary condition at the stellar surface, for the optical depth,
$\tau_{*}=2/3$. In addition, it is worth noting that only some combinations of
all parameters used in Hydwind converged to a physical stationary solution,
i.e., we obtained $141\,067$ $\delta$-slow solutions from our initial set
(about a 2% of our initial input). In the $T_{\rm{eff}}$–$\log\,g$ plane, see
Fig. 1, we show in green dots all converged models, whereas blue dots
indicates that no $\delta$-slow solution was achieved for the given
combination of parameters. Furthermore, the number of converged $\delta$-slow
solutions, in the $T_{\rm{eff}}$–$\log\,g$ plane, shows that most of the
models are concentrated in the region of $\log\,g\geq 1.65$, with a peak
around $T_{\rm{eff}}=14$ kK and $\log\,g=2.4$. Also, few models converged with
values of $\delta\leq 0.28$ and $\alpha\geq 0.57$. This behavior must be
considered at the moment to define the limits of our approximation for
$\delta$-slow solutions.

Finally, for each hydrodynamic model we fitted (Least Squares) the m-CAK line
acceleration ($\hat{g}^{\rm{line}}$) with our proposed line acceleration
expression (Eq. 8) in order to obtain the corresponding new line acceleration
parameters ($\hat{g}_{0}$, $\delta_{1}$, $\delta_{2}$, and $\gamma$).

Table 1: Ranges of the m-CAK line-force parameters for the grid of models.  Parameter | Range
---|---
$\alpha$ | 0.45 – 0.69 (step size of 0.02)
$k$ | 0.05 – 1.00 (step size of 0.05)
$\delta$ | 0.26 – 0.35 (step size of 0.01)

###  4.2 Multivariate Multiple Regression

To derive the relationship for the new line acceleration parameters
($\hat{g}_{0}$, $\delta_{1}$, $\delta_{2}$ and $\gamma$) as function of
stellar ($T_{\mathrm{eff}}$, $\log\mathrm{g}$, $R_{*}/R_{\sun}$) and m-CAK
line-force parameters ($k$, $\alpha$, $\delta$) a MMR is applied to our grid
of models.

A multiple multivariate regression model is:

$Y=XB+Z$ (13)

where $Y$ is a $n\times p$ matrix of data in the $p$ dependent variables, $X$
is a $n\times(1+q)$ matrix of regression: a first column of $1$’s and in the
remaining columns the data of the $q$ independent variables, $B$ is a
$(1+q)\times p$ matrix of parameters (the intercept and $q$ parameters, one
for each of the $q$ independent variables), and $Z$ is a $n\times p$ matrix of
measurement error.

The model is the same for each dependent variable ($y_{i}$, $i$=$1,\dots,p$),
but with different coefficients ($\beta_{ij}$, $i$=$0,\dots,p$;
$j$=$0,\dots,q$), i.e.,

$\displaystyle y_{i}$ $\displaystyle=$
$\displaystyle\beta_{i0}+\beta_{i1}\,T_{\mathrm{eff}}+\beta_{i2}\,\log\mathrm{g}+\beta_{i3}\,R_{*}/R_{\sun}+$
$\displaystyle\beta_{i4}\,k+\beta_{i5}\,\alpha+\beta_{i6}\,\delta+z\hskip
28.45274pt{\rm for}\,\,\,i=1,\dots,p$

where $z$ represents the measurement errors. Each row of $Y$ represents an
observation of each of the $p$ measured response variable. Additional
assumptions in the model are that the expectation of $Y$ is given by $E(Y)=XB$
or $E(Z)=0$, and the covariance matrix of the vectors in the rows of $Y$ is
$\Sigma$, that is, the columns in $Y$ can be correlated. Also, there is an
assumption of normality about the response variables that allows to perform
the hypothesis testing in regression.

For our problem, the dependent variables are $\hat{g}_{0}$, $\delta_{1}$,
$\delta_{2}$ and $\gamma$, and the independent variables are
$T_{\mathrm{eff}}$, $\log\mathrm{g}$, $R_{*}/R_{\sun}$, $k$, $\alpha$ and
$\delta$. The database has $n=141\,067$ records.

A data transformation is necessary to obtain a good fit of the linear model.
Thus, a Box-Cox transformation (Seber & Lee, 2012) is applied to each
dependent variable. This application is performed with the public domain
software R Core Team (2013). The transformations are:
$\hat{g}_{0}\rightarrow\hat{g}_{0}^{0.27}$,
$\delta_{1}\rightarrow(\delta_{1}+1)^{5.3}$,
$\delta_{2}\rightarrow\delta_{2}^{0.45}$,
$\gamma\rightarrow(\gamma+1)^{-3.56}$.

Finally, the estimated parameters are:

$\displaystyle\hat{g}_{0}^{0.27}$ $\displaystyle=$
$\displaystyle-4.548-1.890\times 10^{-4}\,T_{\mathrm{eff}}+$ $\displaystyle
4.393\,\log\mathrm{g}+3.026\times 10^{-2}R_{*}/R_{\sun}-$ $\displaystyle
4.802\times 10^{-3}\,k+3.781\,\alpha-3.212\,\delta,$

$\displaystyle(\delta_{1}+1)^{5.3}$ $\displaystyle=$
$\displaystyle-4.623-3.743\times 10^{-4}\,T_{\mathrm{eff}}+$ $\displaystyle
1.489\times 10^{1}\,\log\mathrm{g}+1.148\times 10^{-1}R_{*}/R_{\sun}+$
$\displaystyle 2.415\,k+9.553\times 10^{1}\,\alpha-1.320\times
10^{2}\,\delta,$

$\displaystyle\delta_{2}^{0.45}$ $\displaystyle=$ $\displaystyle
5.359+8.262\times 10^{-5}\,T_{\mathrm{eff}}-$ $\displaystyle
1.327\,\log\mathrm{g}-8.327\times 10^{-3}R_{*}/R_{\sun}+$ $\displaystyle
2.181\times 10^{-1}\,k+9.618\times 10^{-1}\,\alpha-2.296\,\delta$

and

$\displaystyle(\gamma+1)^{-3.56}$ $\displaystyle=$
$\displaystyle-1.031+7.254\times 10^{-6}\,T_{\mathrm{eff}}+$ $\displaystyle
2.994\times 10^{-1}\,\log\mathrm{g}+3.097\times 10^{-3}R_{*}/R_{\sun}+$
$\displaystyle 1.836\times 10^{-1}\,k-4.828\times
10^{-1}\,\alpha+1.254\,\delta,$

with $R^{2}$ values (proportion of variability of the dependent variable
explained by the regression) given in Table 2. Therefore, the regression
explains almost all the variability of $\hat{g}_{0}^{0.27}$, a large amount of
the variability of $(\delta_{1}+1)^{5.3}$ and $(\gamma+1)^{-3.56}$, and a
minor proportion of $\delta_{2}^{0.45}$.

Table 2: Coefficient of determination ($R^{2}$) of the estimated models.  Model | $R^{2}$
---|---
$\hat{g}_{0}^{0.27}$ | $0.9443$
$(\delta_{1}+1)^{5.3}$ | $0.6016$
$\delta_{2}^{0.45}$ | $0.3408$
$(\gamma+1)^{-3.56}$ | $0.7122$

After fitting the MMR, the estimated values for each dependent variable,
$\hat{g}_{0}^{0.27}$, $(\delta_{1}+1)^{5.3}$, $\delta_{2}^{0.45}$,
$(\gamma+1)^{-3.56}$, are obtained and later transformed into $\hat{g}_{0}$,
$\delta_{1}$, $\delta_{2}$, and $\gamma$ through their respective inverse
functions.

This new relationship for the line acceleration parameters ($\hat{g}_{0}$,
$\delta_{1}$, $\delta_{2}$ and $\gamma$) as function of stellar and m-CAK line
force parameters is valid only for $\delta$-slow solutions, specifically for
values of $\delta$ between 0.29 and 0.35. Therefore, it cannot be compared or
used with others parametrizations obtained using an approximation for the
velocity profile of fast solution (see e.g. Muijres et al., 2012).

##  5 The Approximative Solution

Once we know the relationship (estimated model) between the line acceleration
parameters as a function of the stellar and m-CAK line-force parameters, we
can use Eq. 10 to obtain the velocity profile of the $\delta$-slow wind in
terms of the Lambert W-function.

We point out that considering the number of converged models for some values
of $\alpha$ and $\delta$, we limit our approximation to values of $\alpha$
between $0.45$ and $0.55$, and values of $\delta$ between $0.29$ and $0.35$.
In addition, we could expect a lower precision for values of $\log\,g$ lower
than $1.65$.

In the following of this section, we discuss the accuracy of the terminal
velocities and the derivation of mass-loss rates obtained using this
analytical treatment.

###  5.1 Terminal Velocity

To measure the goodness of fit of the estimated model, the terminal velocity
obtained by Hydwind is compared with our formulated solution.

We consider two terminal velocity vectors: $v_{\infty}^{\mathrm{H}}$ defined
as the terminal velocity calculated with Hydwind (hereafter “true terminal
velocity") and $v_{\infty}^{\mathrm{A}}$ as the terminal velocity obtained
from the our solution at $\hat{r}=r/R_{*}=100$, i.e.,

$v_{\infty}^{\mathrm{A}}=a\,\hat{v}_{\infty}^{\mathrm{A}}=a\,\sqrt{-W_{-1}(x(100))}.$
(19)

The relative error of the estimated terminal velocity
$v_{\infty}^{\mathrm{A}}$ with respect to the true terminal velocity is
calculated by:

${\rm
Relative\,Error}[\%]=100\times\frac{|v_{\infty}^{\mathrm{H}}-v_{\infty}^{\mathrm{A}}|}{v_{\infty}^{\mathrm{H}}},$
(20)

We obtain that the $0.90$ quantile of the distribution of the relative error
are below $21\%$, and the $0.95$ quantile of them are below $27\%$
($q_{0.95}=27.32$).

###  5.2 Mass-loss Rate

Although our solution is developed to obtain a wind velocity profile, we can
derive a recipe to obtain a mass-loss rate. This recipe is based on the m-CAK
theory, specifically the work of Curé (2004), where the velocity profile is
described by our proposed solution (Eq. 10). The full procedure is explained
in Appendix A.

Then, similar to the procedure performed for the terminal velocity, we measure
the goodness of fit of the estimated mass loss rates of the models by
comparing the values (vector) calculated with Hydwind, $\dot{M}^{\mathrm{H}}$,
and the ones obtained with our solution, $\dot{M}^{\mathrm{A}}$.

The relative error of the estimated mass-loss rate $\dot{M}^{\mathrm{A}}$ with
respect to the true mass-loss rate $\dot{M}^{\mathrm{H}}$ was calculated
analogously to the velocity error (Eq. 20). In comparison with the terminal
velocities, the mass-loss rates have slightly higher relative errors. We
observe that most of the data are below $\sim 63\%$. The $0.90$ quantile of
the distribution of the relative error is below $39\%$ and the $0.95$ quantile
is about $46\%$ ($q_{0.95}=46.40$).

Finally, the recipe for the calculation of the estimated mass-loss rate as a
function of stellar and m-CAK line-force parameters ($T_{\mathrm{eff}}$, $\log
g$, $R_{*}/R_{\sun}$, $k$, $\alpha$ and $\delta$) is the following:

  1. 1.

Compute $\hat{g}_{0}$, $\delta_{1}$, $\delta_{2}$ and $\gamma$ from Eqs. 4.2
through 4.2, calculating their respective inverse functions.

  2. 2.

Calculate $v(r)$ from the analytical expression given in Eq. 10.

  3. 3.

Obtain $\dot{M}^{\mathrm{A}}$ using $v(r)$ from (ii) and its gradient in the
m-CAK theory (see Appendix A).

##  6 Discussion and Conclusions

In the frame of the $\delta$-slow wind regime, we have proposed a new
approximate expression for the line force based on the MV08 methodology. This
new expression is a pure function of the radial coordinate and depends on the
following parameters: $\hat{g}_{0}$, $\delta_{1}$, $\delta_{2}$, and $\gamma$.
With this line-force we derived an analytical expressions for the velocity
profile, terminal velocity and a recipe for mass-loss rate (based on m-CAK
theory and our velocity approximation). Furthermore, after generating a grid
of hydrodynamic models, we apply a multivariate multiple regression to obtain
a relationship among these new line-force parameters with the stellar
($T_{\rm{eff}}$, $\log g$, and $R_{*}/R_{\sun}$) and m-CAK line-force
parameters ($\alpha$, $k$, and $\delta$).

The m-CAK line force parameters should be in principle self-consistently
calculated coupling the hydrodynamics with the contribution to the line-
acceleration from hundreds of thousand spectral lines (Lattimer & Cranmer,
2021; Gormaz-Matamala et al., 2019; Pauldrach, 2003, and references therein).
This type of calculations has not been performed for the $\delta$-slow regime,
so far. Nevertheless, based on preliminary line-profile fittings, using the
$\delta$-slow solution, Cidale et al. (2017) found that the value of $\alpha$
is in the same range as in the fast regime, while $k$ is a factor $2$-$3$
lower.

Notwithstanding we can perform a test of our solution using $\alpha$ and $k$
parameters from the fast regime. To this purpose, we consider stellar and wind
parameters from the work of Curé et al. (2011), where they explore the
influence of ionization changes throughout the wind in the velocity profile
for theoretical models of A-type supergiant stars. Thus, we select the models
that match our grid extension (dismissing the region where $\delta\leq 0.28$
and $\alpha\geq 0.57$) in order to compare it to our expression. In addition,
with purpose to test the full range of our work, we also consider the
parameters from models of Venero et al. (2016), where they perform a numerical
study of hydrodynamic solutions within the $\delta$-slow domain, based on
fundamental parameters of typical B supergiants stars.

The stellar and wind parameters from the mentioned works are listed in Table
3. This table also gives the values of the mass-loss rate and terminal
velocity obtained from our analytical solution together with those values
calculated from hydrodynamic results (Hydwind code). All hydrodynamic models
are calculated without stellar rotation.

Table 3: Comparison of the wind parameters obtained via the new analytical solutions ($v_{\infty}^{\mathrm{{\tiny A}}}$, $\dot{M}^{\mathrm{{\tiny A}}}$) with hydrodynamic calculations from Hydwind ($v_{\infty}^{\mathrm{{\tiny H}}}$, $\dot{M}^{\mathrm{{\tiny H}}}$). The models with prefix R and T are from Curé et al. (2011) and Venero et al. (2016), respectively. Model | $T_{\mathrm{eff}}$ | $\log\,g$ | $R_{*}$ | $k$ | $\alpha$ | $\delta$ | $v_{\infty}^{\mathrm{{H}}}$ | $v_{\infty}^{\mathrm{{A}}}$ | $\dot{M}^{\mathrm{{H}}}$ | $\dot{M}^{\mathrm{{A}}}$
---|---|---|---|---|---|---|---|---|---|---
| (kK) | (dex) | ($R_{\odot}$) |  |  |  | (km s-1) | (km s-1) | ($10^{-6}\,M_{\odot}$ yr-1) | ($10^{-6}\,M_{\odot}$ yr-1)
R01 | 11.0 | 2.0 | 70 | 0.37 | 0.49 | 0.29 | 210 | 188 | 0.0052 | 0.0048
R02 | 11.0 | 2.0 | 70 | 0.86 | 0.49 | 0.33 | 201 | 179 | 0.20 | 0.19
R05 | 11.0 | 2.0 | 60 | 0.86 | 0.49 | 0.34 | 185 | 148 | 0.15 | 0.16
R07 | 10.0 | 2.0 | 60 | 0.37 | 0.49 | 0.30 | 207 | 161 | 0.00051 | 0.00042
R08 | 10.0 | 2.0 | 60 | 0.86 | 0.49 | 0.33 | 187 | 155 | 0.017 | 0.017
R11 | 10.0 | 1.7 | 80 | 0.37 | 0.49 | 0.30 | 157 | 116 | 0.0092 | 0.0091
R12 | 10.0 | 1.7 | 80 | 0.86 | 0.49 | 0.34 | 152 | 106 | 0.52 | 0.61
R15 | 9.5 | 2.0 | 60 | 0.37 | 0.49 | 0.30 | 193 | 162 | 0.00015 | 0.00014
R16 | 9.5 | 2.0 | 60 | 0.86 | 0.49 | 0.33 | 136 | 157 | 0.0048 | 0.0047
R19 | 9.5 | 1.7 | 100 | 0.37 | 0.49 | 0.30 | 175 | 185 | 0.0038 | 0.0031
R20 | 9.5 | 1.7 | 100 | 0.86 | 0.49 | 0.34 | 168 | 178 | 0.15 | 0.11
R23 | 9.0 | 1.7 | 100 | 0.37 | 0.49 | 0.33 | 167 | 180 | 0.00025 | 0.00019
R24 | 9.0 | 1.7 | 100 | 0.86 | 0.49 | 0.33 | 171 | 179 | 0.047 | 0.037
T15a | 15.0 | 2.11 | 52 | 0.32 | 0.50 | 0.30 | 200 | 158 | 0.90 | 0.92
T15b | 15.0 | 2.11 | 52 | 0.32 | 0.50 | 0.33 | 191 | 150 | 0.84 | 0.92
T15c | 15.0 | 2.11 | 52 | 0.32 | 0.50 | 0.35 | 186 | 144 | 0.78 | 1.00
T17a | 17.0 | 2.24 | 56 | 0.34 | 0.50 | 0.30 | 236 | 202 | 6.1 | 6.2
T17b | 17.0 | 2.24 | 56 | 0.34 | 0.50 | 0.33 | 225 | 192 | 7.5 | 8.0
T17c | 17.0 | 2.24 | 56 | 0.34 | 0.50 | 0.35 | 220 | 186 | 9.0 | 10.0
T19a | 19.0 | 2.50 | 40 | 0.32 | 0.50 | 0.30 | 270 | 233 | 2.8 | 2.7
T19b | 19.0 | 2.50 | 40 | 0.32 | 0.50 | 0.33 | 257 | 222 | 3.3 | 3.3
T19c | 19.0 | 2.50 | 40 | 0.32 | 0.50 | 0.35 | 251 | 216 | 3.8 | 4.1

The accuracy of our approach is reflected in the low relative errors for the
mass-loss rate and terminal velocity obtained from our solution and the
hydrodynamical code. For the terminal velocity we obtain a relative error mean
and median of $15.6\%$ and $15\%$, respectively. In the case of the mass-loss
rate a relative error mean and median of $10.4\%$ and $7.8\%$ are obtained,
respectively.

The use of approximate expressions that describe closely the hydrodynamics of
stellar winds give the advantage of solving the radiative transfer problem for
moving media in an easy way. In particular, this new expression might properly
describe the winds of late B- and A-type supergiants, without considering a
$\beta$-law with high values ($\beta$ $\gtrsim$ $3$) that lack of any physical
justification in the frame of m-CAK fast solution.

The new expressions for the $\delta$-slow solutions together with the
previously derived expression for the fast solutions (Araya et al., 2014)
provide an easy-to-use procedure to calculate m-CAK wind hydrodynamics.

Furthermore, it is important to remark that these expressions that represent
the hydrodynamics of the wind can be also applied to stellar evolution codes,
where mass loss rates are necessary to estimate the evolutionary phases of a
star.

In future we plan to consider the stellar rotation into our expressions and,
in addition, compare the synthetic line profiles calculated from wind velocity
profiles using a hydrodynamic code and our solutions.

## Acknowledgements

The authors would like to thank the referee, Achim Feldmeier, for his
constructive comments. I.A. thanks the support from FONDECYT project 11190147.
I.A., M.C. & C.A. are grateful with the support from FONDECYT project 1190485.
A.C. & M.C. also acknowledge support from Centro Interdiciplinario de Estudios
Atmosféricos y Astroestadística, Universidad de Valparaíso. C.A., A.G.-M. &
M.C. acknowledge support from Centro de Astrofísica de Valparaíso. C.A. thanks
to FONDECYT project 11190945. R.O.J.V. and L.C. acknowledge financial support
from the Agencia de Promoción Científica y Tecnológica (Préstamo BID PICT
2016/1971), CONICET (PIP 0177), and the Universidad Nacional de La Plata
(Programa de Incentivos 11/G162 and 11/G160, respectively). This project has
received funding from the European Union’s Framework Programme for Research
and Innovation Horizon 2020 (2014-2020) under the Marie Skłodowska-Curie Grant
Agreement No. 823734. We thanks Graeme Candlish and Omar Cuervas for allowing
us to use their computer facilities for the calculation of our grid.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abbott (1982) Abbott D. C., 1982, ApJ, 259, 282
  * Araya et al. (2014) Araya I., Curé M., Cidale L. S., 2014, ApJ, 795, 81
  * Castor et al. (1975) Castor J. I., Abbott D. C., Klein R. I., 1975, ApJ, 195, 157
  * Cidale et al. (2017) Cidale L., Venero R. O. J., Curé M., Haucke M., 2017, in Miroshnichenko A., Zharikov S., Korčáková D., Wolf M., eds, Astronomical Society of the Pacific Conference Series Vol. 508, The B[e] Phenomenon: Forty Years of Studies. p. 11
  * Corless et al. (1993) Corless R. M., Gonnet G. H., Hare D. E. G., Jeffrey D. J., 1993, Maple Technical Newsletter, 9, 12
  * Corless et al. (1996) Corless R. M., Gonnet G. H., Hare D. E. G., Jeffrey D. J., Knuth D. E., 1996, Adv. Comput. Math., 5, 329
  * Cranmer (2004) Cranmer S. R., 2004, American Journal of Physics, 72, 1397
  * Crowther et al. (2006) Crowther P. A., Lennon D. J., Walborn N. R., 2006, A&A, 446, 279
  * Curé (2004) Curé M., 2004, ApJ, 614, 929
  * Curé et al. (2011) Curé M., Cidale L., Granada A., 2011, ApJ, 737, 18
  * Ekström et al. (2008) Ekström S., Meynet G., Maeder A., Barblan F., 2008, A&A, 478, 467
  * Feldmeier (1998) Feldmeier A., 1998, A&A, 332, 245
  * Friend & Abbott (1986) Friend D. B., Abbott D. C., 1986, ApJ, 311, 701
  * Gormaz-Matamala et al. (2019) Gormaz-Matamala A. C., Curé M., Cidale L. S., Venero R. O. J., 2019, ApJ, 873, 131
  * Haucke et al. (2018) Haucke M., Cidale L. S., Venero R., Curé M., Kraus M., Kanaan S., Arcos C., 2018, A&A in press.
  * Kudritzki (2002) Kudritzki R. P., 2002, ApJ, 577, 389
  * Kudritzki et al. (2003) Kudritzki R. P., Bresolin F., Przybilla N., 2003, ApJ, 582, L83
  * Kudritzki et al. (2008) Kudritzki R.-P., Urbaneja M. A., Bresolin F., Przybilla N., Gieren W., Pietrzyński G., 2008, ApJ, 681, 269
  * Lamers & Cassinelli (1999) Lamers H. J. G. L. M., Cassinelli J. P., 1999, in Introduction to Stellar Winds, p. 219
  * Lamers & Rogerson (1978) Lamers H. J. G. L. M., Rogerson Jr. J. B., 1978, A&A, 66, 417
  * Lattimer & Cranmer (2021) Lattimer A. S., Cranmer S. R., 2021, arXiv e-prints, p. arXiv:2101.10375
  * Lefever et al. (2007) Lefever K., Puls J., Aerts C., 2007, A&A, 463, 1093
  * Lucy & Solomon (1970) Lucy L. B., Solomon P. M., 1970, ApJ, 159, 879
  * Mardia et al. (1980) Mardia K. V., Kent J. T., Bibby J. M., 1980, Multivariate Analysis (Probability and Mathematical Statistics). Academic Press
  * Markova & Puls (2008) Markova N., Puls J., 2008, A&A, 478, 823
  * Muijres et al. (2012) Muijres L. E., Vink J. S., de Koter A., Müller P. E., Langer N., 2012, A&A, 537, A37
  * Müller & Vink (2008) Müller P. E., Vink J. S., 2008, A&A, 492, 493
  * Noebauer & Sim (2015) Noebauer U. M., Sim S. A., 2015, MNRAS, 453, 3120
  * Pauldrach (2003) Pauldrach A. W. A., 2003, in Schielicke R. E., ed., Reviews in Modern Astronomy Vol. 16, Reviews in Modern Astronomy. p. 133 (arXiv:astro-ph/0212226), doi:10.1002/9783527617647.ch6
  * Pauldrach et al. (1986) Pauldrach A., Puls J., Kudritzki R. P., 1986, A&A, 164, 86
  * Puls et al. (2000) Puls J., Springmann U., Lennon M., 2000, A&AS, 141, 23
  * R Core Team (2013) R Core Team 2013, R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, http://www.R-project.org/
  * Rencher & Christensen (2012) Rencher A., Christensen W., 2012, Methods of Multivariate Analysis. Wiley Series in Probability and Statistics, Wiley, https://books.google.cl/books?id=0g-PAuKub3QC
  * Searle et al. (2008) Searle S. C., Prinja R. K., Massa D., Ryans R., 2008, A&A, 481, 777
  * Seber & Lee (2012) Seber G., Lee A., 2012, Linear Regression Analysis. Wiley Series in Probability and Statistics, Wiley, https://books.google.cl/books?id=X2Y6OkXl8ysC
  * Stahl et al. (1991) Stahl O., Wolf B., Aab O., Smolinski J., 1991, A&A, 252, 693
  * Venero et al. (2016) Venero R. O. J., Curé M., Cidale L. S., Araya I., 2016, ApJ, 822, 28
  * Verdugo et al. (1999) Verdugo E., Talavera A., Gómez de Castro A. I., 1999, A&AS, 137, 351
  * Villata (1992) Villata M., 1992, A&A, 257, 677
  * Vink et al. (1999) Vink J. S., de Koter A., Lamers H. J. G. L. M., 1999, A&A, 350, 181
  * de Koter et al. (1997) de Koter A., Heap S. R., Hubeny I., 1997, ApJ, 477, 792

##  Appendix A Calculation of the Mass-loss Rate

The calculation of the mass-loss rate $\dot{M}$ is obtained trough the m-CAK
theory, considering the general expression for the line force and the study of
Curé (2004). From this work, we can obtain the location of the singular point
and the mass-loss rate using the singularity and regularity conditions
(expressed with a set of new variables). In our case, the variables related to
velocity are obtained from our proposed solution. It is important to note that
this singular point is the m-CAK one and not the critical point that can be
obtained from Eq. 9 that corresponds to the sonic point.
The change of variables introduced are:

$u=\frac{-R_{*}}{r},\,\,\hat{v}=\frac{v}{a},\,\,{\rm{and}}\,\,\hat{v}^{\prime}=\frac{d\hat{v}}{du}.$
(21)

Considering these new variables, the equation of motion reads:

$\displaystyle F(u,\hat{v},\hat{v}^{\prime})$ $\displaystyle\equiv$
$\displaystyle\left(1-\frac{1}{\hat{v}^{2}}\right)\hat{v}\frac{d\hat{v}}{du}+A+\frac{2}{u}$
(22)
$\displaystyle-C^{\prime}\,FC\,g(u)\,(\hat{v})^{-\delta}\left(\hat{v}\frac{d\hat{v}}{du}\right)^{\alpha}=0,$

where

$A=\frac{G\,M(1-\Gamma)}{a^{2}R_{*}}=\frac{v^{2}_{\mathrm{esc}}}{2a^{2}},$
(23)
$C^{\prime}=C\left(\frac{\dot{M}D}{2\pi}\frac{10^{-11}}{a\,R_{*}^{2}}\right)^{\delta}(a^{2}R_{*})^{(\alpha-1)},$
(24) $C=\Gamma
GMk\left(\frac{4\pi}{\sigma_{E\;}v_{th}\;\dot{M}}\right)^{\alpha}\;\left(\frac{D\dot{M}}{2\pi}\right)^{\delta},$
(25)

and

$g(u)=\left(\frac{u^{2}}{1-\sqrt{1-u^{2}}}\right)^{\delta}.$ (26)

The constant $D$ is defined as:

$D=\frac{(1+Z_{\mathrm{He}}Y_{\mathrm{He}})}{(1+4\,Y_{\mathrm{He}})}\frac{1}{m_{p}},$
(27)

where $m_{\mathrm{p}}$ is the mass of the proton, $Y_{\mathrm{He}}$ is the
helium abundance relative to hydrogen
$\left(n_{\mathrm{He}}/n_{\mathrm{H}}\right)$ and $Z_{\mathrm{He}}$ is the
number of free electrons provided by helium.

To calculate the location of the singular point $u_{c}$, and the eigenvalue,
$C^{\prime}$, it is necessary to satisfy simultaneously, the singularity
condition,

$\frac{\partial}{\partial\hat{v}^{\prime}}F(u,\hat{v},\hat{v}^{\prime})=0\,,$
(28)

and the regularity condition,

$\frac{d}{du}F(u,\hat{v},\hat{v}^{\prime})=\frac{\partial F}{\partial
u}+\frac{\partial F}{\partial\hat{v}}\hat{v}^{\prime}=0.$ (29)

Now, utilizing the change of variables

$Y=\hat{v}\,\hat{v}^{\prime},\,\,{\rm{and}}\,\,\,Z=\frac{\hat{v}}{\hat{v}^{\prime}},$
(30)

Eqs. 22, 28 and 29 are expressed, respectively, as:

$\displaystyle\left(1-\frac{1}{YZ}\right)Y+A+\frac{2}{u}-C^{\prime}\,f_{1}(u,Z)g(u)Z^{-\delta/2}Y^{\alpha-\delta/2}=0,$
(31)
$\displaystyle\left(1-\frac{1}{YZ}\right)Y-C^{\prime}\,f_{2}(u,Z)g(u)Z^{-\delta/2}Y^{\alpha-\delta/2}=0$
(32)
$\displaystyle\left(1+\frac{1}{YZ}\right)Y-\frac{2Z}{u^{2}}-C^{\prime}\,f_{3}(u,Z)g(u)Z^{-\delta/2}Y^{\alpha-\delta/2}=0,$
(33)

See Curé (2004) for the definition of $f_{1}$, $f_{2}$ and $f_{3}$. The set of
equations 31 to 33 are valid for all known solutions from m-CAK theory.

Variables $Y$ and $Z$ are known from our proposed solution (Eq. 10). Now from
Eqs. 32 and 33 we can solve the singular point location, $u=u_{c}$. Note that
$u_{c}\gtrsim 0.1$ to assure a $\delta$-slow solution (Curé et al., 2011).

Finally, the mass loss rate is solved from the eigenvalue, $C^{\prime}$, when
the singular point is replaced in Eq. 22.