# Imaging of the CO snow line in a solar nebula analog

Chunhua Qi1∗, Karin I. Öberg2∗, David J. Wilner1, Paola d’Alessio3,
Edwin Bergin4, Sean M. Andrews1, Geoffrey A. Blake5,
Michiel R. Hogerheijde6, Ewine F. van Dishoeck6,7

1Harvard-Smithsonian Center for Astrophysics
2Departments of Chemistry and Astronomy, University of Virginia
3Centro de Radioastronom a y Astrofisica, Universidad Nacional Autonoma de
Mexico
4Department of Astronomy, University of Michigan
5Division of Geological and Planetary Sciences, California Institute of
Technology
6Leiden Observatory, Leiden University
7Max Planck Institute for Extraterrestrial Physics

∗Contributed equally to this manuscript

> Planets form in the disks around young stars. Their formation efficiency and
> composition are intimately linked to the protoplanetary disk locations of
> “snow lines” of abundant volatiles. We present chemical imaging of the CO
> snow line in the disk around TW Hya, an analog of the solar nebula, using
> high spatial and spectral resolution Atacama Large Millimeter/Submillimeter
> Array (ALMA) observations of N2H+, a reactive ion present in large abundance
> only where CO is frozen out. The N2H+ emission is distributed in a large
> ring, with an inner radius that matches CO snow line model predictions. The
> extracted CO snow line radius of $\sim 30$ AU helps to assess models of the
> formation dynamics of the Solar System, when combined with measurements of
> the bulk composition of planets and comets.

Condensation fronts in protoplanetary disks, where abundant volatiles deplete
out of the gas phase and are incorporated into solids, are believed to have
played a critical role in the formation of planets in the Solar System (?, ?),
and similar “snow lines” in the disks around young stars should affect the
ongoing formation of exoplanets. Snow lines can enhance particle growth and
thus planet formation efficiencies because of 1) substantial increases in
solid mass surface densities exterior to snow line locations, 2) continuous
freeze-out of gas diffusing across the snow line (cold-head effect), 3) pile-
up of dust just inside of the snow line in pressure traps, and 4) an increased
stickiness of icy grains compared to bare ones, which favors dust coagulation
(?, ?, ?, ?, ?). Experiments and theory on these processes have been focused
on the H2O snow line, but the results should be generally applicable to snow
lines of abundant volatiles, with the exception that the “stickiness” of
different icy grain mantles varies. The locations of snow lines of the most
abundant volatiles – H2O, CO2, and CO – with respect to the planet-forming
zone may also regulate the bulk composition of planets (?). Determining snow
line locations is therefore key to probing grain growth, and thus planetesimal
and planet formation efficiencies, and elemental and molecular compositions of
planetesimals and planets forming in protoplanetary disks, including the solar
nebula.

Based on the Solar System composition and disk theory, the H2O snow line
developed at $\sim$3 Astronomical Units (AU, where 1 AU is the distance from
the Earth to the Sun) from the early Sun during the epoch of chondrite
assembly (?). In other protoplanetary disks the snow line locations are
determined by the disk midplane temperature structures, set by a time
dependent combination of the luminosity of the central star, the presence of
other heating sources, the efficiencies of dust and gas cooling, and the
intrinsic condensation temperatures of different volatiles. Because of the low
condensation temperature of CO, the CO snow line occurs at radii of 10’s of AU
around Solar-type stars: this larger size scale makes the CO transition zone
the most accessible to direct observations. The CO snow line is also important
in its own right, because CO ice is a starting point for a complex, prebiotic
chemistry (?). Also without incorporating an enhanced grain growth efficiency
beyond that expected for bare silicate dust, observations of centimeter sized
dust grains in disks, including in TW Hya (?), are difficult to reproduce in
the outer disk. Condensation of CO is very efficient below the CO freeze-out
temperature, with a sticking efficiency close to unity based on experiments
(?), and a CO condensation-based dust growth mode may thus be key to
explaining these observations.

Protoplanetary disks have evolving radial and vertical temperature gradients,
with a warmer surface where CO remains in the gaseous state throughout the
disk, even as it is frozen in the cold, dense region beyond the midplane snow
line (?). This means that the midplane snow line important for planet
formation constitutes a smaller portion of a larger condensation surface.
Because the bulk of the CO emission comes from the disk surface layers, this
presents a challenge for locating the CO midplane snow line. Its location has
been observationally identified toward only one system, the disk around HD
163296, based on (sub-)millimeter interferometric observations of multiple CO
rotational transitions and isotopologues at high spatial resolution,
interpreted through detailed modeling of the disk dust and gas physical
structure (?). An alternative approach to constrain the CO snow line,
suggested in (?) and pursued here, is to image molecular emission from a
species that is abundant only where CO is highly depleted from the gas phase.

N2H+ emission is expected to be a robust tracer of CO depletion because the
presence of gas phase CO both slows down N2H+ formation and speeds up N2H+
destruction. N2H+ forms through reactions between N2 and H${}_{3}^{+}$, but
most H${}_{3}^{+}$ will instead transfer a proton to CO as long as the more
abundant CO remains in the gas phase. The most important destruction mechanism
for N2H+ is proton transfer to a CO molecule, whereas in the absence of CO,
N2H+ is destroyed through a much slower dissociative recombination reaction.
These simple astrochemical considerations predict a correlation between N2H+
and CO depletion, or equivalently an anti-correlation between N2H+ and gas-
phase CO. The latter has been observed in many pre-stellar and protostellar
environments, confirming the basic theory (?, ?). In disks N2H+ should
therefore be present at large abundances only inside the vertical and
horizontal thermal layers where CO vapor is condensing, i.e., beyond the CO
snow line. Molecular line surveys of disks have shown that N2H+ is only
present in disks cold enough to entertain CO freeze-out (?), and marginally
resolved observations hint at a N2H+ emission offset from the stellar position
(?), in agreement with the models of disk chemistry (?). Detailed imaging of
N2H+ emission in protoplanetary disks at the scales needed to directly reveal
CO snow lines with sufficient sensitivity has previously been out of reach.

We used ALMA to obtain images of emission from the 372 GHz dust continuum and
the N2H+ $J=4-3$ line from the disk around TW Hya (Fig. 1, S1) (?). TW Hya is
the closest (54$\pm$6 pc) and as such the most intensively-studied pre-main-
sequence star with a gas-rich circumstellar disk (?, ?). Based on previous
observations of dust and CO emission, and the recent detection of HD line
emission (?), this 3–10 million year old, 0.8 M⊙ T Tauri star (spectral type
K7) is known to be surrounded by an almost face-on ($\sim$6$\degree$
inclination) massive $\sim$0.04 M⊙ gas-rich disk. The disk size in millimeter
dust is $\sim$60 AU, with a more extended ($>$ 100 AU) disk in gas and
micrometer-sized dust (?). Both the disk mass and size conforms well with
solar nebula estimates – the minimum mass of the solar nebula is 0.01 M⊙ based
on planet masses and compositions (?) – and the disk around TW Hya may thus
serve as a template for planet formation in the solar nebula. Our images show
that N2H+ emission is distributed in a ring with an inner diameter of 0.8 to
1.2 arcsec (based on visual inspection), corresponding to a physical inner
radius of 21 to 32 AU. By contrast, CO emission is detected down to radial
scales $\sim$2 AU (?). The clear difference in morphology between the N2H+ and
CO emission can be simply explained by the presence of a CO midplane snow line
at the observed inner edge of the N2H+ emission ring. The different
morphologies cannot be explained by a lack of ions in the inner disk based on
previous spatially and spectrally resolved observations of another ion, HCO+
(?). These HCO+ observations had lower sensitivity and angular resolution than
the N2H+ observations, but they are sufficient to exclude a central hole
comparable in size to that seen in N2H+.

To associate the inner edge radius of the N2H+ emission with a midplane
temperature, and thus a CO freeze-out temperature, requires a model of the
disk density and temperature structure. We adopted the model presented in (?),
updated to conform with recent observations of the accretion rate and grain
settling (Fig. S2–S4, Table S1) (?). In the context of this disk structure
model, the N2H+ inner edge location implies that N2H+ becomes abundant where
the midplane temperature drops to 16–20 K. This is in agreement with
expectations for the CO freeze-out temperature based on the outcome of the
laboratory experiments and desorption modeling by (?), who found CO
condensation/sublimation temperatures of 16–18 K under interstellar
conditions, assuming heat-up rates of 1 K per 102 to 106 years. In outer disk
midplanes, condensation temperatures are expected to at most a few degrees
higher because of a weak dependence on density (?). If CO condenses onto H2O
ice rather than existing CO ice, the condensation temperature will increase
further, but this will only affect the first few monolayers of ice and is not
expected to change the location where the majority of CO freezes out. Some CO
may also remain in the gas phase below the CO freeze-out temperature in the
presence of efficient non-thermal desorption, especially UV photodesorption
(?), but this is expected to be negligible in the disk midplane at 30 AU,
because of UV shielding by upper disk layers. UV photodesorption may affect
the vertical CO snow surface location, however, and it may thus not be
possible to describe the radial and vertical condensation surfaces by a single
freeze-out temperature.

To locate the inner edge of the N2H+ ring more quantitatively, we simulated
the N2H+ emission with a power-law column density distribution and compared
with the data. We assumed the disk material orbits the central star in
Keplerian motion, and fixed the geometric and kinematic parameters of the disk
that affect its observed spatio-kinematic behavior (?). We used the same,
updated density and temperature disk structure model (?), and assumed that the
N2H+ column density structure could be approximated as a radial power-law with
inner and outer edges, while vertically the abundance was taken to be constant
between the lower (toward midplane) and upper (toward surface) boundaries.
This approach crudely mimics the results of detailed astrochemical modeling of
disks, which shows that molecules are predominantly present in well-defined
vertical layers (?, ?), and has been used to constrain molecular abundance
structures in a number of previous studies (?, ?). The inner and outer radii,
power-law index, and column density at 100 AU were treated as free parameters.
We calculated a grid of synthetic N2H+ visibility datasets using the RATRAN
code (?) to determine the radiative transfer and molecular excitation, and
compared with the N2H+ observations. We obtained the best-fit model by
minimizing $\chi^{2}$, the weighted difference between the data and the model
with the real and imaginary part of the complex visibility measured in the
($u,v$)-plane sampled by the ALMA observations of N2H+.

Fig. 2a demonstrates that the inner radius is well constrained to 28–31 AU
(3$\sigma$). This edge determination was aided by the nearly face-on viewing
geometry, because this minimizes the impact of the detailed vertical structure
on the disk modeling outcome. Furthermore, the Keplerian kinematics of the gas
help to constrain the size scale at a level finer than the spatial resolution
implied by the synthesized beam size. As a result, the fitted inner radius is
robust to the details of the density and temperature model (?) (Table S2). In
the context of this model, the best-fit N2H+ inner radius corresponds to a CO
midplane snow line at a temperature of 17 K. Fig. 2b presents the best-fit
N2H+ column density profile together with the best-fit 13CO profile, assuming
a CO freeze-out temperature of 17 K (?) (Fig. S5, Table S3). We fit 13CO
emission (obtained with the Submillimeter Array (?)) because the main
isotopologue CO lines are optically thick. The N2H+ column density contrast
across the CO snow line is at least an order of magnitude (?). Fig. 2c shows
simulated ALMA observations of the best-fit N2H+ $J=4-3$ model, demonstrating
the excellent agreement.

Our quantitative analysis thus confirms the predictions that N2H+ traces the
snow line of the abundant volatile, CO. Furthermore, the agreement between the
quantitative analysis and the visual estimate of the N2H+ inner radius
demonstrate that N2H+ imaging is a powerful tool to determine the CO snow line
radii in disks, whose density and temperature structures have not been modeled
in detail. N2H+ imaging with ALMA may therefore be used to provide statistics
on how snow line locations depend on parameters of interest for planet
formation theory, such as the evolutionary stage of the disks.

The locations of snow lines in solar nebula analogs like TW Hya are also
important to understand the formation dynamics of the Solar System. The H2O
snow line is key to the formation of Jupiter and Saturn (?), while CH4 and CO
freeze-out enhanced the solid surface density further out in the solar nebula,
which may have contributed to the feeding zones of Uranus and Neptune (?),
depending on exactly where these ice giants formed. In the popular Nice model
for the dynamics of the young Solar System, Uranus formed at the largest
radius of all planets, at $\sim 17$ AU (?), and most comets and Kuiper Belt
objects formed further out, to $\sim 35$ AU. The plausibility of this scenario
can be assessed using the bulk compositions of these bodies together with
knowledge of the CO snow line location. In particular, Kuiper Belt objects
contain CO and the even more volatile N2 (?, ?), which implies that they must
have formed beyond the CO snow line. Comets exhibit a range of CO abundances,
some of which seem to be primordial, which suggest the CO snow line was
located in the outer part of their formation region of 15–35 AU (?). This is
consistent with the CO snow line radius that we have determined in the TW Hya
disk. However, in the context of the Nice model, this CO snow line radius is
too large for the ice giants, and suggests that their observed carbon
enrichment has a different origin than the accretion of CO ice (?). A caveat
is that H2O ice can trap CO, though this process is unlikely to be efficient
enough to explain the observations. In either case, the CO snow line locations
in solar nebula analogs like TW Hya offers independent constraints on the
early history of the Solar System.

Fig. 1: Observed images of dust, CO and N2H+ emission toward TW Hya. Left:
ALMA 372 GHz continuum map, extracted from the line free channels of the N2H+
observations. Contours mark [5, 10, 20, 40, 80, 160, 320] mJy beam-1 and the
rms is 0.2 mJy beam-1. Center: image of CO $J=3-2$ emission acquired with the
SMA (?). Contours mark [1, 2, 3, 4, 5] Jy km s-1 beam-1 and the rms is 0.1 Jy
km s-1 beam-1. Right: ALMA image of N2H+ $J=4-3$ integrated emission with a
single contour at 150 mJy km s-1 beam-1 and the rms is 10 mJy km s-1 beam-1.
The synthesized beam sizes are shown in the bottom left corner of each panel.
The red dashed circle marks the best-fit inner radius of the N2H+ ring from a
modeling of the visibilities. This inner edge traces the onset of CO freeze-
out according to astrochemical theory, and thus marks the CO snow line in the
disk midplane.

Fig. 2: Model results for the N2H+ abundance structure toward TW Hya. Upper
left: The $\chi^{2}$ fit surface for the power law index and inner radius of
the N2H+ abundance profile. Contours correspond to the 1–5 $\sigma$ errors and
the blue contour marks 3$\sigma$. Upper right: The best fit N2H+ column
density structure, shown together with the total gas column density and the
best-fit 13CO column density for CO freeze-out at 17 K. The shaded region
marks the N2H+ 1$\sigma$ detection limit. Lower panel: N2H+ observations,
simulated observations of the best-fit N2H+ model, and the imaged residuals,
calculated from the visibilities.

## References and Notes

  * 1. J. S. Lewis, Science 186, 440 (1974).
  * 2. D. J. Stevenson, J. I. Lunine, Icarus 75, 146 (1988).
  * 3. F. J. Ciesla, J. N. Cuzzi, Icarus 181, 178 (2006).
  * 4. A. Johansen, et al., Nature 448, 1022 (2007).
  * 5. E. Chiang, A. N. Youdin, Annual Review of Earth and Planetary Sciences 38, 493 (2010).
  * 6. B. Gundlach, S. Kilias, E. Beitz, J. Blum, Icarus 214, 717 (2011).
  * 7. K. Ros, A. Johansen, A&A 552, A137 (2013).
  * 8. K. I. Öberg, R. Murray-Clay, E. A. Bergin, ApJL 743, L16 (2011).
  * 9. C. Hayashi, Progress of Theoretical Physics Supplement 70, 35 (1981).
  * 10. E. Herbst, E. F. van Dishoeck, ARA&A 47, 427 (2009).
  * 11. D. J. Wilner, P. D’Alessio, N. Calvet, M. J. Claussen, L. Hartmann, ApJL 626, L109 (2005).
  * 12. S. E. Bisschop, H. J. Fraser, K. I. Öberg, E. F. van Dishoeck, S. Schlemmer, A&A 449, 1297 (2006).
  * 13. Y. Aikawa, E. Herbst, A&A 351, 233 (1999).
  * 14. C. Qi, et al., ApJ 740, 84 (2011).
  * 15. C. Qi, K. I. Öberg, D. J. Wilner, ApJ 765, 34 (2013).
  * 16. E. A. Bergin, J. Alves, T. Huard, C. J. Lada, ApJL 570, L101 (2002).
  * 17. J. K. Jørgensen, A&A 424, 589 (2004).
  * 18. K. I. Öberg, et al., ApJ 734, 98 (2011).
  * 19. C. Walsh, H. Nomura, T. J. Millar, Y. Aikawa, ApJ 747, 114 (2012).
  * 20. Materials and methods are available as supporting material on Science Online.
  * 21. J. H. Kastner, B. Zuckerman, D. A. Weintraub, T. Forveille, Science 277, 67 (1997).
  * 22. C. Qi, D. J. Wilner, Y. Aikawa, G. A. Blake, M. R. Hogerheijde, ApJ 681, 1396 (2008).
  * 23. E. A. Bergin, et al., Nature 493, 644 (2013).
  * 24. S. M. Andrews, et al., ApJ 744, 162 (2012).
  * 25. K. A. Rosenfeld, et al., ApJ 757, 129 (2012).
  * 26. D. Hollenbach, M. J. Kaufman, E. A. Bergin, G. J. Melnick, ApJ 690, 1497 (2009).
  * 27. K. Willacy, ApJ 660, 441 (2007).
  * 28. M. R. Hogerheijde, F. F. S. van der Tak, A&A 362, 697 (2000).
  * 29. M. Lecar, M. Podolak, D. Sasselov, E. Chiang, ApJ 640, 1115 (2006).
  * 30. S. E. Dodson-Robinson, K. Willacy, P. Bodenheimer, N. J. Turner, C. A. Beichman, Icarus 200, 672 (2009).
  * 31. K. Tsiganis, R. Gomes, A. Morbidelli, H. F. Levison, Nature 435, 459 (2005).
  * 32. T. C. Owen, et al., Science 261, 745 (1993).
  * 33. S. C. Tegler, et al., ApJ 751, 76 (2012).
  * 34. M. J. Mumma, S. B. Charnley, ARA&A 49, 471 (2011).
  * 35. P. D’Alessio, J. Canto, N. Calvet, S. Lizano, ApJ 500, 411 (1998).
  * 36. P. D’Alessio, N. Calvet, L. Hartmann, S. Lizano, J. Cantó, ApJ 527, 893 (1999).
  * 37. P. D’Alessio, N. Calvet, L. Hartmann, ApJ 553, 321 (2001).
  * 38. P. D’Alessio, N. Calvet, L. Hartmann, R. Franco-Hernández, H. Servín, ApJ 638, 314 (2006).
  * 39. N. I. Shakura, R. A. Sunyaev, MNRAS 175, 613 (1976).
  * 40. N. Calvet, et al., ApJL 630, L185 (2005).
  * 41. C. Espaillat, et al., ApJL 670, L135 (2007).
  * 42. C. Espaillat, et al., ApJ 717, 441 (2010).
  * 43. L. Hartmann, N. Calvet, E. Gullbring, P. D’Alessio, ApJ 495, 385 (1998).
  * 44. A. M. Hughes, D. J. Wilner, C. Qi, M. R. Hogerheijde, ApJ 678, 1119 (2008).
  * 45. N. Calvet, et al., ApJ 568, 1008 (2002).
  * 46. K. I. Uchida, et al., ApJS 154, 439 (2004).
  * 47. A. M. Hughes, D. J. Wilner, S. M. Andrews, C. Qi, M. R. Hogerheijde, ApJ 727, 85 (2011).
  * 48. G. J. Herczeg, B. E. Wood, J. L. Linsky, J. A. Valenti, C. M. Johns-Krull, ApJ 607, 369 (2004).
  * 49. F. L. Schöier, F. F. S. van der Tak, E. F. van Dishoeck, J. H. Black, A&A 432, 369 (2005).
  * 50. P. T. P. Ho, J. M. Moran, K. Y. Lo, ApJL 616, L1 (2004).
  * 51. C. Qi, et al., ApJL 636, L157 (2006).
  * 52. C. Qi, et al., ApJL 616, L11 (2004).

Acknowledgments
We are grateful to S. Schnee for data calibration and reduction assistance.
C.Q. would like to thank the SAO Radio Telescope Data Center (RTDC) staff for
their generous computational support. C.Q., K.I.O. and D.J.W. acknowledges a
grant from NASA Origins of Solar Systems grant No. NNX11AK63. P.D.
acknowledges a grant from PAPIIT-UNAM. E.A.B. acknowledges support from NSF
Grant#1008800. This Report makes use of the following ALMA data:
ADSJAO.ALMA#2011.0.00340.S. ALMA is a partnership of ESO (representing its
member states), NSF (USA), and NINS (Japan), together with NRC (Canada) and
NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint
ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. We also make use of
the Submillimeter Array (SMA) data: project #2004-214 (PI: C. Qi). The SMA is
a joint project between the Smithsonian Astrophysical Observatory and the
Academia Sinica Institute of Astronomy and Astrophysics and is funded by the
Smithsonian Institution and the Academia Sinica.

Supporting Online Material
www.sciencemag.org/cgi/content/full/science.1239560/DC1
Materials and Methods
Table S1–S3
Figs S1–S5
References (35–52)

## Materials and Methods

### Observational details

Continuum and N2H+ line observations toward TW Hya were carried out in ALMA
band 7 (PI: C. Qi) on 19 November, 2012, with 23 to 26 antennas in the Cycle 0
compact configuration. The correlator was configured to observe four windows
with a channel spacing of $\delta\nu$= 61.04 kHz and a bandwidth of 234.375
MHz each. The windows were centered at 372.672 GHz (SPW#1), the rest frequency
of the N2H+ $J=4-3$ line, 372.421 GHz (SPW#0), 358.606 GHz (SPW#2), and
357.892 GHz (SPW#3) . The nearby quasar J1037-295 was used for phase and gain
calibration and 3C279 and J0522-364 were used as bandpass calibrators. The
primary calibrator Ceres provided a mean flux density of 0.61 Jy for the gain
calibrator J1037-295. The visibility data were reduced and calibrated in CASA
3.4.

The atmospheric transmission in the upper (372 GHz) and lower (358 GHz)
sidebands is very different due to strong absorption near 370 GHz. Therefore,
we reduced the data separately for both sidebands. The continuum visibilities
were extracted by averaging the line-free channels in SPW# 0, 1 (upper
sideband) and 2,3 (lower sideband), respectively. We carried out self-
calibration procedures on the continuum as demonstrated in the TW Hya Science
Verification Band 7 CASA Guides, which are available online
(https://almascience.nrao.edu/alma-data/science-verification/tw-hya). The
synthesized beam and RMS for the continuum maps are $0.^{\prime\prime}63\times
0.^{\prime\prime}60$ ($PA=3\degree$), 0.61 mJy beam-1 (upper sideband) and
$0.^{\prime\prime}57\times 0.^{\prime\prime}55$ ($PA=17\degree$), 0.24 mJy
beam-1 (lower sideband). The continuum peak flux densities are determined to
be 2.0097$\pm$0.0062 Jy at 372 GHz and 1.7814$\pm$ 0.0044 Jy at 358 GHz, which
agrees with previous SMA observations and ALMA science verification data (?,
?). We applied the upper sideband continuum self-calibration correction to the
N2H+ $J=4-3$ line data and subtracted the continuum emission in the visibility
domain.

The resulting synthesized beam for the N2H+ $J=4-3$ data cube is
$0.^{\prime\prime}63\times 0.^{\prime\prime}59$ ($PA=-18\degree$), and the
1$\sigma$ rms is 30 mJy beam-1 in 0.1 km s-1 velocity intervals or 8.1 mJy
beam-1 km s-1, which corresponds to a column density detection limit of
2$\times$1011 cm-2 at 17 K. Fig. S1 shows the resulting channel maps for N2H+
$J=4-3$.

### Physical model

The physical model used to interpret the ALMA observations is a steady viscous
accretion disk, heated by irradiation from the central star and mechanical
energy generated by viscous dissipation near the disk midplane (?, ?, ?, ?).
The disk model is axisymmetric, in vertical hydrostatic equilibrium, and the
viscosity follows the $\alpha$ prescription (?). Energy is transported in the
disk by radiation, convection (in regions where the Schwarzschild stability
criterion is not satisfied) and a turbulent energy flux. The penetration of
the stellar and shock generated radiation is calculated using the two first
moments of the radiative transfer equation, taking into account scattering and
absorption by dust grains. This model framework has been used to successfully
reproduce observed disk structures towards several T Tauri and Herbig Ae stars
(?, ?, ?, ?, ?, ?). Following (?, ?, ?), the model includes a tapered
exponential edge to simulate viscous spreading:
$\Sigma\sim\dot{M}\Omega_{k}/\alpha T_{c}$, with
$\alpha(R)=\alpha_{0}exp(-R/R_{c})$. Using the parameter values listed in
Table S1, this physical model provides a good fit to the broadband spectral
energy distribution (SED) of the disk of TW Hya (Fig. S2), except for the mid-
and far-IR wavelengths, which is complicated by contributions of the inner
disk wall, around 3.5–4 AU and optically thin hot dust from the inner hole,
according to (?, ?). The $\alpha_{0}$ adopted in this model is also consistent
with the upper limit on the turbulent line widths, ($<$ 40 m s-1) at $\sim$1–2
scale heights (?), and an ever lower turbulence is expected in the midplane.

Following (?), given a measured disk mass accretion rate, 2$\times$10-9
$M_{\odot}$ yr-1 (?), and the viscous $\alpha$ parameter as formulated above,
the vertical temperature and density structures are mainly regulated by the
degree of grain growth and settling. In the present model, we introduce this
effect in a parametric way. The dust is assumed to consist of two populations
of grains with different size distribution functions (?), with
$a_{max}^{small}=0.25\mu$m (as in the interstellar medium) and
$a_{max}^{big}=1$ mm (consistent with the SED slope at mm wavelengths), and
different spatial distributions such that the abundance of the large grains
increases towards the midplane. The small dust grain to gas mass ratio is
parameterized by $\zeta_{\rm small}$, which is lower than the ISM value
because a large fraction of the dust mass is contained in larger grains. The
amount of dust that is in small dust grains is parameterized by
$\epsilon=\zeta_{\rm small}/\zeta_{\rm ISM}<1$. The amount of large dust
grains (parameterized by $\epsilon_{\rm big}$) is calculated so that the total
dust mass at each radius is conserved, and $\epsilon$ is constrained by the
slope of the SED in the far-IR and sub-mm wavelength range. The maximum grain
size in the model is set to 1 mm, despite observational evidence for the
existence of cm-sized grains (?), because we found a power-law size
distribution of grains with a maximum size in cm could not fit the mm SED
slope. This is probably caused by a radial distribution difference for the mm
and cm-sized dust grains due to differential dust migration. The exclusion of
cm-sized grains should have no effect on the conclusions of this study since
their impact on the midplane temperature structure is minimal.

The surface that separates the regions where the small and large dust grain
populations dominate is parameterized by $z_{\rm big}(R)$ in terms of the
local gas scale heights, $H$. Different $z_{\rm big}$ values are expected to
result in disk structures with different thermal profiles because of changes
in the shape of the irradiation surface, which determines the fraction of
stellar emission intercepted and reprocessed by the disk. We vary $z_{\rm
big}$ between 2$H$ and 3.5$H$, which results in the different vertical disk
temperature structures seen in Fig. S3. Around the observed inner edge of the
N2H+ ring, changing $z_{\rm big}$ also changes the midplane temperature by up
to 2 K (Fig. S4). Despite the importance of $z_{\rm big}$ for vertical
temperature and density structure, Fig. S2 shows that models with different
$z_{\rm big}$ values present small variations in the SEDs, and all models fit
the observed SED satisfactorily. This is generally true for SED modeling
because of the degeneracy of the dust data with the parameter $z_{\rm big}$
(?) and the nearly face-on disk geometry makes SED modeling especially
challenging. The details of the vertical structure in the TW Hya disk are
therefore uncertain, and below we analyze the N2H+ line emission using the
full range of $z_{\rm big}$ values to explore its effect on our conclusions.

### N2H+ line modeling

We adopted the molecular abundance model introduced by (?, ?), and assumed
that the N2H+ emission originates in a vertical layer with a constant
abundance between the surface ($\sigma_{s}$) and midplane ($\sigma_{m}$)
boundaries which are represented by vertically integrated hydrogen column
densities measured from the disk surface in units of 1.59$\times$1021 cm-2. We
fix the vertical surface boundary $\sigma_{s}$ to 3.2 and the midplane
boundary $\sigma_{m}$ to 100, which simulates an emission layer close to the
midplane, in accordance with model predictions. Fitting these boundaries would
require a combination of a well constrained vertical temperature structure and
multiple N2H+ transitions and is thus outside of the scope of this study. To
test the importance of the assumed vertical structure on the inferred radial
distribution of N2H+ we simulate N2H+ visibilities for disk structure models
with the range of $z_{\rm big}$ values and disk structures shown in Figs.
S2-S4.

We model the radial distribution of N2H+ as a power law
N${}_{100}\times(r/100)^{p}$ with an inner radius $R_{in}$ and outer radius
$R_{out}$, where $N_{100}$ is the column density at 100 AU in cm-2, $r$ is the
distance from the star in AU, and $p$ is the power-law index. For each $z_{\rm
big}$ structure model, we compute a grid of synthetic N2H+ $J=4-3$ visibility
datasets over a range of $R_{out}$, $R_{in}$, $p$ and $N_{100}$ values and
compare with the observations. The best-fit model is obtained by minimizing
$\chi^{2}$, the weighted difference between the data and the model with the
real and imaginary part of the complex visibility measured in the
$(u,v)$-plane sampled by the ALMA observations. We use the two-dimensional
Monte Carlo model RATRAN (?) to calculate the radiative transfer and molecular
excitation. The collisional rate coefficients are taken from the Leiden Atomic
and Molecular Database (?).

Table S2 gives the best-fit N2H+ distribution parameters for each $z_{\rm
big}$ model, as well as the corresponding midplane temperature $T_{c}$ at the
N2H+ inner edge. The power law index of the surface density varies between
$2.4$ and $-3.6$, and the column density at the inner edge between
4$\times$1012 and 2$\times$1015 cm-2. Given the 1$\sigma$ detection limit of
2$\times$1011 cm-2, the column density contrast at the inner edge of the N2H+
ring is at least 20 and could be much larger. Across this range of models, the
inner radius varies by less than 5 AU, and the midplane temperature at the
inner radius varies by less than 1 K. The inner radius is thus well
constrained, which implies that the key feature of the N2H+ distribution
needed to constrain the CO snow line is robust with respect to the details of
the physical model assumptions. Channel maps for the best-fit model using the
fiducial $z_{\rm big}=3H$ are shown in Fig. S1 together with the observed data
and imaged residuals, demonstrating the excellent agreement.

### CO distribution

To derive the CO distribution and test the self-consistency of the fiducial
best-fit model, we also modeled two CO isotopologues which we observed with
the Submillimter Array (SMA) (?) in 2005 February 27 and April 10. The main
isotopologue was also observed, but is optically thick and therefore not
included in this analysis. The SMA receivers operated in a double-sideband
mode with an intermediate frequency (IF) band of 4–6 GHz from the local
oscillator frequency, sent over fiber optic transmission lines to 24
overlapping “chunks” of the digital correlator. The correlator was configured
to include CO, 13CO and C18O $J=2-1$ in one setting: the tuning was centered
on the CO $J=2-1$ line at 230.538 GHz in chunk S15, while the 13CO $J=2-1$ at
220.399 GHz and C18O $J=2-1$ at 219.560 GHz were simultaneously observed in
chunk 12 and 22, respectively (?). Combinations of two array configuration
(compact and extended) were used to obtain projected baselines ranging from 6
to 180 m. The observing loops used J1037-295 as the gain calibrator. The
bandpass response was calibrated using observations of 3C279. Flux calibration
was done using observations of Titan and Callisto. Routine calibration tasks
were performed using the MIR software package
(http://www.cfa.harvard.edu/$\sim$cqi/mircook.html), and imaging and
deconvolution were accomplished in MIRIAD. The integrated fluxes are reported
in Table S3. Fig. S5 shows the spatially integrated spectra of 13CO and C18O
$J=2-1$ extracted from the SMA channel maps in 8′′ square boxes centered on TW
Hya.

Following (?), CO is assumed to be present in the disk between a lower
boundary set by the CO freeze-out temperature derived from the N2H+ modeling,
and an upper boundary set by photodissociation, though the choice of upper
boundary in this case has a very small effect on the modeled emission profiles
of 13CO and C18O. The CO abundance structure was optimized using the same
procedure as for N2H+ above. Fig. S5 shows that the best-fit CO abundance
distribution fits the CO isotopologue observations well when using the
fiducial disk structure and assuming standard isotope ratios and CO freeze-out
at the N2H+ inner edge temperature of 17 K. In contrast the models with much
smaller $z_{\rm big}$ cannot reproduce the relative C18O and 13CO fluxes
without order of magnitude deviations from the cosmic isotope ratios. More
data with better sensitivity and resolution from the emission of CO and its
isotopologues and a detailed surface heating model (as suggested by (?)) are
needed to constrain the temperature structure and the $z_{\rm big}$ value in
the disk of TW Hya.

Table S1: Physical model for the disk of TW Hya  Parameters | Values
---|---
Stellar and accretion properties
Spectral type | K7
Effective temperature: $T_{*}$ (K) | 4110
Estimated distance: $d$ (pc) | 54
Stellar radius: $R_{*}$ (R⊙) | 1.04
Stellar mass: $M_{*}$ (M⊙) | 0.8
Accretion rate: $\dot{M}$ ($M_{\odot}$ yr-1) | 2$\times$10-9
Disk structure properties
Disk mass: $M_{d}$ (M⊙) | 0.04
Characteristic radius: $R_{c}$ (AU) | 60
Viscosity coefficient: $\alpha_{0}$ | 0.0007
Depletion factor of the atmospheric small grains: $\epsilon^{a}$ | 0.01
$z_{\rm big}$a (Hb) | 3.0
Disk geometric and kinematic propertiesc
Inclination: $i$ (deg) | 6
Systemic velocity: $V_{\rm LSR}$ (km s-1) | 2.86
Turbulent line width: $\delta$vturb (km s-1) | 0.05
Position angle: $P.A.$ (deg) | 155
aSee definition in paper. |
bGas scale height. |
cParameters adopted from (?, ?, ?, ?). |
Table S2: N2H+ $J=4-3$ fitting resultsa $z_{\rm big}$ ($H$) |  $R_{in}$ (AU) |  $T_{\rm c}$b (K) | $p$ |  $N_{100}$ (cm-2) |  $R_{out}$ (AU)
---|---|---|---|---|---
2.0 | 25${}^{+4}_{-6}$ | 15–18 | 2.4${}^{+0.6}_{-0.3}$ | (1.4$\pm$0.2) $\times$1014 | 150$\pm$10
2.5 | 30${}^{+1}_{-3}$ | 15–16 | 0.4${}^{+0.6}_{-0.4}$ | (2.9$\pm$0.5) $\times$1014 | 140$\pm$10
3.0 | 30${}^{+1}_{-2}$ | 16–17 |  $-$2.0${}^{+0.5}_{-0.7}$ | (1.6$\pm$0.3) $\times$1014 | 140$\pm$10
3.5 | 30${}^{+1}_{-4}$ | 16–18 |  $-$3.6${}^{+0.6}_{-0.8}$ | (2.5$\pm$0.4) $\times$1013 | 140$\pm$10
aErrors within 3$\sigma$. |  |  |  |  |
bTemperature range based on Fig. S4. |  |  |  |  |
Table S3: TW Hya CO isotopologue observation results. Transition | Frequency (GHz) | Beam /$P.A.$ |  $\int Fdv$ (Jy km s-1)
---|---|---|---
13CO $J=2-1$ | 220.399 |  $2.^{\prime\prime}7\times 1.^{\prime\prime}8$ / $-$3.0° | 2.72[0.18]
C18O $J=2-1$ | 219.560 |  $2.^{\prime\prime}8\times 1.^{\prime\prime}9$ / $-$1.3° | 0.68[0.18]

Fig. S1:  Channel maps of the N2H+ $J=4-3$ line emission observed by ALMA from
the disk around TW Hya. The LSR velocity is indicated in the upper right of
each channel, while the synthesized beam size and orientation
($0.^{\prime\prime}63\times 0.^{\prime\prime}59$ at a position angle of
$-$18.1∘) is indicated in the lower left panel. The contours are 0.03
(1$\sigma$) $\times[3,6,9,12,15,18]$ Jy beam-1 . Fig. S2: The TW Hya SED and
the model results for $\dot{M}=2\times 10^{-9}M_{\odot}/yr$,
$\alpha_{0}=0.0007$, $\epsilon=0.01$, and $R_{c}=60$ AU. See the SED
references in (?). The new mm/submm fluxes are from the ALMA science
verification data and this paper (marked by green triangles). The different
model SED lines correspond to $z_{\rm big}/H=2-3.5$, with the fiducial 3H
model shown with a solid line, demonstrating that the SED modeling does not
provide strong constraints on this parameter.  Fig. S3: Vertical temperature
and density profiles at $R=30$ AU for disk models with $\dot{M}=2\times
10^{-9}\ M_{\odot}/yr$, $\alpha_{0}=0.0007$, $\epsilon=0.01$, $R_{c}=60$ AU,
and different values of $z_{\rm big}$ values. Upper panel: Temperature versus
height at $R=30$ AU for disk models with different values of $z_{\rm big}/H=$
2.0, 2.5, 3.0, 3.5 (from left to right). The fiducial model, with $z_{\rm
big}=3H$ is shown with a solid line. Lower panel: Density versus height at
$R=30$ AU for the same models. The lines from top to bottom correspond to
models with $z_{\rm big}/H=$ 2.0, 2.5, 3.0, 3.5.  Fig. S4: Midplane
temperature profiles for the TW Hya disk for different $z_{\rm big}$ values
showing the effect of $z_{\rm big}$ on the midplane temperature around the
N2H+ inner edge.  Fig. S5: CO isotopologue lines toward TW Hya, observed with
the SMA (grey) and the best-fit CO modeling results using the fiducial disk
structure developed to interpret the N2H+ observations. The dashed line marks
the $V_{\rm LSR}$ toward TW Hya.


# Tunneling time from spin fluctuations in Larmor clock

Durmuş Demir

###### Abstract

Tunneling time, time needed for a quantum particle to tunnel through a
potential energy barrier, can be measured by a duration marker. One such
marker is spin reorientation due to Larmor precession. With a weak magnetic
field in $z$ direction, the Larmor clock reads two times, $\tau_{y}$ and
$\tau_{z}$, for a potential energy barrier along the $y$ axis. The problem is
to determine the actual tunneling time (ATT). Büttiker defines
$\sqrt{\tau_{y}^{2}+\tau_{z}^{2}}$ to be the ATT. Steinberg and others, on the
other hand, identify $\tau_{y}$ with the ATT. The Büttiker and Steinberg times
are based on average spin components but in non-commuting spin system average
of one component requires the other two to fluctuate. In the present work, we
study the effects of spin fluctuations and show that the ATT can well be
$\tau_{y}+\frac{\tau_{z}^{2}}{\tau_{y}}$. We analyze the ATT candidates and
reveal that the fluctuation-based ATT acts as a transmission time in all of
the low-barrier, high-barrier, thick-barrier and classical dynamics limits. We
extract this new ATT using the most recent experimental data by the Steinberg
group. The new ATT qualifies as a viable tunneling time formula.

###### keywords:

Larmor time , Spin uncertainty , Tunneling time

††journal: Physics Letters A

[inst1] organization=Sabancı University, addressline=Faculty of Engineering
and Natural Sciences, city=Tuzla, postcode=34956, state=Istanbul,
country=Turkey

##  1 Introduction

Tunneling is the transmission of quantum particles through potential energy
barriers which exceed their total energy [1, 2]. It is a pure quantum effect
underlying various natural phenomena and technological devices. Its duration,
the tunneling time, has been formulated by utilizing various observables with
various methods [3, 4, 5, 6]. The reason for the existence of various
tunneling time formulas is that the role of time in quantum theory is a
controversial issue. There is the notion of observable time (flight time)
hampered by Pauli’s famous theorem [7, 8, 9]. There is also the notion of
external time related to surroundings of the quantum system [9, 10]. There is
again the notion of intrinsic time that can be marked by an observable
intrinsic to the quantum system. It is this intrinsic time which will be
studied in the present work. By definition, march of the intrinsic time can be
parametrized by using an appropriate marker namely an observable
representative of the quantum system. The marker could be phase shift [12],
spin precession [13, 14], entropy production [15] or any other change
pertinent to the tunneling particle. The spin precession across feebly-
magnetized potential barriers, for instance, meters the tunneling time without
altering energetics of the particle. This marker, the Larmor clock [16, 17,
18], reads two times, $\tau_{y}$ and $\tau_{z}$, with a feeble magnetic field
in the $z$ direction and a potential energy barrier along the $y$ axis. These
times have recently been measured by Steinberg group using the rubidium atoms
[19]. There is, in general, no telling what combination of $\tau_{y}$ and
$\tau_{z}$ forms the actual tunneling time (ATT) [16, 20, 21, 8]. Büttiker
defines $(\tau_{y}^{2}+\tau_{z}^{2})^{1/2}$ to be the ATT. Steinberg and
others, on the other hand, identify $\tau_{y}$ with the ATT after attributing
$\tau_{z}$ to measurement back-action [22, 20, 19]. The key point about
$\tau_{y}$ and $\tau_{z}$ is that they are based on average spin components,
and do not therefore involve possible contributions from spin fluctuations.

The Larmor precession times $\tau_{y}$ and $\tau_{z}$ are set by the average
spin components. They do not carry information on how the transmitted spins
are distributed about their averages. Their distributions are important
because different spin orientations imply different precession times. In this
sense, variances of the spins normalized to their average values (namely the
Fano factors [23, 24] for spin components) act as a measure of how dispersed
the spin components get with respect to their average values in course of the
tunneling. Fano factor for a given spin component is equal to $\hbar/2$ if
spins are Poisson-distributed, larger if spins are clustered, and smaller if
spins are uniformly-distributed [25]. Needless to say, each spin distribution
leads to a precession time distribution of its own, and, in this sense, the
Fano factor is expected to lead to a factual determination of the ATT.

In the present work, we study the effects of spin fluctuations on Larmor
clock. We find that spin uncertainty, whose relevance is hinted by the spin
component in the direction of the magnetic field, gives rise to a new
tunneling time formula. This fluctuation-induced ATT applies directly to
experimental data [19], and differs from the Büttiker and Steinberg ATT
definitions in terms of the transmission speeds involved. In Sec. 2 below we
derive the ATT from spin fluctuations, and contrast it with the Büttiker and
Steinberg definitions in asymptotic limits. In Sec. 3 we analyze the
fluctuation-induced ATT we constructed using the available experimental data
[19], and discuss its physics implications considering the relevant limit
values. In Sec. 4 we conclude with a discussion of the fluctuation-induced ATT
in regard to its potential applications like the quantum annealing [26, 27].

##  2 Actual Tunneling Time from Spin Fluctuations

Let us consider an ensemble of identical particles each having mass $m$,
magnetic moment $\mu$ and gyromagnetic ratio $g$. The ensemble is endowed with
average energy $E$ and average spin components

$\displaystyle\langle S_{x}\rangle=\hbar/2\,,\;\langle
S_{y}\rangle=0\,,\;\langle S_{z}\rangle=0$ (1)

as sharply-peaked distributions with small variance-to-mean ratios. This
ensemble moves along the $y$-axis with average momentum $\sqrt{2mE}$ to
scatter at a potential energy profile $V(y)$, which confines a uniform
magnetic field $B$ in the $z$-direction. The particles get eventually either
reflected or transmitted. Under the effect of the magnetic field, the
transmitted particles attain the average spin components [16, 17, 18]

$\displaystyle\langle S_{x}\rangle=\hbar/2\,,\;\langle
S_{y}\rangle=(\hbar/2)\omega_{L}\tau_{y}\,,\;\langle
S_{z}\rangle=(\hbar/2)\omega_{L}\tau_{z}$ (2)

up to ${\cal{O}}(\omega_{L}^{2})$ accuracy in the Larmor frequency
$\hbar\omega_{L}=g\mu B$ in a feeble magnetic field $(\hbar\omega_{L}\ll E)$.
This transmitted spin differs from the incident spin in (1) by its nonzero
$y$\- and $z$-components, and gives rise therefore to two distinct precession
times, $\tau_{y}$ and $\tau_{z}$. The time $\tau_{y}$ is a direct consequence
of the magnetic field because it remains nonzero even when the potential is
absent. The time $\tau_{z}$, on the other hand, results from the potential in
that it vanishes in the limit of vanishing potential. There is, nevertheless,
no telling what combination of $\tau_{y}$ and $\tau_{z}$ is the ATT. In fact,
in his seminal work [16], Büttiker defines the ATT to be

$\displaystyle({\rm ATT})_{B}=\sqrt{\tau_{y}^{2}+\tau_{z}^{2}}$ (3)

as the quadrature sum of the $\tau_{y}$ and $\tau_{z}$. In the weak
measurement formalism [22], on the other hand, Steinberg takes the ATT to be

$\displaystyle({\rm ATT})_{S}=\tau_{y}$ (4)

by attributing $\tau_{z}$ to measurement back-action [20, 21]. There are thus
two proposed formulas for the duration of tunneling through a given barrier.
Below, we construct a third formula based on spin fluctuations.

The definitions of $\tau_{y}$ and $\tau_{z}$ in (2) make it clear that both
$({\rm ATT})_{B}$ and $({\rm ATT})_{S}$ are based on the average spin
components $\langle S_{y,z}\rangle$. It is, however, not possible to measure
any one of the spin components without disrupting the other two. This follows
from non-commuting nature of the spin components. It is in this sense that the
spin uncertainty comes to the fore, and puts $\langle S_{z}\rangle$ at the
center stage in agreement with its role as the measurement back-action in the
weak measurement formalism [20, 21]. In fact, $\langle S_{z}\rangle$ is a
sensitive probe of the potential and, at the same time, underpins the
uncertainly relation

$\displaystyle(\Delta S_{x})^{2}(\Delta
S_{y})^{2}\geq\frac{\hbar^{2}}{4}\langle S_{z}\rangle^{2}$ (5)

according to which the $S_{x}$ and $S_{y}$ variances ($i=x,y$)

$\displaystyle(\Delta S_{i})^{2}=\langle(S_{i}-\langle
S_{i}\rangle)^{2}\rangle=\langle S_{i}^{2}\rangle-\langle S_{i}\rangle^{2}$
(6)

are related to each other seesawically, with a pivot at $\langle
S_{z}\rangle$. In view of these variances, the ratio $(\Delta
S_{x,y})^{2}/\langle S_{x,y}\rangle$, which is known as the Fano factor [23,
24], acts a measure of how dispersed the spin components get with respect to
their average values in course of the transmission (tunneling). It is equal to
$\hbar/2$ if spins are Poisson-distributed, larger if spins are clustered, and
smaller if spins are uniformly-distributed [25].

The uncertainty product $(\Delta S_{x})^{2}(\Delta S_{y})^{2}$ in (5) is
expected to have a strong correlation with the potential barrier. This follows
from sensitivity of $\langle S_{z}\rangle$ to the potential profile, and comes
to mean that the time delay caused by the potential, the sought-for ATT
itself, must actually be encoded in the product $(\Delta S_{x})^{2}(\Delta
S_{y})^{2}$ via $\langle S_{z}\rangle$. In this regard, we first reorganize
$(\Delta S_{x})^{2}(\Delta S_{y})^{2}$ as normalized Fano factors. Next, we
observe that the sought-for ATT (transmission time), irrespective of how it is
formulated, must be measured in terms of the inverse Larmor frequency
$\omega_{L}^{-1}$ because it is the only time scale in the problem. In view of
these features, spin fluctuations can be conjectured to give rise to a new ATT

$\displaystyle({\rm ATT})_{F}=\omega_{L}^{-1}\times\frac{(\Delta
S_{x})^{2}}{\frac{\hbar}{2}\langle S_{x}\rangle}\times\frac{(\Delta
S_{y})^{2}}{\frac{\hbar}{2}\langle S_{y}\rangle}$ (7)

in which the subscript “F” is a reminder that this new ATT is based on spin
“fluctuations”. Up to ${\mathcal{O}}(\omega_{L}^{2})$ accuracy under a feeble
magnetic field $(\hbar\omega_{L}\ll E)$, the spin variances in (7) take the
forms [16]

$\displaystyle(\Delta S_{x})^{2}$ $\displaystyle=$
$\displaystyle\frac{\hbar^{2}}{4}\omega_{L}^{2}(\tau_{y}^{2}+\tau_{z}^{2})$
(8) $\displaystyle(\Delta S_{y})^{2}$ $\displaystyle=$
$\displaystyle\frac{\hbar^{2}}{4}\left(1-\omega_{L}^{2}\tau_{y}^{2}\right)$
(9)

after using the spin averages in (2) in the formula (6) for spin variances by
taking into account the fact that $\langle S_{i}^{2}\rangle={\hbar^{2}}/{4}$
(for each $i=x,y,z$) by the property of the Pauli spin matrices. Now, using
these spin variances in (7) along with the spin averages in (2) one is led to
the fluctuation-induced transmission time formula

$\displaystyle({\rm ATT})_{F}=\tau_{y}+\frac{\tau_{z}^{2}}{\tau_{y}}$ (10)

under a feeble magnetic field $(\hbar\omega_{L}\ll E)$ for which
${\cal{O}}(\omega_{L}^{2})$ terms in (7) can all be dropped, and average spin
components in (2) can be measured without disrupting the state of the particle
[22, 20, 21]. (It is necessary to suppress ${\mathcal{O}}(\omega_{L}^{2})$
contributions as they enhance also orbital effects.) It differs from the
earlier transmission time definitions [16, 22, 20] not only by its functional
form but also by its base rock of spin uncertainty. It has general validity in
that it is independent of what shape [17, 18] the potential has, how wide or
high the tunneling region $(V(y)>E)$ is, and if the potential is time-
dependent [28] or not [17, 18]. It can be extracted directly from $\tau_{y}$
and $\tau_{z}$ data irrespective of if these two precession times are measured
by experiment or modeled by theory.

Table 1: The three ATT candidates in the low-barrier, high-barrier, thick-barrier and classical dynamics limits. | $\tau_{y}$ | $\tau_{z}$ | $({\rm ATT})_{B}$ | $({\rm ATT})_{S}$ | $({\rm ATT})_{F}$
---|---|---|---|---|---
| low-barrier: $V_{0}\ll E$
---
(fixed $E$)
$\tau_{c}(0,E)$ | $0$ | $\tau_{c}(0,E)$ | $\tau_{c}(0,E)$ | $\tau_{c}(0,E)$
| high-barrier: $E\ll V_{0}$
---
(fixed $V_{0}$)
$0$ | $\tau_{c}(V_{0},0)$ | $\tau_{c}(V_{0},0)$ | $0$ | $\infty$
| thick-barrier: $L^{2}\gg\frac{\hbar}{m}\tau_{c}(V_{0},E)$
---
(fixed $V_{0},E$)
$\frac{\hbar}{V_{0}}\frac{\tau_{c}(V_{0},E)}{\tau_{c}(0,E)}$ | $\infty$ | $\infty$ | $\frac{\hbar}{V_{0}}\frac{\tau_{c}(V_{0},E)}{\tau_{c}(0,E)}$ | $\infty$
| classical dynamics: $\hbar\rightarrow 0$
---
(fixed $V_{0},E,L$)
$0$ | $\tau_{c}(V_{0},E)$ | $\tau_{c}(V_{0},E)$ | $0$ | $\infty$

There are three ATT candidates $({\rm ATT})_{B}$, $({\rm ATT})_{S}$ and $({\rm
ATT})_{F}$. The question of which ATT is realized in nature will eventually be
answered by experiment. In the absence of experimental data, all one can do is
to contrast these ATTs in physically discriminate configurations. To this end,
a rectangular potential barrier of height $V_{0}$ and width $L$ proves
particularly eligible. Indeed, this configuration admits a complete analytic
illustration and, as a result, one finds that [16, 17]

$\displaystyle\tau_{y}=\frac{\hbar}{2}\left(\frac{E}{V_{0}-E}\right)^{1/2}\frac{2\lambda(V_{0}-2E)+V_{0}\sinh
2\lambda}{4E(V_{0}-E)+V_{0}^{2}\sinh^{2}\lambda}$ (11)

and

$\displaystyle\tau_{z}=\frac{\hbar}{4}\frac{V_{0}}{V_{0}-E}\frac{2(V_{0}-2E)\sinh^{2}\lambda+V_{0}\lambda\sinh
2\lambda}{4E(V_{0}-E)+V_{0}^{2}\sinh^{2}\lambda}$ (12)

with the barrier opacity $\hbar\lambda=L(2m(V_{0}-E))^{1/2}$. The ATT
candidates $({\rm ATT})_{B}$, $({\rm ATT})_{S}$ and $({\rm ATT})_{F}$ derive
from $\tau_{y}$ and $\tau_{z}$ via their definitions in (3), (4) and (10).
Their physics implications are best revealed by contrasting them in their
asymptotic forms. Indeed, their forms in the low-barrier, high-barrier, thick-
barrier and classical dynamics limits can shed light on their capabilities as
transmission times. In this regard, they are compared in Table 1 for each of
these limits. In the low-barrier limit ($E/V_{0}\gg 1$), they all agree and
give a unique ATT

$\displaystyle\tau_{c}(V_{0},E)=\frac{mL}{\sqrt{2m|V_{0}-E|}}$ (13)

designating the time it takes for a classical particle to traverse a distance
$L$ with effective energy $|V_{0}-E|$ [16, 20, 29]. This agreement of the
three ATT candidates, as given in Table 1, stems from vanishing of $\tau_{z}$
for low-potentials. It is in this sense that $\tau_{z}$ is a sensitive probe
of the potential. It qualifies to be a viable measure of the march of time in
the presence of the potential barrier.

The potential barrier becomes opaque ($\lambda\gg 1$) in the limit of thick
and high barriers. It becomes opaque in the limit also of the a classical
dynamics. In this opaque limit, the three ATT candidates take the asymptotic
forms

$\displaystyle\lim_{\lambda\rightarrow\infty}({\rm ATT})_{B}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\hbar^{2}}{V^{2}_{0}}\frac{\tau^{2}_{c}(V_{0},E)}{\tau^{2}_{c}(0,E)}+\tau^{2}_{c}(V_{0},E)},$
$\displaystyle\lim_{\lambda\rightarrow\infty}({\rm ATT})_{S}$ $\displaystyle=$
$\displaystyle\frac{\hbar}{V_{0}}\frac{\tau_{c}(V_{0},E)}{\tau_{c}(0,E)}\,,$
(14) $\displaystyle\lim_{\lambda\rightarrow\infty}({\rm ATT})_{F}$
$\displaystyle=$
$\displaystyle\frac{\hbar}{V_{0}}\frac{\tau_{c}(V_{0},E)}{\tau_{c}(0,E)}+\frac{V_{0}\tau_{c}(V_{0},E)\tau_{c}(0,E)}{\hbar},$

in which $\tau_{c}(0,E)\equiv\tau_{c}(V_{0}=0,E)$ where $\tau_{c}(V_{0},E)$ is
defined in (13). These high-opacity limits encompass each of the classical,
thick-barrier and high-barrier limits, as elaborated below:

  1. 1.

In the high-barrier limit ($E/V_{0}\ll 1$), $({\rm ATT})_{B}$ and $({\rm
ATT})_{S}$, as given in Table 1, exhibit non-quantum behavior because they
both remain finite for tunneling through the barrier. This is not expected in
quantum mechanics since wavefunction is diminished in the barrier region, and
this diminishing suppresses the transmission exponentially to give cause to a
long transmission (tunneling) time [16, 17, 18]. In contrast to $({\rm
ATT})_{B}$ and $({\rm ATT})_{S}$, the fluctuation-based tunneling time $({\rm
ATT})_{F}$ becomes infinitely long and this is in perfect agreement with the
quantum expectations.

  2. 2.

In the thick-barrier limit ($L^{2}\gg\frac{\hbar}{m}\tau_{c}(V_{0},E)$), as
given in Table 1, $({\rm ATT})_{S}$ becomes width-independent, and $({\rm
ATT})_{B}$ and $({\rm ATT})_{F}$ both tend to infinity. It is clear that
$({\rm ATT})_{S}$ exhibits again non-quantum behavior. The longevity of $({\rm
ATT})_{B}$ and $({\rm ATT})_{F}$, on the other hand, is precisely what is
expected of transmission (tunneling) time in view again of the suppression of
the transmission probability under the barrier [16, 17, 18]. It is possible to
infer that $({\rm ATT})_{S}$ can pertain to group velocity since group
velocity can be superluminal. The other two, $({\rm ATT})_{B}$ and $({\rm
ATT})_{F}$, on the other hand, could result from wave-front velocity (the
subluminal speed with which waves carry information) [30, 31, 32].

  3. 3.

In the classical particle limit ($\hbar\rightarrow 0$), as given in Table 1,
$({\rm ATT})_{B}$ and $({\rm ATT})_{S}$ exhibit non-classical behavior as they
both remain finite. This is not expected in classical mechanics in which
tunneling through barriers is simply impossible. The fluctuation-induced
tunneling time $({\rm ATT})_{F}$ becomes infinitely long. This longevity is
precisely what is expected of transmission (tunneling) time for a classical
particle [1, 4, 8].

These asymptotic behaviors lead to an unambiguous conclusion which is that the
$({\rm ATT})_{F}$ possesses all the features expected of a transmission
(tunneling) time. On physical grounds, therefore, $({\rm ATT})_{F}$ is
expected to be a serious ATT candidate, possibly pertaining to wave-front
speed [30, 31, 32].

##  3 Extracting Actual Tunneling Time from Experimental Data

In this section, we shall extract $({\rm ATT})_{F}$ from the available
experimental data, and discuss its asymptotics. Experiments with inert gases
have already given strong evidence for finite tunneling time [33, 34, 35].
These atomic evidences are hardly eligible for model building since it is not
possible to know the start time of the electron tunneling [15]. The
experiments with cold atoms have turned these evidences into accurate
measurements. Indeed, $\tau_{y}$ and $\tau_{z}$ have recently been measured by
Steinberg group [19] on an ensemble of some 8,000 Bose-condensed 87Rb
isotopes. In the experiment, the ensemble is staged [36, 37, 38] to scatter at
an optical Gaussian barrier of peak value $V_{0}=135\ {\rm nK}$ and width
$L_{0}=1.3\ {\rm\mu m}$, and spin components of the transmitted isotopes are
measured via their absorption spectra. The experimental data agrees with basic
theoretical predictions (Table 1) in that, for low potentials
$(E/V_{0}\rightarrow\infty)$, $\tau_{z}$ vanishes $(\tau_{z}\rightarrow 0)$
but $\tau_{y}$ remains finite $(\tau_{y}\rightarrow mL_{0}/\sqrt{2mE})$. The
data does not extend to slower isotopes but still it exhibits an indistinct
trend that, for high potentials $(E/V_{0}\rightarrow 0)$, $\tau_{y}$ vanishes
$(\tau_{y}\rightarrow 0)$ but $\tau_{z}$ remains finite. The actual
transmission duration of the 87Rb atoms, modeled by the $({\rm ATT})_{F}$ in
(10) as a fluctuation-induced tunneling time, lies above both $\tau_{y}$ and
$\tau_{z}$ on the average, and has error bars of similar size as theirs (Fig.
1). The $({\rm ATT})_{F}$ approaches to $\tau_{y}$ at low potentials
$(E/V_{0}\rightarrow\infty)$ and shows a slight trend that it may blow up at
high potentials $(E/V_{0}\rightarrow 0)$.

Figure 1: The $({\rm ATT})_{F}$ (green/triangle) as derived via (10) from the
experimental data [19] on precession times $\tau_{y}$ (red/bullet) and
$\tau_{z}$ (blue/diamond) of the 87Rb isotopes transmitted with average energy
$E$ through a Gaussian barrier of peak value $V_{0}=135\ {\rm nK}$ and width
$L_{0}=1.3\ {\rm\mu m}$. There is an indistinct trend that, in the deep
tunneling regime $(E\ll V_{0})$, $\tau_{y}$ $(\tau_{z})$ can decrease
(increase) and, as a result, the $({\rm ATT})_{F}$ can take large values.
Figure 2: The $({\rm ATT})_{F}$ (green/full curve) as derived via (10) with
$\tau_{y}$ in (11) (red/dashed curve) and $\tau_{z}$ in (12) (blue/dotted
curve) evaluated for 87Rb isotopes transmitted with average energy $E$ through
a rectangular barrier of height $V_{0}=135\ {\rm nK}$ and width $L_{0}=1.3\
{\rm\mu m}$. Also depicted is the classical time $\tau_{c}$ in (13)
(black/dot-dashed curve). For high barriers $(E/V_{0}\rightarrow 0)$,
$\tau_{y}$ vanishes, $\tau_{z}$ decreases but remains finite (Table 1) and, as
a result, the $({\rm ATT})_{F}$ blows up, showing that it takes forever to
tunnel though a sky-high barrier. In contrast to $({\rm ATT})_{F}$, the
classical time $\tau_{c}$ blows up as $E/V_{0}\rightarrow 1$. Figure 3: The
$({\rm ATT})_{F}$ (green/full curve) as derived via (10) with $\tau_{y}$ in
(11) (red/dashed curve) and $\tau_{z}$ in (12) (blue/dotted curve) evaluated
for 87Rb isotopes transmitted with average energy $E=V_{0}/2$ through a
rectangular barrier of height $V_{0}=135\ {\rm nK}$ and width $L$ (in units of
$L_{0}=1.3\ {\rm\mu m}$). Also depicted is the classical time $\tau_{c}$ in
(13) (black/dot-dashed curve). For wide barriers $(L/L_{0}\rightarrow\infty)$,
$\tau_{y}$ becomes width-independent, $\tau_{z}$ increases slightly (Table 1),
$\tau_{c}$ grows linearly and yet the $({\rm ATT})_{F}$ grows exponentially,
showing that it takes forever to tunnel though a worlds-wide barrier.

The way the $({\rm ATT})_{F}$, $\tau_{y}$ and $\tau_{z}$ approach to the high
and wide barrier limits in (14) are important not only for testing the models
but also for guiding future experiments (Table 1). Indeed, transcribing the
lab specs [19] of $V_{0}=135\ {\rm nK}$ and $L_{0}=1.3\ {\rm\mu m}$ to the
rectangular potential under concern, it is found that the $({\rm ATT})_{F}$
gradually increases and eventually diverges as $E/V_{0}\rightarrow 0$, and
ensures this way that high barriers are hard to tunnel through (Fig. 2). The
experimental $({\rm ATT})_{F}$ (Fig. 1), too, is expected to exhibit a similar
divergence if future experiments probe the deep tunneling regime $(E\ll
V_{0})$ beyond the available data [19]. In contrast to the $({\rm ATT})_{F}$,
$\tau_{y}$, the transmission time $({\rm ATT})_{S}$ [20, 21] in weak
measurement formalism [22], makes the unphysical prediction that tunneling
must take zero-time for macroscopic bodies and high barriers. In contrast, the
classical time $\tau_{c}$ blows up as $E/V_{0}\rightarrow 1$ kinematically
(Table 1 and Fig. 2).

The wide barrier limit (Table 1) and the way it is approached (Fig. 3), yet to
be probed by experiments like neutron tunneling through isotopic
nanostructures [39], reveal that the $({\rm ATT})_{F}$ involves slower-than-
light transmission, and is thus free of the Hartman effect [5, 40]. The
classical time $\tau_{c}$ grows linearly with the barrier width. In contrast,
however, the precession time $\tau_{y}$, identified with actual transmission
time in weak measurement approach [22, 20, 21], implies faster-than-light
transmission (Table 1), and exhibits therefore the Hartman effect [40]. The
resolution [30] of this unphysical behavior is that $\tau_{y}$ pertains to
group delay rather than wave-front delay [31, 32] and, since actual
transmission time must be set by the wave-front speed [32] faster-than-light
transmission implied by $\tau_{y}$ is of no physical consequence. Likewise,
faster-than-light transmission observed [41, 42, 43, 44] in evanescent
electromagnetic waves refers to group delay [30, 31, 32]. Thus, in view of its
physical behavior at the extremes of the potential (Table 1), and in view also
of its longevity compared to $\tau_{y}$ and $\tau_{z}$, the fluctuation-
induced time $({\rm ATT})_{F}$ possesses every reason to qualify as wave-front
delay [32].

##  4 Conclusion

In conclusion, we have shown that spin uncertainty, whose relevance is hinted
by the spin component along the magnetic field, shadows forth a new tunneling
time $({\rm ATT})_{F}$. This new time is a slower-than-light combination of
the two Larmor times $\tau_{y}$ and $\tau_{z}$. It remains valid for general
potential profiles, and can be extracted directly from the Larmor clock
readings. It is concordant with the wavefront delay, and respects therefore
fundamentals of the relativity and quantum. Needless to say, $({\rm ATT})_{F}$
gives the durations of tunneling-enabled phenomena in nature, which range from
nuclear fusion [1] in physics to smell perception [45] in chemistry to DNA
mutation [46, 47, 48, 49] in biology. Likewise, it sets operation speeds of
all tunneling-driven processes such as quantum annealing [50, 51] – a
tunneling-enabled quantum computation mechanism [52, 53, 26] which has already
been implemented on an 1,800-qubit chip [27] by D-Wave Systems. Quantum
annealing works in search spaces having numerous local minima [53, 26], and
getting out of each minimum costs a time delay given by the $({\rm ATT})_{F}$.
In this sense, the $({\rm ATT})_{F}$ pops up as an additional factor to be
taken into account in assessing the performance of quantum computers. This
example of quantum annealing can be extended to operation speeds of various
other tunneling-driven technological devices.

## Acknowledgements

This work is supported in part by the IPS Project B.A.CF-20-02239 at Sabancı
University. The author is grateful Ramón Ramos for sharing the experimental
data, G. Demir for computational help, and O. Sargın for discussions.

## References

  * [1] authorG. Gamow,  journalZ. Phys. volume51, pages204 (year1928).
  * [2] authorR. W. Gurney,  journalNature volume122, pages439 (year1928).
  * [3] authorL. A. MacColl,  journalPhys. Rev. volume40, pages621 (year1932).
  * [4] authorE. H. Hauge and authorJ. A. Støvneng,  journalRev. Mod. Phys. volume61, pages917 (year1989).
  * [5] authorV. S. Olkhovsky and authorE. Recami,  journalPhysics Reports volume214, pages339 (year1992).
  * [6] authorR. Landauer and authorT. Martin,  journalRev. Mod. Phys. volume66, pages217 (year1994).
  * [7] authorW. Pauli,  title _General Principles of Quantum Mechanics_ book(Springer, New York, 2012).
  * [8] authorG. Field,  journalPhilSci Archive http://philsci-archive.pitt.edu/ 18446/ (year2020).
  * [9] authorP. Busch,  journal“The time–energy uncertainty relation.” In Time in quantum mechanics, pp. 73-105 [arXiv:quant-ph/0105049] (yearSpringer, 2008).
  * [10] authorJ.Hilgevoord,  journalStud. Hist. Phil. Sci. B: Stud. Hist. Phil. Mod. Phys. volume36, pages29 (year2005).
  * [11] authorM. B. Altaie, authorD. Hodgson, and authorA. Beige,  journalFront. in Phys. volume10, pages897305 [arXiv:2203.12564 [quant-ph]] (year2022).
  * [12] authorE. P. Wigner,  journalPhys. Rev. volume98, pages145 (year1955).
  * [13] authorA. Baz’,  journalSov. J. Nucl. Phys. volume4, pages182 (year1966).
  * [14] authorV. Rybachenko,  journalSov. J. Nucl. Phys. volume5, pages635 (year1967).
  * [15] authorD. Demir and authorT. Güner,  journalAnnals of Physics volume386, pages291 (year2017).
  * [16] authorM. Büttiker,  journalPhys. Rev. B volume27, pages6178 (year1983).
  * [17] authorJ. P. Falck and authorH. D. Hauge,  journalPhys. Rev. B volume38, pages3287 (year1988).
  * [18] authorH. M. Krenzlin, authorJ. Budczies, and authorK. W. Kehr,  journalPhys. Rev. A volume53, pages3749 (year1996).
  * [19] authorR. Ramos, authorD. Spierings, authorI. Racicot, and authorA. M. Steinberg,  journalNature volume583, pages529 (year2020).
  * [20] authorA. M. Steinberg,  journalPhys. Rev. Lett. volume74, pages2405 (year1995).
  * [21] authorA. M. Steinberg,  journalPhys. Rev. A volume52, pages32 (year1995).
  * [22] authorY. Aharonov, authorD. Z. Albert, and authorL. Vaidman,  journalPhys. Rev. Lett. volume60, pages2325 (year1988).
  * [23] authorU. Fano,  journalPhys. Rev. volume72, pages26 (year1947).
  * [24] authorK. Rajdl, authorP. Lansky, and authorL. Kostal,  journalFront. Comput. Neurosci. volume14, pages569059 (year2020).
  * [25] authorD. R. Cox,  journalJ. Royal Stat. Soc. volumeXVII, pages129 (year1955).
  * [26] authorC. McGeoch,  journalTheor. Comp. Sci. volume816, pages169 (year2020).
  * [27] authorA. D. King, authorJ. Carrasquilla, authorJ. Raymond, authorI. Ozfidan, authorE. Andriyash, authorA. Berkley, authorM. Reis, authorT. Lanting, authorR. Harris, authorF. Altomare, authorK. Boothby, authorP. I. Bunyk, authorC. Enderud, authorA. Fréchette, authorE. Hoskinson, authorN. Ladizinsky, authorT. Oh, authorG. Poulin-Lamarre, authorC. Rich, authorY. Sato, author A. Y. Smirnov, authorL. J. Swenson, authorM. H. Volkmann, authorJ. Whittaker, authorJ. Yao, authorE. Ladizinsky, authorM. W. Johnson, authorJ. Hilton, and authorM. H. Amin,  journalNature volume560, pages456 (year2018).
  * [28] authorA. Elçi and authorH. P. Hjalmarson,  journalJ. Math. Phys. volume50, pages102101 (year2009).
  * [29] authorM. Büttiker and authorR. Landauer,  journalPhys. Rev. Lett. volume49, pages1739 (year1982).
  * [30] authorR. Chiao, _Superluminal phase and group velocities (….)_ ,  arXiv:1111.2402 [physics.class-ph] (2011).
  * [31] authorL. Brillouin and authorH. S. W. Massey, _Wave propagation and group velocity_ book(Elsevier, Saint Louis, 2014).
  * [32] authorJ. L. Leander,  journalJ. Acous. Soc. Am. volume100, pages3503 (year1996).
  * [33] authorA. S. Landsman, authorM. Weger, authorJ. Maurer, authorR. Boge, authorS. Ludwig, authorC. Cirelli, authorL. Gallmann, and authorU. Keller,  journalOptica volume1, pages343 (year2014).
  * [34] authorN. Camus, authorE. Yakaboylu, authorL. Fechner, authorM. Klaiber, authorM. Laux, authorY. Mi, author K. Z. Hatsagortsyan, authorT. Pfeifer, author C. H. Keitel, and authorR. Moshammer,  journalPhys. Rev. Lett. volume119, pages023201 (year2017).
  * [35] authorU. S. Sainadh, authorH. Xu, authorX. Wang, authorA. Atia-Tul-Noor, authorW. C. Wallace, authorN. Douguet, authorA. Bray, authorI. Ivanov, authorK. Bartschat, authorA. Kheifets, authorR. T. Sang, and authorI. V. Litvinyuk,  journalNature volume568, pages75 (year2019).
  * [36] authorS. Potnis, authorR. Ramos, authorK. Maeda, authorL. D. Carr, and authorA. M. Steinberg,  journalPhys. Rev. Lett. volume118, pages060402 (year2017).
  * [37] authorX. Zhao, authorD. A. Alcala, authorM. A. McLain, authorK. Maeda, authorS. Potnis, authorR. Ramos, authorA. M. Steinberg, and authorL. D. Carr,  journalPhys. Rev. A. volume96, pages063601 (year2017).
  * [38] authorR. Ramos, authorD. Spierings, authorS. Potnis, and authorA. M. Steinberg,  journalPhys. Rev. A volume98, pages23611 (year2018).
  * [39] authorA. Matiwane, authorJ. Sackey, authorM. L. Lekala, and authorM. Maaza,  journalMRS Advances volume3, pages2609 (year2018).
  * [40] authorT. E. Hartman,  journalJ. Appl. Phys. volume33, pages3427 (year1962).
  * [41] authorA. Ranfagni, authorD. Mugnai, authorP. Fabeni, and authorG. P. Pazzi,  journalAppl. Phys. Lett. volume58, pages774 (year1991).
  * [42] authorA. Enders and authorG. Nimtz,  journalJ. Phys. I volume2, pages1693 (year1992).
  * [43] authorA. M. Steinberg, authorP. G. Kwiat, and authorR. Y. Chiao,  journalPhys. Rev. Lett. volume71, pages708 (year1993).
  * [44] authorC. Spielmann, authorR. Szipöcs, authorA. Stingl, and authorF. Krausz,  journalPhys. Rev. Lett. volume73, pages2308 (year1994).
  * [45] authorE. R. Bittner, authorA. Madalan, authorA. Czader, and authorG. Roman,  journalJ. Chem. Phys. volume137, pages551 (year2012).
  * [46] authorP. O. Löwdin,  journalRev. Mod. Phys. volume35, pages724 (year1963).
  * [47] authorA. D. Godbeer, authorJ. S. Al-Khalili, and authorP. D. Stevenson,  journalPhys. Chem. Chem. Phys. volume17, pages13034 (year2015).
  * [48] authorG. Çelebi, authorE. Özçelik, authorE. Vardar, and authorD. Demir,  journalProg. Biophys. Mol. Biol. volume167, pages96 (year2021).
  * [49] authorE. Özçelik, authorE. Akar, authorS. Zaman, and authorD. Demir,  journalProg. Biophys. Mol. Biol. pageshttps://doi.org/10.1016/j.pbiomolbio.2022.05.009 (year2022).
  * [50] authorB. Apolloni, authorM. C. Carvalho, and authorD. De Falco,  journalStoc. Proc. Appl. volume33, pages233 (year1989).
  * [51] authorT. Kadowaki and authorH. Nishimori,  journalPhys. Rev. E volume58, pages5355 (year1998).
  * [52] authorE. Farhi, authorJ. Goldstone, authorA. Ludgren, and authorD. Preda,  journalScience volume292, pages472 (year2001).
  * [53] authorB. Heim, authorT. Rønnow, authorS. V. Isakov, and authorM. Troyer,  journalScience volume348, pages215 (year2015).