# On the Masses, Age, & Architecture of the VHS J1256$-$1257AB b System

Trent J. Dupuy1 , Michael C. Liu2 , Elise L. Evans1, William M. J. Best3 ,
Logan A. Pearce4 , Aniket Sanghi2,3 , Mark W. Phillips2 , and Daniella C.
Bardalez Gagliuffi5
1Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh, EH9 3HJ, UK
2Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu,
HI 96822, USA
3Department of Astronomy, The University of Texas at Austin, 2515 Speedway
C1400, Austin, TX 78712, USA
4Steward Observatory, University of Arizona, Tucson, AZ 85721, USA
5American Museum of Natural History, 200 Central Park West, New York, NY
10024, USA E-mail<EMAIL_ADDRESS>Graduate Research Fellow
0000-0001-9823-1445 0000-0003-2232-7664 0000-0003-0562-1511
0000-0003-3904-7378 0000-0002-1838-4757 0000-0001-6041-7092
0000-0001-8170-7072

###### Abstract

VHS J1256$-$1257 AB is an ultracool dwarf binary that hosts a wide-separation
planetary-mass companion that is a key target of the JWST Exoplanet Early
Release Science (ERS) program. Using Keck adaptive optics imaging and aperture
masking interferometry, we have determined the host binary’s orbit ($a=1.96\pm
0.03$ au, $P=7.31\pm 0.02$ yr, $e=0.8826^{+0.0025}_{-0.0024}$) and measured
its dynamical total mass ($0.141\pm 0.008$ $M_{\sun}$). This total mass is
consistent with VHS J1256$-$1257 AB being a brown dwarf binary or pair of very
low-mass stars. In addition, we measured the orbital motion of VHS
J1256$-$1257 b with respect to the barycenter of VHS J1256$-$1257 AB, finding
that the wide companion’s orbit is also eccentric ($e=0.73^{+0.09}_{-0.10}$),
with a mutual inclination of $116\degree\pm 16\degree$ with respect to the
central binary. This orbital architecture is consistent with VHS J1256$-$1257
b attaining a significant mutual inclination through dynamical scattering and
thereafter driving Kozai-Lidov cycles to pump the eccentricity of VHS
J1256$-$1257 AB. We derive a cooling age of $140\pm 20$ Myr for VHS
J1256$-$1257 AB from low-mass stellar/substellar evolutionary models. At this
age, the luminosity of VHS J1256$-$1257 b is consistent with both deuterium-
inert and deuterium-fusing evolutionary tracks. We thus find a bimodal
probability distribution for the mass of VHS J1256$-$1257 b, either $11.8\pm
0.2$ $M_{\rm Jup}$ or $16\pm 1$ $M_{\rm Jup}$, from Saumon & Marley (2008)
hybrid models. Future spectroscopic data to measure isotopologues such as HDO
and CH3D could break this degeneracy and provide a strong test of substellar
models at the deuterium-fusion mass boundary.

###### keywords:

astrometry – binaries: visual – planetary systems

††pubyear: 2022††pagerange: On the Masses, Age, & Architecture of the VHS
J1256$-$1257AB b System–On the Masses, Age, & Architecture of the VHS
J1256$-$1257AB b System

##  1 Introduction

Two fundamental parameters govern the bulk properties of gas-giant planets and
brown dwarfs: mass and age. Masses are difficult to measure directly for
imaged planets because of their long orbital periods, although there has been
progress for a few planets inside of 20 au (e.g., Snellen & Brown, 2018; Dupuy
et al., 2019; Dupuy et al., 2022; Brandt et al., 2021b, d). Ages for imaged
planets have so far been largely reliant on an object belonging to a well-
studied young association because of the limitations of determining precise
ages for field stars. One notable exception is the Y-dwarf exoplanet WD
0806-661b (Luhman et al., 2011) where its age of $1.5^{+0.5}_{-0.3}$ Gyr is
determined by the cooling age of its white-dwarf host star.

There are also a small handful of gas-giant companions that could have
precisely measured substellar cooling ages. Such ages are similar to white-
dwarf cooling ages but without the need to estimate a stellar progenitor’s
lifetime. Because both stars and brown dwarfs begin with an initial entropy
that is related to their mass, the age of such an object can be determined by
measuring its mass and present-day luminosity (e.g., Liu et al., 2008;
Burgasser, 2009). Gas giants in systems where the host stars are binary brown
dwarfs or pre-main–sequence stars are amenable to such cooling-age
measurements. VHS J125601.92$-$125723.9 (hereinafter VHS J1256$-$1257) is one
of the few such systems.

Gauza et al. (2015) used data from the Visible and Infrared Survey Telescope
for Astronomy (VISTA) Hemisphere Survey (VHS) to discover that VHS
J125601.58$-$125730.3 (hereinafter VHS J1256$-$1257 b) is a common-proper
motion companion to VHS J1256$-$1257 A, at a projected separation of
$8\aas@@fstack{\prime\prime}06$ (178 au). They derived spectral types and
gravity classifications for the host and companion of M$7.5\pm 0.5$ int-g and
L$7.0\pm 1.5$ vl-g, respectively. They also measured a parallax of $78.8\pm
6.4$ mas, which placed the companion in the same location as HR 8799 b on the
color-magnitude diagram. The parallax has since been updated, first by the
Hawai‘i Infrared Parallax Program ($45.0\pm 2.4$ mas; Dupuy et al., 2020) and
most recently by Gaia EDR3 ($47.27\pm 0.47$ mas = $21.14\pm 0.22$ pc; Gaia
Collaboration et al., 2016, 2021). The system is therefore more distant than
originally thought and no longer a direct analog to HR 8799 b, but given the
primary’s young age (150–300 Myr; Dupuy et al., 2020), the companion is still
potentially planetary mass. Its cool temperature and wide separation make VHS
J1256$-$1257 b appealing for direct studies, including being the primary
spectroscopy target for JWST’s Exoplanet Early Release Science program
(Hinkley et al., 2022).

Adaptive optics (AO) imaging has revealed that the host VHS J1256$-$1257 A is
itself a tight binary (Stone et al., 2016; Rich et al., 2016), making VHS
J1256$-$1257 a rare triple system potentially composed entirely of substellar
objects. The inner binary’s projected separation at discovery was $2.62\pm
0.03$ au (using the latest distance), which would correspond to a
$\sim$10-year orbital period. The likelihood of obtaining dynamical masses on
such a relatively short time scale motivated us to begin an orbit monitoring
campaign, as we have done for other substellar binaries (e.g., Dupuy & Liu,
2017). We present here the dynamical masses for the binary components and a
corresponding cooling age that suggests the directly imaged companion VHS
J1256$-$1257 b may be below the deuterium-fusing mass limit.

Table 1: Keck/NIRC2 Relative Astrometry of VHS J1256$-$1257 AB. Epoch | Sep (mas) | PA (°) | Corr |  $\Delta{K_{\rm MKO}}$ (mag)
---|---|---|---|---
2016.059 | $128.21\pm 0.14$ | $168.07\pm 0.05$ | $0.01$ | $0.021\pm 0.006$
2017.050 | $139.7\pm 2.8$ | $163.5\pm 1.9$ | $-0.20$ |  $-0.06\pm 0.16$
2017.220 | $140.5\pm 0.7$ | $162.9\pm 0.4$ | $0.35$ | $0.031\pm 0.023$
2018.017 | $135.84\pm 0.29$ | $158.37\pm 0.15$ | $0.83$ | $0.026\pm 0.026$
2019.259 | $109.6\pm 0.4$ | $150.36\pm 0.11$ | $0.32$ | $0.032\pm 0.028$
2021.018 |  $35.02\pm 0.26$ | $112.5\pm 0.6$ | $0.00$ | $0.033\pm 0.005$
2022.066 |  $70.0\pm 1.8$ | $181.0\pm 1.7$ | $0.64$ |  $0.19\pm 0.04$
2022.271 |  $85.6\pm 0.7$ | $177.2\pm 0.4$ | $-0.52$ |  $0.25\pm 0.05$

##  2 Observations

We obtained astrometry for the VHS J1256$-$1257 system from the Keck II
Telescope using the facility AO system. We began monitoring on 2016 Jan 22 UT
using the Maunakea Observatories $K$-band filter (Tokunaga et al., 2002) and
NIRC2’s narrow-camera mode, which has a pixel scale of $9.971\pm 0.004$ mas
pixel-1(Service et al., 2016) and field-of-view of
$10\aas@@fstack{\prime}2\times 10\aas@@fstack{\prime}2$. Most of our
measurements were made using the standard laser guide star (LGS) AO system
(Wizinowich et al., 2006), which uses a Shack-Hartmann wavefront sensor to
measure the LGS and a separate red-optical sensor observing VHS J1256$-$1257
AB itself as a tip-tilt reference. At one epoch, 2022 Jan 24 UT, we instead
used the infrared pyramid wavefront sensor (Bond et al., 2020) for natural
guide star AO correction using VHS J1256$-$1257 AB.

Our first imaging from 2016 was obtained less than a year after the discovery
imaging using MagAO/Clio2 and NIRC2 from Stone et al. (2016), and it was
consistent with their measurement of increasing projected separation. In the
following, we use the earliest data that comes from MagAO along with our own
NIRC2 data. By 2018 the projected separation began decreasing and eventually
the binary was unresolved on 2021 Jan 5 UT. We obtained aperture masking
interferometry data the following night, using the 9-hole mask (Ireland et
al., 2008), and successfully resolved VHS J1256$-$1257 AB at 35 mas
separation. In addition, starting with our first observation in 2016, we
regularly obtained our imaging in a way that captured both VHS J1256$-$1257 AB
as well as VHS J1256$-$1257 b. This allowed us to measure astrometry for VHS
J1256$-$1257 AB and the orbital motion of VHS J1256$-$1257 b.

Our methodology for reducing NIRC2 imaging and masking data is described
extensively in our previous work (e.g., Dupuy & Liu, 2017). Briefly, we
perform standard calibrations (dark subtraction and flat-fielding with dome
flats) and then measure the separation, position angle (PA), and flux ratio in
individual images. This is done using StarFinder (Diolaiti et al., 2000) when
possible, but when this fails at the closest separations we use an analytical,
multi-component Gaussian PSF model optimized using the Levenberg-Marquardt
algorithm implemented in IDL by the mpfit routine (Markwardt, 2009). We
correct our measured pixel positions for NIRC2’s distortion using the Service
et al. (2016) astrometric calibration, which also provides the pixel scale and
orientation of the images. Final measurements at an epoch are the mean and
standard deviation of values from individual images. For masking data, we
obtain binary parameters by fitting the closure phases using the Sydney
pipeline (Ireland et al., 2008).111https://github.com/mikeireland/idlnrm

Table 1 presents our relative astrometry for VHS J1256$-$1257 AB, including
the linear Pearson correlation coefficient (Corr) for separation and PA. Table
2 presents astrometry of VHS J1256$-$1257 b we derived from imaging epochs
where VHS J1256$-$1257 AB was sufficiently well resolved for StarFinder
analysis and that conformed to a standard configuration with the NIRC2
$y$-axis ${\rm PA}\approx 0$° and VHS J1256$-$1257 AB at NIRC2
$(x,y)\approx(250,800)$ pix. By keeping all three components in approximately
the same location on NIRC2, the astrometry should be minimally impacted by the
$\approx$1 mas uncertainty in the distortion solution. We follow convention in
referring to relative declination as $\Delta\delta$ and right ascension as
$\Delta\alpha^{*}\equiv\Delta\alpha\cos\delta$.

Table 2: Keck/NIRC2 Relative Astrometry of VHS J1256$-$1257 b.

Epoch |  $\Delta{\alpha^{*}}_{\rm b-A}$ (mas) |  $\Delta{\delta}_{\rm b-A}$ (mas)
---|---|---
2016.059 | $-4974.0\pm 2.8$ | $-6409.8\pm 2.3$
2017.220 | $-4977.6\pm 2.2$ | $-6412.8\pm 2.1$
2018.017 | $-4982.9\pm 2.5$ | $-6406.5\pm 2.6$
2022.271 | $-5049.0\pm 1.3$ | $-6387.5\pm 1.1$

  * Note.

– Astrometry relative to VHS J1256$-$1257 A not the barycenter of AB. As
described in Section 4, we determine the position of VHS J1256$-$1257 b
relative to the VHS J1256$-$1257 AB barycenter to be
$\Delta(\alpha^{*},\delta)_{\rm b-AB}=(-5025.6\pm 2.0,-6350\pm 9)$ mas at the
mean epoch 2019.853.

Our multi-epoch NIRC2 data precisely constrain the $K$-band flux ratio of VHS
J1256$-$1257 AB. We used only StarFinder and masking results, as these should
be less prone to systematic errors (epochs 2017.05, 2018.02, 2019.26, and
2021.02). The flux ratios are in excellent agreement, with $\chi^{2}=0.44$ and
3 degrees of freedom (dof), so we adopt the weighted average $\Delta{K_{\rm
MKO}}=0.033\pm 0.004$ mag.

The integrated-light spectrum of VHS J1256$-$1257 AB was obtained with
IRTF/SpeX in prism mode on 2016 Feb 19 UT as part of NASA IRTF program
2016A079 (PI: Bardalez Gagliuffi). The target was observed at airmass 1.19
with the 0$\aas@@fstack{\prime\prime}$5 slit and $6\times 60$ s exposures. The
A0 star HD 112304 was observed immediately after the target for flux
calibration and telluric correction. Internal flat fields and argon arc lamps
followed the standard observations for wavelength calibration. All data were
reduced with the IDL package SpeXtool v4.1. Further details on the
observations and instrument settings can be found in Bardalez Gagliuffi et al.
(2014) and Burgasser et al. (2010).

##  3 Luminosities

We computed the combined-light bolometric luminosity of VHS J1256$-$1257 AB by
direct integration of its unresolved optical to mid-infrared (MIR) spectral
energy distribution (SED). Our assembled SED consists of available Pan-
STARRS-1 (PS1; Chambers et al., 2016) optical photometry ($g$, $r$, $y$), the
near-infrared (NIR) IRTF/SpeX prism spectrum from Section 2, NIR photometry
from 2MASS (Cutri et al., 2003), and MIR photometry from the CatWISE catalog
($W1$ and $W2$ bands; Eisenhardt et al., 2020; Marocco et al., 2021) and
AllWISE catalog ($W3$ and $W4$ bands; Cutri et al., 2013). We began by flux-
calibrating the SpeX spectrum using the weighted average of calibrations
derived from PS1 $y$ and 2MASS $JHKs$ photometry, assuming a systematic noise
floor of 0.01 mag for all the filters. We then integrated the flux-calibrated
SpeX spectrum to determine the NIR contribution to the bolometric flux, with
an error that accounts for the uncertainties in the spectral data points and
the overall flux calibration. We determined the optical and MIR contributions
to the bolometric flux by simultaneously fitting BT-Settl model atmospheres
(CIFIST2011/2015; Allard et al., 2012; Baraffe et al., 2015) to the PS1 and
WISE photometry (computing synthetic photometry from the models) and the SpeX
spectrum (with the models degraded to the non-linear spectral resolution of
the 0$\aas@@fstack{\prime\prime}$5 slit). We found the best-fitting BT-Settl
model had $\mbox{$T_{\rm eff}$}=2700$ K and $\mbox{$\log(g)$}=5.0$ dex. Our
final bolometric flux was found by adding the NIR contribution to the
integration of the model outside the wavelength range of the SpeX spectrum.
The uncertainty in the optical+MIR contribution was obtained from the standard
deviation of the corresponding measurements derived using the four model
spectra adjacent in $T_{\rm eff}$ and $\log(g)$ to the best-fitting model. Our
final bolometric flux of VHS J1256$-$1257 AB is $1.47\pm 0.04\times 10^{-13}$
W m-2. Using its parallactic distance of $21.14\pm 0.22$ pc, we calculated a
bolometric luminosity $\log((\mbox{$L_{\rm bol}$}_{\rm,A}+\mbox{$L_{\rm
bol}$}_{\rm,B})/\mbox{$L_{\sun}$})=-2.687\pm 0.021$ dex.

To derive component luminosities for VHS J1256$-$1257 AB, we used the
empirical relation between $K_{\rm 2MASS}$-band absolute magnitude and $L_{\rm
bol}$ from Dupuy & Liu (2017). We assume that $\Delta{K_{\rm
MKO}}=\Delta{K_{\rm 2MASS}}$ here because of the near-unity flux ratio. Using
a Monte Carlo method, we drew random absolute magnitudes representative of VHS
J1256$-$1257 A ($9.1\pm 0.2$ mag, truncated at 8.7 mag, the upper limit of the
empirical relation). We then simulated VHS J1256$-$1257 B by adding $0.033\pm
0.004$ mag to this absolute magnitude and computed the difference in derived
$L_{\rm bol}$ values from the relation. We found $\log(\mbox{$L_{\rm
bol}$}_{\rm,B}/\mbox{$L_{\rm bol}$}_{\rm,A})=-0.012\pm 0.002$ dex. We
therefore calculated component luminosities of $\log(\mbox{$L_{\rm
bol}$}_{\rm,A}/\mbox{$L_{\sun}$})=-2.982\pm 0.021$ dex and $\log(\mbox{$L_{\rm
bol}$}_{\rm,B}/\mbox{$L_{\sun}$})=-2.994\pm 0.021$ dex.

In the following, for the luminosity of VHS J1256$-$1257 b, we use
$\log(\mbox{$L_{\rm bol}$}_{\rm,b}/\mbox{$L_{\sun}$})=-4.63\pm 0.03$ dex,
which is the value derived by the atmosphere model analysis of Hoch et al.
(2022) and adjusted to the more precise Gaia parallax.

##  4 Orbit analysis

For the three-body system of VHS J1256$-$1257 AB b, we separate our orbital
analysis into the relative orbit of the host binary ($\sim$2 au) and the
orbital motion of the wide companion ($\sim$170 au) relative to the AB
barycenter. Dynamical interactions are negligible for such a wide, low mass-
ratio ($M_{\rm comp}/M_{\rm host}<0.1$) system.

To fit the relative orbit of VHS J1256$-$1257 AB we used orvara (v1.0.4;
Brandt et al., 2021c). orvara utilizes a novel, highly-efficient eccentric
anomaly solver and determines posteriors of orbital parameters using the
affine-invariant (Goodman & Weare, 2010) Markov-Chain Monte Carlo (MCMC)
sampler emcee (Foreman-Mackey et al., 2013) with parallel-tempering (Vousden
et al., 2016). We provide our orvara configuration files as supplementary data
here, but briefly, we fitted all eight standard parameters for a relative
astrometric fit with their default priors (linear-flat in eccentricity $e$ and
viewing angles, except inclination $p(i)\propto\sin{i}$, and log-flat in mass
and semimajor axis $a$). Relative astrometry only constrains the total mass,
$\mbox{$M_{\rm tot}$}\equiv M_{\rm A}+M_{\rm B}$, so we placed no limiting
priors on the component masses. Thus, $M_{\rm A}$ and $M_{\rm B}$ varied
freely in the orvara MCMC analysis, but were always constrained implicitly to
follow a consistent $M_{\rm tot}$. Our results are based on a run with 100
walkers and $10^{5}$ steps for the MCMC and 5 temperatures for parallel
tempering. We thinned our chains, retaining every 50th step, and discarded the
first 50% as burn-in, yielding $10^{5}$ final samples in our posterior.

Figure 1: Relative astrometry of VHS J1256$-$1257 AB from our Keck/NIRC2 data
(blue hexagons; Table 1) and MagAO/Clio2 discovery imaging (red circles).
Left: Positions of VHS J1256$-$1257 B (colored symbols; error bars are smaller
than the symbols) relative to VHS J1256$-$1257 A (black star). The maximum-
likelihood orbit (thick black line) is shown with its line of nodes (dotted
grey line) and time of periastron passage (solid grey line and arrow)
indicated. Thin lines color-coded by the total mass show 50 randomly drawn
orbits drawn from our MCMC posterior. Right: The same measurements and orbits
shown as separation (top) and PA (bottom) as a function of time with subpanels
showing the residuals with respect to the maximum-likelihood orbit.  Table 3:
VHS J1256$-$1257 AB orbital parameters derived from orvara.

Property | Median $\pm$1$\sigma$ | 95.4% c.i.
---|---|---
Total mass, $M_{\rm tot}$ ($M_{\sun}$) | $0.141\pm 0.008$ | 0.125, 0.157
Semimajor axis, $a$ (au) | $1.96\pm 0.03$ | 1.89, 2.03
Eccentricity, $e$ | $0.8826_{-0.0024}^{+0.0025}$ | 0.8776, 0.8875
Inclination, $i$ (°) |  $118.7\pm 1.0$ | 116.7, 120.7
PA of ascending node, $\Omega$ (°) | $4.4\pm 0.5$ | 3.5, 5.3
Argument of periastron, $\omega$ (°) |  $44.9\pm 1.0$ | 42.8, 46.9
Mean longitude at $t_{\rm ref}$, $\lambda_{\rm ref}$ (°) |  $-163.5\pm 1.6$ |  $-$166.7, $-$160.1
Period, $P$ (yr) | $7.307_{-0.024}^{+0.023}$ | 7.262, 7.357
Time of periastron, $t_{p}$ (yr) |  $2021.537_{-0.014}^{+0.015}$ | 2021.507, 2021.566
$\tau\equiv(t_{p}-t_{\rm ref})/P$ | $0.579\pm 0.007$ | 0.565, 0.592

  * *

Reference epoch $t_{\rm ref}=2010.0$ (55197 MJD).

  * Note.

Free parameters in the MCMC are shown in the top section. These were used to
compute the parameters in the bottom section. All posterior distributions are
nearly Gaussian.

Our measured total mass of $0.141\pm 0.008$ $M_{\sun}$ suggests that the
components of VHS J1256$-$1257 AB are possibly brown dwarfs with masses of
$74\pm 4$ $M_{\rm Jup}$ if their mass ratio is near unity. Their orbital
eccentricity of $0.8826^{+0.0025}_{-0.0024}$ is the highest ever measured for
a very low-mass binary (e.g., see Dupuy & Liu, 2011, 2017).

Independent of the orbit analysis of VHS J1256$-$1257 AB we measured the
orbital motion of VHS J1256$-$1257 b relative to its host’s barycenter
(denoted b–AB). As mentioned in Section 2, at some epochs we obtained imaging
of all three objects in individual NIRC2 images. The relative positions of A
and B are typically measured $\sim$10$\times$ more precisely than A or B to
the companion b. This allowed us to approximate the errors in VHS J1256$-$1257
AB relative astrometry to be negligible compared to those of VHS J1256$-$1257
b. Under this assumption, the position of VHS J1256$-$1257 b relative to VHS
J1256$-$1257 A can be written as

$\hfill\Delta{\alpha^{*}}_{\rm b-A}=\Delta{\alpha^{*}}_{\rm
b-AB}+(\mu_{\alpha^{*},{\rm b-AB}}\times t)-[(M_{\rm A}/\mbox{$M_{\rm
tot}$})\times\Delta{\alpha^{*}}_{\rm A-B}]\hfill$ (1) $~{}\Delta{\delta}_{\rm
b-A}=\Delta{\delta}_{\rm b-AB}+(\mu_{\delta,{\rm b-AB}}\times t)-[(M_{\rm
A}/\mbox{$M_{\rm tot}$})\times\Delta{\delta}_{\rm A-B}],$ (2)

where the left-hand side corresponds to the measurements in Table 2, the
$\Delta_{\rm A-B}$ values on the far right side can be derived from Table 1,
and the rest are free parameters.

We used mpfit to find the best-fit solution and then used a Monte Carlo
approach to derive the uncertainties in the fit by randomly drawing simulated
measurements from the best-fit model with scatter equal to the individual
input measurements. Unfortunately, the mass ratio is poorly constrained in
this analysis ($M_{\rm A}/\mbox{$M_{\rm tot}$}=0.45\pm 0.08$), likely due to
the eccentric orbit and the fact that the measurements used here (i.e., when
all three components are well resolved) happen to come from a similar phase of
the orbit. In contrast, the orbital motion of the companion relative to the
VHS J1256$-$1257 AB barycenter is well detected at $\mu_{\alpha^{*},{\rm
b-AB}}=-10.7\pm 0.6$ mas yr-1 and $\mu_{\delta,{\rm b-AB}}=-0.6\pm 0.7$ mas
yr-1.

To fit the orbit of VHS J1256$-$1257 b relative to the AB barycenter we used
the python package lofti_gaia (Pearce et al.,
2020).222https://github.com/logan-pearce/lofti_gaia lofti_gaia is based on the
Orbits-For-The-Impatient (OFTI; Blunt et al., 2017) rejection-sampling method
and fits orbital parameters of resolved binaries in Gaia using their proper
motions and radial velocities if available. Here we adopted the architecture
of lofti_gaia to use our measured proper motion of b relative to the AB
barycenter at the mean observation epoch 2019.853, rather than Gaia astrometry
at the mean Gaia epoch. We fitted six orbital parameters: semimajor axis
($a_{\rm b}$), eccentricity ($e_{\rm b}$), inclination ($i_{\rm b}$), argument
of periastron ($\omega_{\rm b}$), longitude of ascending node ($\Omega_{\rm
b}$), and time of periastron passage ($t_{\rm p,b}$). Total system mass and
distance were drawn from normal distributions of $0.152\pm 0.010$ $M_{\sun}$
and $21.14\pm 0.22$ pc. This system mass is based on our measured mass for VHS
J1256$-$1257 AB and an estimated mass of $0.011\pm 0.006$ $M_{\sun}$ for VHS
J1256$-$1257 b from our evolutionary model analysis in Section 5. OFTI
rejection sampling generates trial orbits by drawing random values for four
orbital parameters from priors in $e_{\rm b}$: Uniform on [0,1); $\cos(i_{\rm
b})$: Uniform on [-1,1]; $\omega_{\rm b}$: Uniform on [0,2$\pi$]; orbit phase,
$(t_{\rm p,b}-2019.853)/P_{\rm b}$: Uniform on [0,1]. OFTI then scales the
semimajor axis and rotates $\Omega_{\rm b}$ to match the input data and
determines whether to accept or reject a trial by comparing its proper motion
in RA and Dec to our measured values. There is no prior on $a_{\rm b}$ or
$\Omega_{\rm b}$.

We ran lofti_gaia on our measured proper motions until $10^{5}$ trial orbits
were accepted. Table 4 reports the output probability distributions of orbital
parameters of VHS J1256$-$1257 b around its host, and Figure 2 shows these
orbits on the sky.

Figure 2: The orbit of VHS J1256$-$1257 b (orange circle) relative to its host
binary (brown circles; separation not to scale). Shown are 50 randomly drawn
orbits from our lofti_gaia analysis, color-coded by eccentricity.  Table 4:
VHS J1256$-$1257 b orbital parameters derived from lofti_gaia.

Property | Median $\pm$1$\sigma$ | 95.4% c.i.
---|---|---
Semimajor axis $a_{\rm b}$ (au) | $360_{-150}^{+110}$ |  170, 1040
Eccentricity $e_{\rm b}$ | $0.73_{-0.10}^{+0.09}$ | 0.56, 0.93
Inclination $i_{\rm b}$ (°) | $26_{-16}^{+12}$ |  4, 51
Argument of periastron $\omega_{\rm b}$ (°) | $180_{-140}^{+80}$ |  0, 330
PA of the ascending node $\Omega_{\rm b}$ (°) | $40_{-160}^{+50}$ |  $-$140, 200
Time of periastron $t_{p},{\rm b}$ (yr) |  $1250\pm 80$ |  990, 1430
Period $P_{\rm b}$ (kyr) | $17_{-11}^{+8}$ |  5, 85
$\tau_{\rm b}\equiv(t_{p}-2019.853)/P$ | $0.046_{-0.036}^{+0.022}$ | 0.000, 0.108

  * Note.

Our lofti_gaia analysis adopted a system mass of $M_{\rm A}+M_{\rm B}+M_{\rm
b}=0.152\pm 0.010$ $M_{\sun}$. Free parameters from the fit are shown in the
top section. These were used to compute the parameters in the bottom section.

##  5 Evolutionary model analysis

Substellar objects with well-determined luminosities enable precise
evolutionary model-derived cooling ages (when mass is known) and masses (when
age is known). Some key aspects of evolutionary models are quite uncertain,
such as the treatment of clouds. The relatively sparse tests of models with
objects of known mass, age, and luminosity have found a mixed bag of agreement
and potential problems (e.g., Dupuy et al., 2014; Beatty et al., 2018; Brandt
et al., 2021a), so we note that any mass or age derived from evolutionary
models should be treated with corresponding uncertainty. In the following, we
use our dynamical mass measurement of VHS J1256$-$1257 AB to determine a
substellar cooling age for the system and then use this cooling age to
estimate the mass of VHS J1256$-$1257 b.

For VHS J1256$-$1257 AB, the most appropriate evolutionary models are from
Baraffe et al. (2015). As in our previous work (e.g., Dupuy & Liu, 2017), we
used a Monte Carlo rejection-sampling approach to derive an age probability
distribution from input luminosity and mass prior distributions. We assumed a
linear-flat prior in age and a log-flat prior in $M_{\rm A}$, drawing random,
uniformly distributed values, while simultaneously drawing random values of
$M_{\rm tot}$ from our MCMC posterior. We calculated $M_{\rm B}$ as the
difference between $M_{\rm tot}$ and $M_{\rm A}$. For each age-mass pair, for
each component, we computed a model luminosity from bilinear interpolation of
the model grid. The probability of any sample being accepted was
$p=e^{-(\chi^{2}-\chi_{\rm min}^{2})/2}$, where $\chi^{2}$ was computed as the
sum of comparing our measured luminosities to the model-calculated ones, and
$\chi^{2}_{\rm min}$ was the lowest value among the ensemble of trial values.
A sample was accepted if a randomly drawn number $0<u<1$ for a given trial
satisfied $p>u$. We then computed other model-derived properties, such as
$T_{\rm eff}$, using the accepted mass and age samples.

We found a cooling age of $140\pm 20$ Myr for VHS J1256$-$1257 AB, with an
approximately Gaussian probability distribution, which is consistent with the
nondetection of lithium in its spectrum. Extreme lithium depletion of
$>10^{-4}$ corresponds to the older part of the age range at $<2$$\sigma$
according to the Baraffe et al. (2015) models.

We used VHS J1256$-$1257 AB’s age posterior to perform a rejection-sampling
analysis of the companion’s properties. The only evolutionary model grid that
reaches VHS J1256$-$1257 b’s luminosity and accounts for cloud evolution is
the “hybrid” grid of Saumon & Marley (2008). As seen in Figure 3, Saumon &
Marley (2008) models predict that objects with the luminosity and age of VHS
J1256$-$1257 b should be rare because they fall in a gap between low-mass
objects that cannot fuse deuterium and more massive objects that are either
fusing deuterium now (and thus more luminous at this age) or have already
fused their deuterium (and thus are older at this luminosity). Our rejection
sampling analysis correspondingly results in a bimodal posterior distribution.
The more probable outcome, with 70% of the posterior, is that VHS J1256$-$1257
b is a deuterium-bearing object of $11.8\pm 0.2$ $M_{\rm Jup}$. The less
probable outcome is that VHS J1256$-$1257 b has already depleted its deuterium
and is an object of $16\pm 1$ $M_{\rm Jup}$.

Figure 3: Left: evolutionary models from Saumon & Marley (2008) predict
$L_{\rm bol}$ as a function of age for objects of varying masses. (Light and
dark colored lines alternate between intervals of 0.01 $M_{\sun}$, and solid
and dashed lines alternate between even and odd masses.) Orange contours
indicate the 1$\sigma$ and 2$\sigma$ joint constraints on the luminosity and
age of VHS J1256$-$1257 b, which is just to the right of the 0.011 $M_{\sun}$
(11.5 $M_{\rm Jup}$) track. It is also consistent with more massive tracks
that are no longer D fusing. Objects with masses 0.012–0.013 $M_{\sun}$ are
predicted to be much more luminous at this age due to a D-fusion spike. Right:
posterior distribution of the mass of VHS J1256$-$1257 b from our rejection
sampling Monte Carlo analysis of these models. The distribution is bimodal,
with a 70/30 split, respectively, between a D-inert mass of $11.8\pm 0.2$
$M_{\rm Jup}$ and D-fusing mass of $16\pm 1$ $M_{\rm Jup}$. The lower mass
corresponds to a younger age and thus larger radius, lower $T_{\rm eff}$, and
lower $\log(g)$.

The radius of VHS J1256$-$1257 b, according to these models, is 1.27 $R_{\rm
Jup}$ in the lower-mass scenario and 1.19 $R_{\rm Jup}$ in the higher-mass
scenario. This translates into slightly different effective temperatures of
$1122\pm 16$ K and $1171\pm 17$ K, respectively, as well as surface gravities
of $4.275\pm 0.006$ dex and $4.45\pm 0.03$ dex.

The bimodality in VHS J1256$-$1257 b’s mass drives the properties we derived,
so alternative model assumptions could potentially shift the balance
significantly in favor of one or the other possibilities. The Saumon & Marley
(2008) models, in particular, have such a wide gap in the $L_{\rm bol}$–age
diagram because the onset of the L/T transition, which these models assume
begins at $\mbox{$T_{\rm eff}$}=1400$ K, happens to occur at nearly the same
age as deuterium fusion for objects near the deuterium-fusion mass boundary.
The L/T transition slows cooling, so objects stay luminous both because of
cloud disappearance and deuterium fusion. The L/T transition probably occurs
at lower temperatures for low-gravity objects like VHS J1256$-$1257 b (e.g.,
Metchev & Hillenbrand, 2006; Dupuy et al., 2009; Filippazzo et al., 2015),
which may significantly impact the size and shape of the deuterium-fusing gap
in $L_{\rm bol}$-age space. Naively, such a delayed and lower-$T_{\rm eff}$
L/T transition might be expected to make even lower-mass isochrones have
higher luminosities in Figure 3, which would in turn make it more likely that
VHS J1256$-$1257 b is indeed below the deuterium-fusion mass boundary.

##  6 Orbital architecture & Origins

We have astrometrically determined the three-dimensional orbits of both the
inner host binary (A–B) and the outer companion about its barycenter (AB–b).
This allows us to constrain the orbital architecture of the system and thus,
potentially, shed light on its origin. One crucial measurement that our orbit
determinations enables is the true mutual inclination of the A–B and AB–b
orbital planes,

$\hfill\cos{i_{\rm AB-b}}=\cos{i_{\rm AB}}\cos{i_{\rm b}}+\sin{i_{\rm
AB}}\sin{i_{\rm b}}\cos(\Omega_{\rm AB}-\Omega_{\rm b}).\hfill$ (3)

Propagating all measurement uncertainties from orvara and lofti analyses, we
find $i_{\rm AB-b}=116\pm 16$°. This reveals that the angular momentum vectors
of the two orbital planes are misaligned (7$\sigma$) and also possibly
pointing in opposite directions (1.6$\sigma$).

VHS J1256$-$1257 AB’s orbit is highly eccentric, and one possible explanation
for this is that it has been pumped up by Kozai-Lidov cycles (Kozai, 1962;
Lidov, 1962). The observed mutual inclination is consistent with the range of
critical inclinations for this mechanism to operate, $39.2\degree<i_{\rm
AB-b}<140.8\degree$. The masses, eccentricities, and orbital periods also
imply Kozai-Lidov oscillation periods less than the age of the system,
$\log(\tau_{\rm KL}/{\rm Myr})=1.6^{+0.4}_{-0.5}$ dex (Eq. 1; Fabrycky &
Tremaine, 2007). Under the conservative assumption that VHS J1256$-$1257 AB’s
initial eccentricity was zero, its maximum eccentricity attainable from Kozai-
Lidov cycles is $[1-(5/3)\cos^{2}{i_{\rm AB-b,initial}}]^{-1/2}$ (Fabrycky &
Tremaine, 2007). To achieve the observed $e_{\rm AB}=0.883$ would thus have
required an initial misalignment of $>$68.7°or $<$111.3°, the latter of which
is in excellent agreement with our measured mutual inclination.

Whatever might have caused an initial misalignment between the two orbital
planes may also be responsible for the unusual configuration of this system.
The companion mass ratio relative to the inner binary is quite low ($M_{\rm
b}/M_{\rm AB}=0.08$–0.10), especially for a multiple system with total mass
$<$0.2 $M_{\sun}$ (e.g., Burgasser et al., 2007). At higher masses, the
formation of such systems has been suggested to be due to the disintegration
of high-order multiples at young ages (e.g., Stamatellos & Whitworth, 2009;
Reipurth & Mikkola, 2015), although such systems should be quite rare (e.g.,
Bate, 2012).

One final clue to the origins of the VHS J1256$-$1257 system is the
eccentricity of the wide companion’s orbit. Its periastron distance of
$90^{+20}_{-30}$ au is consistent with a more compact initial configuration
for the system that led to VHS J1256$-$1257 b being scattered onto a wide,
eccentric, and misaligned orbit.

##  7 Summary

We measured high-precision relative astrometry of all three components of the
VHS J1256$-$1257 system. Our orbital analysis yields a total dynamical mass of
the inner binary ($M_{\rm A}+M_{\rm B}=0.141\pm 0.008$ $M_{\sun}$), high
eccentricities for both the inner and outer orbits, and a mutual inclination
of $116\pm 16$° between them. We thus confirmed that the host binary may be a
pair of brown dwarfs, derived their integrated-light luminosity from SED-
fitting, and measured a cooling age of $140\pm 20$ Myr from Saumon & Marley
(2008) hybrid evolutionary models. We found that VHS J1256$-$1257 b has a
sufficiently low luminosity that, at such a young age, it is either below the
deuterium-fusing mass boundary or, marginally less likely, massive enough to
have depleted its deuterium long ago. Regardless of the mass of VHS
J1256$-$1257 b, the orbital architecture implies a dynamical origin, perhaps
from the disintegration of a high-order multiple or scattering within a
protostellar disk.

If VHS J1256$-$1257 b is indeed below the D-fusion mass boundary, then
molecular absorption bands from D-bearing isotopologues of water (HDO) and
methane (CH3D) may be detectable in high-S/N 3–5 µm JWST spectra (e.g., Morley
et al., 2019). We also anticipate that similar observations will be possible
for the other rare triple systems with planetary-mass companions for which
substellar cooling ages are possible (e.g., 2MASS J0249$-$0557; Dupuy et al.,
2018).

## Acknowledgements

This research was funded in part by the Gordon and Betty Moore Foundation
through grant GBMF8550 to M. Liu. A. Sanghi acknowledges support from the
Research Experience for Undergraduate program at the Institute for Astronomy,
University of Hawaii, Manoa funded through NSF grant #2050710. We thank
Spencer Hurt for the BT-Settl models used in the bolometric luminosity
calculation. The data presented herein were obtained at the W.M. Keck
Observatory, which is operated as a partnership between the California
Institute of Technology, the University of California, and NASA. The
Observatory was made possible by the generous financial support of the W.M.
Keck Foundation. This work has made use of data from the European Space Agency
(ESA) mission Gaia, processed by the Gaia Data Processing and Analysis
Consortium (DPAC). Funding for the DPAC has been provided by national
institutions, in particular the institutions participating in the Gaia
Multilateral Agreement. For the purpose of open access, the author has applied
a Creative Commons Attribution (CC BY) licence to any Author Accepted
Manuscript version arising from this submission. Finally, the authors wish to
recognize and acknowledge the very significant cultural role and reverence
that the summit of Maunakea has always had within the indigenous Hawaiian
community. We are most fortunate to have the opportunity to conduct
observations from this mountain.

## Data Availability

All of our NIRC2 data are available on the Keck Observatory Archive (KOA),
which is operated by the W. M. Keck Observatory and the NASA Exoplanet Science
Institute (NExScI), under contract with the National Aeronautics and Space
Administration. We include configuration files for our orbit analysis in the
supplemental data.

## References

  * Allard et al. (2012) Allard F., Homeier D., Freytag B., 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 2765
  * Baraffe et al. (2015) Baraffe I., Homeier D., Allard F., Chabrier G., 2015, A&A, 577, A42
  * Bardalez Gagliuffi et al. (2014) Bardalez Gagliuffi D. C., et al., 2014, ApJ, 794, 143
  * Bate (2012) Bate M. R., 2012, MNRAS, 419, 3115
  * Beatty et al. (2018) Beatty T. G., Morley C. V., Curtis J. L., Burrows A., Davenport J. R. A., Montet B. T., 2018, AJ, 156, 168
  * Blunt et al. (2017) Blunt S., et al., 2017, AJ, 153, 229
  * Bond et al. (2020) Bond C. Z., et al., 2020, Journal of Astronomical Telescopes, Instruments, and Systems, 6, 039003
  * Brandt et al. (2021a) Brandt G. M., et al., 2021a, arXiv e-prints, p. arXiv:2109.07525
  * Brandt et al. (2021b) Brandt G. M., Brandt T. D., Dupuy T. J., Li Y., Michalik D., 2021b, AJ, 161, 179
  * Brandt et al. (2021c) Brandt T. D., Dupuy T. J., Li Y., Brandt G. M., Zeng Y., Michalik D., Bardalez Gagliuffi D. C., Raposo-Pulido V., 2021c, AJ, 162, 186
  * Brandt et al. (2021d) Brandt G. M., Brandt T. D., Dupuy T. J., Michalik D., Marleau G.-D., 2021d, ApJ, 915, L16
  * Burgasser (2009) Burgasser A. J., 2009, in Mamajek E. E., Soderblom D. R., Wyse R. F. G., eds, International Astronomical Union Symposium Series Vol. 258, The Ages of Stars. pp 317–326 (arXiv:0811.2508), doi:10.1017/S1743921309031974
  * Burgasser et al. (2007) Burgasser A. J., Reid I. N., Siegler N., Close L., Allen P., Lowrance P., Gizis J., 2007, in Reipurth B., Jewitt D., Keil K., eds, Protostars and Planets V. p. 427 (arXiv:astro-ph/0602122)
  * Burgasser et al. (2010) Burgasser A. J., Cruz K. L., Cushing M., Gelino C. R., Looper D. L., Faherty J. K., Kirkpatrick J. D., Reid I. N., 2010, ApJ, 710, 1142
  * Chambers et al. (2016) Chambers K. C., et al., 2016, arXiv e-prints, p. arXiv:1612.05560
  * Cutri et al. (2003) Cutri R. M., et al., 2003, 2MASS All Sky Catalog of point sources., doi:10.26131/IRSA2
  * Cutri et al. (2013) Cutri R. M., et al., 2013, Explanatory Supplement to the AllWISE Data Release Products, Explanatory Supplement to the AllWISE Data Release Products, by R. M. Cutri et al.
  * Diolaiti et al. (2000) Diolaiti E., Bendinelli O., Bonaccini D., Close L., Currie D., Parmeggiani G., 2000, A&AS, 147, 335
  * Dupuy & Liu (2011) Dupuy T. J., Liu M. C., 2011, ApJ, 733, 122
  * Dupuy & Liu (2017) Dupuy T. J., Liu M. C., 2017, ApJS, 231, 15
  * Dupuy et al. (2009) Dupuy T. J., Liu M. C., Ireland M. J., 2009, ApJ, 699, 168
  * Dupuy et al. (2014) Dupuy T. J., Liu M. C., Ireland M. J., 2014, ApJ, 790, 133
  * Dupuy et al. (2018) Dupuy T. J., et al., 2018, AJ, 156, 57
  * Dupuy et al. (2019) Dupuy T. J., Brandt T. D., Kratter K. M., Bowler B. P., 2019, ApJ, 871, L4
  * Dupuy et al. (2020) Dupuy T. J., et al., 2020, Research Notes of the American Astronomical Society, 4, 54
  * Dupuy et al. (2022) Dupuy T. J., Brandt G. M., Brandt T. D., 2022, MNRAS, 509, 4411
  * Eisenhardt et al. (2020) Eisenhardt P. R. M., et al., 2020, ApJS, 247, 69
  * Fabrycky & Tremaine (2007) Fabrycky D., Tremaine S., 2007, ApJ, 669, 1298
  * Filippazzo et al. (2015) Filippazzo J. C., Rice E. L., Faherty J., Cruz K. L., Van Gordon M. M., Looper D. L., 2015, ApJ, 810, 158
  * Foreman-Mackey (2016) Foreman-Mackey D., 2016, The Journal of Open Source Software, 1, 24
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2021) Gaia Collaboration et al., 2021, A&A, 649, A1
  * Gauza et al. (2015) Gauza B., Béjar V. J. S., Pérez-Garrido A., Zapatero Osorio M. R., Lodieu N., Rebolo R., Pallé E., Nowak G., 2015, ApJ, 804, 96
  * Goodman & Weare (2010) Goodman J., Weare J., 2010, Communications in Applied Mathematics and Computational Science, 5, 65
  * Hinkley et al. (2022) Hinkley S., et al., 2022, arXiv e-prints, p. arXiv:2205.12972
  * Hoch et al. (2022) Hoch K. K. W., et al., 2022, arXiv e-prints, p. arXiv:2207.03819
  * Ireland et al. (2008) Ireland M. J., Kraus A., Martinache F., Lloyd J. P., Tuthill P. G., 2008, ApJ, 678, 463
  * Kozai (1962) Kozai Y., 1962, AJ, 67, 591
  * Lidov (1962) Lidov M. L., 1962, Planet. Space Sci., 9, 719
  * Liu et al. (2008) Liu M. C., Dupuy T. J., Ireland M. J., 2008, ApJ, 689, 436
  * Luhman et al. (2011) Luhman K. L., Burgasser A. J., Bochanski J. J., 2011, ApJ, 730, L9
  * Markwardt (2009) Markwardt C. B., 2009, in Bohlender D. A., Durand D., Dowler P., eds, Astronomical Society of the Pacific Conference Series Vol. 411, Astronomical Data Analysis Software and Systems XVIII. p. 251 (arXiv:0902.2850)
  * Marocco et al. (2021) Marocco F., et al., 2021, ApJS, 253, 8
  * Metchev & Hillenbrand (2006) Metchev S. A., Hillenbrand L. A., 2006, ApJ, 651, 1166
  * Morley et al. (2019) Morley C. V., Skemer A. J., Miles B. E., Line M. R., Lopez E. D., Brogi M., Freedman R. S., Marley M. S., 2019, ApJ, 882, L29
  * Pearce et al. (2020) Pearce L. A., Kraus A. L., Dupuy T. J., Mann A. W., Newton E. R., Tofflemire B. M., Vanderburg A., 2020, ApJ, 894, 115
  * Reipurth & Mikkola (2015) Reipurth B., Mikkola S., 2015, AJ, 149, 145
  * Rich et al. (2016) Rich E. A., Currie T., Wisniewski J. P., Hashimoto J., Brandt T. D., Carson J. C., Kuzuhara M., Uyama T., 2016, ApJ, 830, 114
  * Saumon & Marley (2008) Saumon D., Marley M. S., 2008, ApJ, 689, 1327
  * Service et al. (2016) Service M., Lu J. R., Campbell R., Sitarski B. N., Ghez A. M., Anderson J., 2016, PASP, 128, 095004
  * Snellen & Brown (2018) Snellen I. A. G., Brown A. G. A., 2018, Nature Astronomy, 2, 883
  * Stamatellos & Whitworth (2009) Stamatellos D., Whitworth A. P., 2009, MNRAS, 392, 413
  * Stone et al. (2016) Stone J. M., et al., 2016, ApJ, 818, L12
  * Tokunaga et al. (2002) Tokunaga A. T., Simons D. A., Vacca W. D., 2002, PASP, 114, 180
  * Vousden et al. (2016) Vousden W. D., Farr W. M., Mandel I., 2016, MNRAS, 455, 1919
  * Wizinowich et al. (2006) Wizinowich P. L., et al., 2006, PASP, 118, 297


# Macroscopic Superposition States in Isolated Quantum Systems

Roman V. Buniy<EMAIL_ADDRESS>Schmid College of Science, Chapman
University, Orange, CA 92866  Stephen D. H. Hsu<EMAIL_ADDRESS>Department of
Physics and Astronomy, Michigan State University

###### Abstract

For any choice of initial state and weak assumptions about the Hamiltonian,
large isolated quantum systems undergoing Schrödinger evolution spend most of
their time in macroscopic superposition states. The result follows from von
Neumann’s 1929 Quantum Ergodic Theorem. As a specific example, we consider a
box containing a solid ball and some gas molecules. Regardless of the initial
state, the system will evolve into a quantum superposition of states with the
ball in macroscopically different positions. Thus, despite their seeming
fragility, macroscopic superposition states are ubiquitous consequences of
quantum evolution. We discuss the connection to many worlds quantum mechanics.

Highly complex superposition states have been realized in the laboratory Nobel
. Despite their seeming fragility, such states are of great importance in
quantum information and computing, as well as to theoretical questions in
quantum foundations. It may come as a surprise that isolated systems with many
degrees of freedom naturally evolve into macroscopic superposition states.
Such states contain orthogonal components which differ in macroscopic
quantities such as the position or momentum of large objects normally thought
to be “classical” in nature.

In what follows we use a specific example (essentially, that of Brownian
motion) to illustrate how this result follows from the 1929 Quantum Ergodic
Theorem (QET) of John von Neumann QET . This theorem was largely forgotten for
over 50 years, until its resurrection in 2009-10 GLMTZ ; GLTZ . The QET
contains insights relevant both to quantum statistical mechanics and the
foundations of quantum mechanics. We conclude with some remarks on the latter
topic.

The QET goes beyond typicality (concentration of measure) results concerning
thermalization in isolated quantum systems typical . The typicality result
says that almost all pure states $\Psi$ of large systems are maximally
entangled, and tracing over all but a small subspace $1$ yields a density
matrix $\rho_{1}$ which is close to the normalized identity, i.e. the
microcanonical density matrix. This implies thermal properties for the small
subspace. The QET is focused specifically on the subspace of macroscopic
observables, as opposed to a generic subset of microscopic degrees of freedom.
Von Neumann proves a stronger result concerning time evolution (ergodicity) of
the system: that all initial states $\Psi_{0}$ will spend most of their time
evolution as typical states with respect to the macro subspace (see equation
(11) below), subject of course to certain assumptions required for the theorem
RS .

The calculation given below illustrates that for any subspace of a large
system (including, for example, a subspace defined by a set of macroscopic
observables), the density operator $\rho_{1}$ induced by tracing over the
other degrees of freedom of a random pure state is overwhelmingly likely to be
close to $\rho_{1}\approx n_{1}^{-1}I_{1}$. Based on this dominance of the
measure, one can argue heuristically that even if the system begins in an
exceptional state where this property is strongly violated, dynamical
evolution will cause it to spend most of its time in a typical state. The QET
provides a rigorous basis for this intuition.

Let $\\{\phi_{1,j_{1}}\\}_{j_{1}=1}^{n_{1}}$ and
$\\{\phi_{2,j_{2}}\\}_{j_{2}=1}^{n_{2}}$ be complete sets of orthonormal
eigenstates of two Hilbert spaces labeled 1 and 2, with identity operators
$I_{1}$ and $I_{2}$. For a state

$\displaystyle\psi=\sum_{j_{1}=1}^{n_{1}}\sum_{j_{2}=1}^{n_{2}}c_{j_{1},j_{2}}\phi_{1,j_{1}}\otimes\phi_{2,j_{2}}$
(1)

the corresponding density operator $\rho=\psi\psi^{\dagger}$ and the partial
density operator $\rho_{1}=\operatorname{tr}_{2}{\rho}$ are

$\displaystyle\rho$
$\displaystyle=\sum_{j_{1},k_{1}=1}^{n_{1}}\sum_{j_{2},k_{2}=1}^{n_{2}}c_{j_{1},j_{2}}c^{*}_{k_{1},k_{2}}\phi_{1,j_{1}}\phi_{1,k_{1}}^{\dagger}\otimes\phi_{2,j_{2}}\phi_{2,k_{2}}^{\dagger},$
(2) $\displaystyle\rho_{1}$
$\displaystyle=\sum_{j_{1},k_{1}=1}^{n_{1}}\sum_{j_{2}=1}^{n_{2}}c_{j_{1},j_{2}}c^{*}_{k_{1},j_{2}}\phi_{1,j_{1}}\phi_{1,k_{1}}^{\dagger}$
(3)

We assume that $\sum_{j_{1}=1}^{n_{1}}\sum_{j_{2}=1}^{n_{2}}\lvert
c_{j_{1},j_{2}}\rvert^{2}=1$, in which case $\lVert\psi\rVert^{2}=1$ and
$\operatorname{tr}{\rho}=\operatorname{tr}_{1}{\rho_{1}}=1$. For the uniform
distribution on $S^{2n_{1}n_{2}-1}$ we have the expected value

$\displaystyle\mathsf{E}{(c_{j_{1},j_{2}}c^{*}_{k_{1},k_{2}})}$
$\displaystyle=(n_{1}n_{2})^{-1}\delta_{j_{1},k_{1}}\delta_{j_{2},k_{2}},$ (4)

which implies $\mathsf{E}{\rho_{1}}=n_{1}^{-1}I_{1}$.

We can show that $\rho_{1}$ is strongly concentrated near its expected value,
with small probability of deviation, $\mathsf{P}{}$. We use Levy’s lemma Levy
; Milman ; Jenkinson , which implies that for any Lipschitz continuous
function $f\colon\mathbb{S}^{2n_{1}n_{2}-1}\to\mathbb{R}$, a uniformly
distributed random variable $c\in\mathbb{S}^{2n_{1}n_{2}-1}$ and any
$\epsilon>\delta$,

$\displaystyle\mathsf{P}{\bigl{(}\lvert
f(c)-\mathsf{E}{f(c)}\rvert\geq\epsilon\bigr{)}}\leq\exp{\biggl{(}-\frac{n_{1}n_{2}(\epsilon-\delta)^{2}}{\lVert
f\rVert_{\textrm{L}}^{2}}\biggr{)}},$ (5)

where

$\displaystyle\lVert
f\rVert_{\textrm{L}}=\sup_{c\in\mathbb{S}^{2n_{1}n_{2}-1}}\lVert\nabla_{c}f(c)\rVert$
(6)

is the Lipschitz norm of $f$ and $\delta=(\frac{\pi}{4n_{1}n_{2}})^{1/2}$. The
specific bound (5) is deduced using Proposition 1.8 in Levy . (See Appendix
for additional details.) The result is a consequence of the extreme
concentration of measure exhibited in the geometry of high dimensional spheres
Levy ; Milman . Given that $\rho_{1}$ is quadratic in $c$, which is itself
bounded in norm, we can bound the gradient $\nabla_{c}\rho_{1}$ to obtain an
exponential concentration in probability:

$\displaystyle\mathsf{P}{\bigl{(}\lvert\operatorname{Re}{(\rho_{1})_{j_{1},j_{1}}}-n_{1}^{-1}\rvert\geq\epsilon\bigr{)}}\leq\exp{\bigl{(}-\tfrac{n_{1}n_{2}}{4}(\epsilon-\delta)^{2}}\bigr{)},$
(7)
$\displaystyle\mathsf{P}{\bigl{(}\lvert\operatorname{Re}{(\rho_{1})_{j_{1},k_{1}}}\rvert\geq\epsilon\bigr{)}}\leq\exp{\bigl{(}-n_{1}n_{2}\epsilon^{2}}\bigr{)},\quad
j_{1}\not=k_{1},$ (8)
$\displaystyle\mathsf{P}{\bigl{(}\lvert\operatorname{Im}{(\rho_{1})_{j_{1},k_{1}}}\rvert\geq\epsilon\bigr{)}}\leq\exp{\bigl{(}-n_{1}n_{2}\epsilon^{2}\bigr{)}},\quad
j_{1}\not=k_{1}.$ (9)

In particular, as $n_{1}n_{2}$ becomes large the density matrix $\rho_{1}$
becomes exponentially concentrated near $\rho_{1}=n_{1}^{-1}I_{1}$. The set of
states $\psi$ for which $\rho_{1}$ deviates from this value approaches measure
zero. (Note neither $n_{1}$ nor $n_{2}$ has to be individually large, as long
as their product is.) Thus it is plausible, and the QET demonstrates, that
under dynamical evolution of $\psi$: $\rho_{1}\approx n_{1}^{-1}I_{1}$ almost
all of the time.

Now let us turn to a specific example, and focus on macroscopic superposition
states. Consider a solid ball and gas molecules isolated in a box. The
Hamiltonian can be divided into three terms: $H_{\textrm{B}}$ describing the
ball, $H_{\textrm{G}}$ for the gas molecules, and $H_{\textrm{BG}}$ describing
interactions between the two. Let

$H_{\textrm{B}}=\frac{P^{2}}{2M}+V(Q)+\sum_{i=1}^{N}\frac{p_{i}^{2}}{2m}+U(q_{1},\dotsc,q_{N}),$
(10)

where $P$ and $Q$ are center of mass momentum and position, and $i$ label
relative coordinates and momenta Messiah  of the $N+1$ individual molecules
making up the ball. Importantly, $P$ and $Q$ commute with the relative degrees
of freedom.

We take the potential $V$ to be non-zero, so that the energy of the ball
varies very slightly with location in the box (i.e., the bottom of the box is
not exactly flat). This assumption is motivated by technical aspects of the
proof of the QET (specifically the assumption that energy levels and gaps are
non-degenerate) but does not affect the main physical aspects of our
discussion. In a realistic setting we expect that the technical assumptions
used in the QET to be satisfied. The details of $H_{\textrm{G}}$ and
$H_{\textrm{BG}}$ are not important, except that we assume that the gas
molecules are excluded from the space occupied by the ball. For example,
$H_{\textrm{BG}}$ may contain repulsive potential terms which depend on $Q$.

The macroscopic quantities we are interested in are the total energy of the
system and the center of mass position $Q$ and momentum $P$ of the ball. The
macroscopic energy operator is a coarse grained version of the (microscopic)
Hamiltonian $H=H_{\textrm{B}}+H_{\textrm{G}}+H_{\textrm{BG}}$. The eigenvalues
of $H$ are grouped into bands (energy shells) whose width is macroscopic
(i.e., there are many eigenvalues in each band; the energy difference between
bands is experimentally accessible). Each of the $H$ eigenstates in a band is
assigned a coarse grained eigenvalue equal to the average energy in the band,
so that the macroscopic energy operator has very large degeneracies (all of
the states in a band share the same eigenvalue).

The definition of the macroscopic energy given above can be generalized: given
a microscopic operator $A$, which has eigenstate expansion
$A=\sum_{a}a|a\rangle\langle a|$, the coarse grained or macroscopic version is
obtained by dividing the sum into bands of nearby eigenvalues, and replacing
the coefficients $a$ with the average within the band.

An important property of macroscopic observables is that they can be measured
simultaneously, so the corresponding operators must commute. The macroscopic
momentum $P$ and position $Q$ of the ball are constructed using wave packet
states which have simultaneously a central value of momentum and of position,
and a spread in each which is small but respects the uncertainty principle.
These states are orthogonalized, resulting in simultaneous eigenstates for
macroscopic $P$ and $Q$. The procedure is discussed in von Neumann’s proof of
the QET.

The full set of mutually commuting macroscopic operators $\\{M_{i}\\}$
includes coarse grained versions of $H,P,Q$. Within a specific energy shell of
fixed macroscopic energy there are many individual phase space cells with
different values of macroscopic $(P,Q)$. In the general case other quantities
may be included in $\\{M_{i}\\}$.

We restrict attention to a specific energy shell of dimension $D$, which is
further partitioned into subspaces with definite values of each of these
macroscopic variables $\\{M_{i}\\}$. Let $\nu$ denote a full set of values for
the macroscopic variables, and $d_{\nu}$ denote the dimensionality of the
subspace with values labeled by $\nu$, so $D=\sum_{\nu}d_{\nu}$.

The Quantum Ergodic Theorem QET ; GLMTZ ; GLTZ  states that: for any initial
state $\Psi_{0}$ and almost any Hamiltonian $H$, in the long run the system
spends almost all of its time in typical states $\Psi$ with the property that

$\displaystyle\langle\Psi\lvert\bm{P}_{\nu}\rvert\Psi\rangle\approx\frac{d_{\nu}}{D},$
(11)

where $\bm{P}_{\nu}$ is a projector onto the macro subspace labeled by $\nu$.
The conditions on $H$ are the absence of degeneracies (which are unlikely, in
the absence of exact symmetries) and a technical condition on the relation
between the energy eigenstate basis and the $\nu$ basis. In the limit of large
number of degrees of freedom the technical condition is violated only for a
vanishingly small subset of macroscopic operators Reimann . Both conditions
can be made to hold in realistic settings, such as the case of Brownian motion
of a macroscopic ball interacting with a gas, for natural choices of $H$. In
particular, if a specific setup does not satisfy the conditions (e.g., choice
of interaction terms, box geometry, or specific coarse grained macroscopic
operator; this would be highly unlikely for the reasons given), a slightly
perturbed version of the setup is almost certain to suffice. The theorem has
been explicitly tested in simulations of specific systems RS .

The QET implies that in the typical state (i.e., most of the time, in the long
run) there is significant probability to find the ball in any of the phase
space cells labeled by $\nu$, as long as the $d_{\nu}$ are similar in
magnitude. In our example this would mean that the ball position could be
roughly anywhere in the box.

We can show that this corresponds to a macroscopic superposition state by
noting that $\Psi$ must be a pure state at all times, which by the Schmidt
decomposition must have the form

$\displaystyle\Psi=\sum_{\nu}c_{\nu}\psi_{\nu}\otimes\phi_{\nu}$ (12)

where $\psi_{\nu}$ describes the (macroscopic) state of the ball and
$\phi_{\nu}$ describes the gas molecules. By construction, the $\nu$ states
are mutually orthogonal: i.e., $\psi_{\nu}$ has zero overlap with
$\psi_{\nu^{\prime}}$ for $\nu\neq\nu^{\prime}$ and similarly with the $\phi$
states. The gas molecule (“environment”) state is entangled with the ball
state, as the gas molecules cannot occupy the same space as the ball. The QET
implies that the magnitudes $\lvert c_{\nu}\rvert$ are comparable as long as
the corresponding $d_{\nu}$ are similar in size. In concrete terms this
assumption is satisfied for $\nu$ corresponding to phase space cells $(P,Q)$
in which the macroscopic energies of the ball are nearly the same, leaving
similar volumes of phase space to be occupied by the gas molecules (subject to
the total energy constraint for the given macroscopic energy shell).

Note the claim is not that the ball spends some of the time at $\nu$ and some
at $\nu^{\prime}$. Our focus is not the time average of its position. The
claim is stronger: that at almost all times the ball is in a macroscopic
superposition state, and the entangled environment (air molecules, or even an
observer whose brain is a macroscopic neural net) is as well. This all follows
from the assumption of unitary (Schrödinger) evolution inside the box, using
the QET. The basis denoted by $\nu$ is the natural basis for a macroscopic
observer of the system.

We can give the following intuitive (quasi-classical) explanation for the
result. The ball is subject to fluctuations arising from collisions with the
gas — i.e., Brownian motion. Over long enough timescales there is a non-
negligible amplitude for fluctuations to cause the ball to random walk from
its starting position to any other position in the box. As Schrödinger
dynamics is linear, $\Psi$ evolves into a superposition of the possible
outcome states. The location of the ball is reflected in the corresponding
state of the gas molecules (ball and gas cannot occupy the same space), which
act as an environment for the ball. After enough time, the total wavefunction
can be decomposed into orthogonal branches which differ macroscopically. This
description suggests that the timescale over which an initially localized ball
state evolves into a typical $\Psi$ satisfying equation (11) is similar to
that for a Brownian random walk to reach all parts of the box.

The example above can be generalized to almost any system whose dynamics
permit multiple macroscopic outcomes. The QET implies that the quantum state
evolves into a superposition over these outcomes. The ubiquity of macroscopic
superposition states under ordinary Schrödinger evolution is of course an
aspect of many worlds, or no collapse, quantum mechanics manyworlds . This
view of quantum mechanics is typically introduced via the measurement process,
which we briefly summarize below.

Let $Q$ be a single qubit and $M$ a macroscopic device which measures the spin
of the qubit along a particular axis. The eigenstates of spin along this axis
are denoted $|\pm\rangle$. We define the operation of $M$ as follows, where
the combined system is $S=Q+M$:

$\displaystyle|+\rangle\otimes|M\rangle\longrightarrow|S_{+}\rangle,$ (13)
$\displaystyle|-\rangle\otimes|M\rangle\longrightarrow|S_{-}\rangle,$ (14)

where $S_{+}$ denotes a state of the total system in which the apparatus $M$
has recorded a $+$ outcome (i.e., is in state $M_{+}$), and similarly with
$S_{-}$. We can then ask what happens to a superposition state
$|\Psi_{Q}\rangle=c_{+}|+\rangle+c_{-}|-\rangle$ which enters the device $M$.
In the conventional formulation, with measurement collapse, one of the two
final states $|S_{+}\rangle$ or $\,|S_{-}\rangle$ is realized, with
probabilities $|c_{+}|^{2}$ and $|c_{-}|^{2}$ respectively.

However, if the combined system $S=Q+M$ evolves according to the Schrödinger
equation (in particular, linearly), we obtain a superposition of measurement
device states:

$\displaystyle\big{(}c_{+}|+\rangle~{}+~{}c_{-}|-\rangle\big{)}\otimes|M\rangle~{}~{}\longrightarrow~{}~{}c_{+}\,|S_{+}\rangle~{}+~{}c_{-}\,|S_{-}\rangle.$
(15)

This seems counter to actual experience: measurements produce a single
outcome, not a superposition state. But an observer in the state $S_{+}$ might
be unaware of the second branch of the wave function in state $S_{-}$ for
dynamical reasons related to a phenomenon called decoherence decoherence . Any
object sufficiently complex to be considered either a measuring device or
observer (i.e., which, in the conventional formulation can be regarded as
semi-classical) will have many degrees of freedom. A measurement can only be
said to have occurred if the states $M_{+}$ and $M_{-}$ are very different:
the outcome of the measurement must be stored in a redundant and
macroscopically accessible way in the device (or, equivalently, in the local
environment). Therefore, the overlap of $M_{+}$ with $M_{-}$ is of order
$\exp(-N)$, where $N$ is a macroscopic number of degrees of freedom, and the
future dynamical evolution of each branch is unlikely to alter this. For all
practical purposes, as John Bell put it Bell , an observer on one branch can
ignore the existence of the other: they are said to have decohered. Each of
the two observers will perceive a collapse to have occurred, although the
evolution of the overall system $S$ has continued to obey the Schrödinger
equation.

The Quantum Ergodic Theorem shows that macroscopic superpositions of outcome
states like $|S_{\pm}\rangle$ are ubiquitous under Schrödinger evolution.
Interactions such as the scattering of particles in Brownian motion lead to
such states, even if observers in the closed system (described by the total
wave function $\Psi$) are unaware of it.

## Appendix A

Writing
$c_{j_{1},j_{2}}=c^{\prime}_{j_{1},j_{2}}+ic^{\prime\prime}_{j_{1},j_{2}}$,
where $c^{\prime}_{j_{1},k_{1}},c^{\prime\prime}_{j_{1},k_{1}}\in\mathbb{R}$,

$\displaystyle\sum_{j_{1}=1}^{n_{1}}\sum_{j_{2}=1}^{n_{2}}(c^{\prime
2}_{1,j_{1},k_{1}}+c^{\prime\prime 2}_{j_{1},k_{1}})=1,$ (16)

and using

$\displaystyle\operatorname{Re}{(\rho_{1})_{j_{1},k_{1}}}=\sum_{j_{2}=1}^{n_{2}}(c^{\prime}_{j_{1},j_{2}}c^{\prime}_{k_{1},j_{2}}+c^{\prime\prime}_{j_{1},j_{2}}c^{\prime\prime}_{k_{1},j_{2}}),$
(17)
$\displaystyle\operatorname{Im}{(\rho_{1})_{j_{1},k_{1}}}=\sum_{j_{2}=1}^{n_{2}}(c^{\prime\prime}_{j_{1},j_{2}}c^{\prime}_{k_{1},j_{2}}-c^{\prime}_{j_{1},j_{2}}c^{\prime\prime}_{k_{1},j_{2}}),$
(18)

we find

$\displaystyle\lVert\nabla_{c}\operatorname{Re}{(\rho_{1})_{j_{1},k_{1}}}\rVert^{2}$
$\displaystyle=\lVert\bigl{(}\delta_{l_{1},j_{1}}c^{\prime}_{k_{1},l_{2}}+c^{\prime}_{j_{1},l_{2}}\delta_{l_{1},k_{1}},\delta_{l_{1},j_{1}}c^{\prime\prime}_{k_{1},l_{2}}+c^{\prime\prime}_{j_{1},l_{2}}\delta_{l_{1},k_{1}}\bigr{)}_{1\leq
l_{1}\leq n_{1},1\leq l_{2}\leq n_{2}}\rVert^{2}$
$\displaystyle=\sum_{l_{2}=1}^{n_{2}}\bigl{(}c^{\prime
2}_{j_{1},l_{2}}+c^{\prime\prime 2}_{j_{1},l_{2}}+c^{\prime
2}_{k_{1},l_{2}}+c^{\prime\prime
2}_{k_{1},l_{2}}\bigr{)}+2\delta_{j_{1},k_{1}}\sum_{l_{2}=1}^{n_{2}}\bigl{(}c^{\prime
2}_{j_{1},l_{2}}+c^{\prime\prime 2}_{j_{1},l_{2}}\bigr{)}$
$\displaystyle=\sum_{l_{2}=1}^{n_{2}}\bigl{(}\lvert
c_{j_{1},l_{2}}\rvert^{2}+\lvert
c_{k_{1},l_{2}}\rvert^{2}\bigr{)}+2\delta_{j_{1},k_{1}}\sum_{l_{2}=1}^{n_{2}}\lvert
c_{j_{1},l_{2}}\rvert^{2}$ (19)
$\displaystyle\lVert\nabla_{c}\operatorname{Im}{(\rho_{1})_{j_{1},k_{1}}}\rVert^{2}$
$\displaystyle=\lVert\bigl{(}c^{\prime\prime}_{j_{1},l_{2}}\delta_{l_{1},k_{1}}-\delta_{l_{1},j_{1}}c^{\prime\prime}_{k_{1},l_{2}},\delta_{l_{1},j_{1}}c^{\prime}_{k_{1},l_{2}}-c^{\prime}_{j_{1},l_{2}}\delta_{l_{1},k_{1}}\bigr{)}_{1\leq
l_{1}\leq n_{1},1\leq l_{2}\leq n_{2}}\rVert^{2}$
$\displaystyle=\sum_{l_{2}=1}^{n_{2}}\bigl{(}c^{\prime
2}_{j_{1},l_{2}}+c^{\prime\prime 2}_{j_{1},l_{2}}+c^{\prime
2}_{k_{1},l_{2}}+c^{\prime\prime
2}_{k_{1},l_{2}}\bigr{)}-2\delta_{j_{1},k_{1}}\sum_{l_{2}=1}^{n_{2}}\bigl{(}c^{\prime
2}_{j_{1},l_{2}}+c^{\prime\prime 2}_{j_{1},l_{2}}\bigr{)}$
$\displaystyle=\sum_{l_{2}=1}^{n_{2}}\bigl{(}\lvert
c_{j_{1},l_{2}}\rvert^{2}+\lvert
c_{k_{1},l_{2}}\rvert^{2}\bigr{)}-2\delta_{j_{1},k_{1}}\sum_{l_{2}=1}^{n_{2}}\lvert
c_{j_{1},l_{2}}\rvert^{2}.$ (20)

As a result,

$\displaystyle\lVert\nabla_{c}\operatorname{Re}{(\rho_{1})_{j_{1},j_{1}}}\rVert^{2}\leq
4,$ (21)
$\displaystyle\lVert\nabla_{c}\operatorname{Re}{(\rho_{1})_{j_{1},k_{1}}}\rVert^{2}\leq
1,\quad j_{1}\not=k_{1},$ (22)
$\displaystyle\lVert\nabla_{c}\operatorname{Im}{(\rho_{1})_{j_{1},k_{1}}}\rVert^{2}\leq
1,\quad j_{1}\not=k_{1}.$ (23)

For bounds (8) and (9) on the off-diagonal terms in $\rho_{1}$ we use a
somewhat stronger inequality as the median $\mathsf{M}{f}$ and expectation
$\mathsf{E}{f}$ coincide: $\mathsf{P}{(\lvert
f-\mathsf{M}{f}\rvert\geq\epsilon)}\leq\exp(-n_{1}n_{2}\epsilon^{2}\lVert
f\rVert^{-2}_{\textrm{L}})$, Proposition A.0.5 in Jenkinson . For probability
of deviation from expectation we use Proposition 1.8 in Levy  and the median
result to obtain

$\displaystyle\mathsf{P}{(\lvert
f-\mathsf{E}{f}\rvert\geq\epsilon)}\leq\exp{\Bigl{(}-\frac{n_{1}n_{2}(\epsilon-\delta)^{2}}{\lVert
f\rVert^{2}_{\textrm{L}}}\Bigr{)}}$ (24)

with $\delta=(\frac{\pi}{4n_{1}n_{2}})^{1/2}\lVert f\rVert_{\textrm{L}}$ and
$\epsilon\geq\delta$. The specific bounds we obtain from concentration of
measure are stronger than those previously given in the literature as
expressions of Levy’s Lemma. They can be applied, for example, to relatively
small quantum systems in the laboratory.

## References

  * (1) D. Wineland, Nobel Lecture: Superposition, Entanglement, and Raising Schrödinger’s Cat, Reviews of Modern Physics 85, 1103 (2013).
  * (2) J. von Neumann, Beweis des Ergodensatzes und des H-Theorems, German original in Zeitschrift fuer Physik 57, 30 (1929); Proof of the Ergodic Theorem and the H-Theorem in Quantum Mechanics, English translation by Roderich Tumulka, European Phys. J. H 35, 201 (2010).
  * (3) S. Goldstein, J. L. Lebowitz, C. Mastrodonato, R. Tumulka, N. Zanghı, Normal Typicality and von Neumann’s Quantum Ergodic Theorem, Proceedings of the Royal Society London A 466 3203 (2010).
  * (4) S. Goldstein, J. L. Lebowitz, R. Tumulka, N. Zanghì, Long-Time Behavior of Macroscopic Quantum Systems: Commentary Accompanying the English Translation of John von Neumann’s 1929 Article on the Quantum Ergodic Theorem, European Phys. J. H 35, 173 (2010).
  * (5) S. Popescu, A. J. Short, A. Winter, Entanglement and the Foundation of Statistical Mechanics, Nature Physics 2, 754 (2006).
  * (6) M. Rigol and M. Srednicki, Alternatives to Eigenstate Thermalization, Phys. Rev. Lett. 108, 110601 (2012).
  * (7) M. Ledoux, The Concentration of Measure Phenomenon, Mathematical Surveys and Monographs 89, American Mathematical Society (2001).
  * (8) V. Milman and G. Schechtman, Asymptotic Theory of Finite Dimensional Normed Spaces, Lecture Notes in Mathematics, Springer (1986).
  * (9) J. Jenkinson, Convex Geometric Connections to Information Theory, Doctoral Dissertation, Case Western Reserve University, (2013).
  * (10) A. Messiah, Quantum Mechanics, Volume 1, North Holland Publishing Company, p.361 (1961).
  * (11) P. Reimann, Phys. Rev. Lett. 115, 010403 (2015).
  * (12) H. Everett, III, Rev. Mod. Phys. 29, 454 (1957); The Theory of the Universal Wavefunction (PhD thesis), reprinted in B. S. DeWitt and R. N. Graham (eds.), The Many-Worlds Interpretation of Quantum Mechanics, Princeton University Press (1973); J. Barrett and P. Byrne, The Everett Interpretation of Quantum Mechanics: Collected Works 1955–1980 with Commentary, Princeton University Press (2012); S. D. H. Hsu, Mod. Phys. Lett. A 27, 1230014 (2012).
  * (13) D. Giulini, C. Kiefer, E. Joos, J. Kupsch, I. O. Stamatescu, H. D. Zeh, Decoherence and the Appearance of a Classical World in Quantum Theory, Springer (2003).
  * (14) J. S. Bell, Phys. World 3, 33 (1990).


# Nano-scale displacement sensing based on Van der Waals interaction

Lin Hu  ICQD/Hefei National Laboratory for Physical Sciences at Microscale,
and Key Laboratory of Strongly-Coupled Quantum Matter Physics, Chinese Academy
of Sciences, and Department of Physics, University of Science and Technology
of China, Hefei, Anhui 230026, China  Synergetic Innovation Center of Quantum
Information $\&$ Quantum Physics, University of Science and Technology of
China, Hefei, Anhui 230026, China  Jin Zhao<EMAIL_ADDRESS>ICQD/Hefei
National Laboratory for Physical Sciences at Microscale, and Key Laboratory of
Strongly-Coupled Quantum Matter Physics, Chinese Academy of Sciences, and
Department of Physics, University of Science and Technology of China, Hefei,
Anhui 230026, China  Synergetic Innovation Center of Quantum Information $\&$
Quantum Physics, University of Science and Technology of China, Hefei, Anhui
230026, China  Jinlong Yang<EMAIL_ADDRESS>ICQD/Hefei National Laboratory
for Physical Sciences at Microscale, and Key Laboratory of Strongly-Coupled
Quantum Matter Physics, Chinese Academy of Sciences, and Department of
Physics, University of Science and Technology of China, Hefei, Anhui 230026,
China  Synergetic Innovation Center of Quantum Information $\&$ Quantum
Physics, University of Science and Technology of China, Hefei, Anhui 230026,
China

###### Abstract

We propose the nano-scale displacement sensor with high resolution for weak-
force systems could be realized based on vertical stacked two-dimensional (2D)
atomic corrugated layer materials bound through Van der Waals (VdW)
interaction. Using first-principles calculations, we found the electronic
structure of bi-layer blue phosphorus (BLBP) varies appreciably to both the
lateral and vertical interlayer displacement. The variation of electronic
structure due to the lateral displacement is attributed to the changing of the
interlayer distance $d_{z}$ led by atomic layer corrugation, which is in a
uniform picture with vertical displacement. Despite different stacking
configurations, the change of in-direct band gap is proportional to
$d_{z}^{-2}$. This stacking configuration independent $d_{z}^{-2}$ law is
found also works for other graphene-like corrugated bi-layer materials, for
example MoS2. By measuring the tunable electronic structure using absorption
spectroscopy, the nano-scale displacement could be detected. BLBP represents a
large family of bi-layer 2D atomic corrugated materials for which the
electronic structure is sensitive to the interlayer vertical and lateral
displacement, thus could be used for nano-scale displacement sensor. Since
this kind of sensor is established on atomic layers coupled through VdW
interaction, it provides unique applications in measurements of nano-scale
displacement induced by tiny external force.

###### pacs:

High-resolution displacement sensing is essential for scientific measurement
with a wide range of applications, among which, the nano-scale positioning
measurement is becoming increasingly important in ultra-precision
manufacturing and metrology in nano-technological applications. On micro-scale
level, many micro-sensors such as microphones, accelerometers and pressure
sensors, ultimately rely on accurately detecting small displacements with
micro-scale resolution 1 ; 2 ; 3 . With the development of nanotechnology, the
transducer fundamental concepts at the micro-scale are extended into the nano-
scale region 4 . Nano-scale displacement sensing with high resolution has
urgent applications for nano-devices.

Conventional displacement sensing are based on optical interferometry 5 ,
capacitance 6 , and piezoresistance 7 ; 8 . For the nano-scale sensing, the
conventional methods encounters different drawbacks. The optical methods are
limited by the drawbacks of non-linearity, fringe effects, proportionality
errors, volume and cost of optical components, and mostly important, the light
wavelength. Capacitive sensing suffers from signal loss and is limited by area
6 . Piezoresistors scaled down to nano-scale leads to strongly increased
resistance and the resistance noise followed 7 ; 8 . To achieve the high
resolution in nano-scale, new materials, for example nanotubes, have been
applied for nano-scale sensing 4 ; 9 ; 10 ; 11 ; 12 ; 13 . Yet despite these
attempts in nano-materials, it is still complex to get high resolution for
displacement induced by tiny force. For example, to measure the displacement
of a quantum mechanical oscillator, a single electron transistor was applied
as a displacement sensor to achieve the quantum-limited sensitivity 1 . Hence
the alternative complementary methods of measuring precise displacement for
weak-force system have urgent need in both scientific research and practical
applications.

Vertical stacked two-dimensional (2D) atomic layers formed by Van der Waals
(VdW) interaction arouse intense research interest recently. The weak VdW
interaction leads to many new physical phenomena including new van Hove
singularities 14 ; 15 ; 16 ; 17 , Fermi velocity renormalization 18 ; 19  and
Hofstadter s butterly pattern 20 ; 21 ; 22 ; 23 ; 24  in graphene bi-layer as
well as the tunable interlayer coupling in bi-layer semiconducting MoS2 25 .
In this Letter, we propose a new nano-scale displacement sensor based on bi-
layer blue phosphorus (BLBP). Blue phosphorus is a new two dimensional (2D)
material predicted by first principles calculations 26 . Similar with other 2D
layered materials, BLBP is coupled by VdW interaction. Using first-principles
calculations within density functional theory (DFT), we found the electronic
structure of BLBP is significantly sensitive to the lateral and vertical
displacement between the two layers. In a uniform picture with vertical
displacement, the sensitivity of electronic structure to the lateral
displacement is due to the variation of the vertical interlayer distance
($d_{z}$) led by atomic layer corrugation. In spite of different stacking
configurations, the change of in-direct band gap is proportional to
$d_{z}^{-2}$. By detecting the change of electronic structure using absorption
spectroscopy, one can measure the nano-scale displacement. This proposed
sensor is based on VdW interaction, therefore, it provides unique applications
in weak-force systems.

First-principles calculations are performed with the Vienna Ab-initio
Simulation Package (VASP) 27 ; 28 . The Perdew-Burke-Ernzerhof (GGA-PBE) 29
functional with VdW correction proposed by Grimme (DFT-D2) 30  is used to
optimize the structures and get the initial electronic structures 31 ; 32 ; 33
; 34 ; 35 ; 36 ; 37 ; 38 ; 39 ; 40 ; 41 ; 42 . An energy cutoff of 600 eV is
selected for the plane wave basis sets. The projector-augmented wave (PAW) 43
; 44  method is used for electron-ion interaction. A vacuum spacing of about
15Å ensures that the interactions between the layers are negligible. The
dipole correction is applied to compensate the dipole interaction in the
direction perpendicular to the plane. The k-point sampling uses the Monkhorst-
Pack scheme 45  and employs a $11\times 11\times 1$ mesh. The criterion of
maximum force during optimization on each atom is less than 0.01 eV/Å and the
convergence criteria for energy is $10^{-5}$ eV. We scan the size of unit cell
for each system to obtain the lattice parameters with the lowest energy. The
spin-orbital coupling effect is included in the calculations. It is well known
that DFT usually underestimated the band gap of semiconductor. To correct this
problem, we chose HSE 46  functional to obtain all the electronic structure
information.

BLBP is a typical VdW bi-layer 2D system. The binding energy calculated by
DFT-D2 is 25 meV/atom with a vertical interlayer distance as 3.22Å for an A-A
stacking bi-layer, which is at the same magnitude with graphene. The novelty
for this bi-layer semiconductor is that the band structure significantly
depends on the stacking configuration of the two BP layers. Fig.1 (b) and (c)
give plottings of the in-direct band gap ($E_{gap}$) and total energy change
($E_{total}$) when we move the top layer BP laterally into different
directions while keeping the bottom layer fixed as shown in Fig.1 (a). One can
see that the band gap changes from 1.75 eV to 2.33 eV while the total energy
changes within 12 meV. In Fig.1(c) the energies of valance band maximum (VBM)
and conduction band minimum (CBM) respected to the vacuum level are plotted.
It is shown that CBM is sensitive to the lateral displacement while VBM almost
remains unchanged.

Figure 1: (color online)The BLBP electronic structure dependence on the
lateral displacement. (a) The schematic diagram of lateral interlayer
displacement in which the bottom layer is fixed and top layer is moved
laterally. (b-d) Three dimensional plotting of $E_{gap}$ (b), $E_{total}$ (c),
$E_{VBM}$ and $E_{CBM}$ (d) when the top layer is moved laterally in $x$ and
$y$ directions. For $E_{VBM}$ and $E_{CBM}$ the energy of vacuum level setting
to zero.

We chose the $x$ direction in which $E_{gap}$ changes the most significantly
to investigate the detailed electronic structure. Fig.2 (a) shows the three
representative stacking configurations marked as A-A, A-B and A-B’ when the
upper layer moves along $x$ direction. $E_{gap}$ changes from 1.75 eV for A-A
stacking to 2.04 eV for A-B stacking and then to 2.33 eV for A-B’ stacking as
shown in Fig.2(b). Consistent with Fig.1(d), the band structure plottings in
Fig.2 (b) also show that the remarkable band gap change is mainly induced by
the change of the conduction band (CB). The valance band (VB) dispersion also
varies however the energy of VBM remains unchanged. To understand this, in
Fig.2 (c) we plot out the spatial probability density distribution of VBM and
CBM of BLBP with different stacking. It is clearly shown that for VBM, there
is almost no interlayer interaction. Yet CBM is shown to be the bonding state
between two layers which contains significant delocalized orbital probability
density in the interlayer region.

Figure 2: (color online) Atomic geometry structures (a), band structures (b)
and oribital probability density distribution of VBM and CBM (c) of A-A, A-B
and A-B’ stacking BLBP. Corresponding $E_{total}$, $E_{gap}$ and $d_{z}$ are
also given. For $E_{total}$, the energy of A-A stacking is set to zero. For
the oribital probability density distribution iso-surface plotting, the iso-
value is 0.004 $e/au^{3}$.

The $\sigma$ bonding character and delocalized probability density
distribution of CBM suggest that the vertical interlayer distance $d_{z}$ to
be a key factor to determine the energy of CBM. In Fig.3 (b) we plot out
$E_{total}$, CBM energy ($E_{CBM}$) and $d_{z}$ corresponding to lateral
displacement in $x$ direction ($d_{x}$) of the top layer. One can see that the
$E_{total}$, $E_{CBM}$ and $d_{z}$ follow the same trend correlated with
$d_{x}$. $d_{z}$ changes remarkably from 3.2 Å for A-A stacking to 4.1 Å for
A-B’ stacking as shown in Fig.2 (a). The strongly configuration dependent
interlayer distance can be understood by repulsive steric effects of the
corrugated atomic layers as discussed by Liu et.al. 25 . The lowest energy A-A
stacking has the smallest dz because the protrudent atoms in the top layer sit
in the vacancies of the bottom layer resulting in reduced repulsion. For A-B’
stacking the protrudent atoms in top and bottom layers are face-to-face,
leading to the largest repulsion and $d_{z}$.

Figure 3: (color online) (a) The change of $E_{total}$, $E_{CBM}$ and $d_{z}$
correlated with the lateral interlayer displacement $d_{x}$ on the route of
A-A C A-B C A-B’ C A-A. (b) The change of $E_{total}$ and $E_{CBM}$ correlated
with $d_{z}$ for the vertical interlayer displacement of A-A, A-B and A-B’
stacking configurations. The fitting of $E_{CBM}=E_{0}-\beta d_{z}^{-2}$ is
also shown.

The distinct variation of $E_{CBM}$ with changing of $d_{z}$ implies the
vertical displacement may affect the electronic structure. In Fig.3 (b) we
plot out $E_{CBM}$ and $E_{total}$ when we move the top layer BP vertically
for different stacking configurations. We found that despite the different
stacking configurations, $E_{CBM}$ follows quadratic curve equation which
could be fit as: $E_{CBM}=E_{0}-\beta d_{z}^{-2}$, where $E_{0}$ and $\beta$
are two parameters (The values are given in Table S1). $E_{0}$ converges to
$E_{CBM}$ for single layer BP when $d_{z}$ is infinite while $\beta$ is
determined by the interlayer orbital hopping integral of CBM. The independence
of $\beta$ on stacking configuration implys the interlayer interaction for
different stacking is similar. This $d_{z}^{-2}$ variation is supported by
tight-binding model presented in ref 47 . Because the energy of VBM keeps as a
constant, it is clear that $E_{gap}$ also follows the same $d_{z}^{-2}$ law
with $E_{CBM}$. The variation of $E_{CBM}$ and $E_{gap}$ is as large as 1.5 eV
within the vertical displacement of 1.8 Å and total energy change of 20 meV.

Figure 4: (color online). (a) The schematic design diagram of the nano-scale
displacement sensor based on BLBP. (b-c) The Absorption intensity (b) and band
structure (c) of BLBP when the top layer is moved laterally from A-A (black
line) through A-B (cyan line) to A-B’ (purple line). (d) The change of
$E_{PeakII}$ correlated with $d_{z}$ for the vertical interlayer displacement
of A-A, A-B and A-B’ stacking configurations. The fitting of
$E_{PeakII}=E_{0}-\beta d_{z}^{-2}$ is also shown.

Because the electronic structure is very sensitive to the interlayer
displacement, BLBP is proposed to be used as a nano-scale displacement sensor
for weak-force systems. We provide a schematic design as shown in Fig.4 (a).
During the measurement, the bottom layer is fixed while the top layer can be
shifted laterally or vertically by a very weak force. The displacement is
determined by measuring the electronic structure change. We propose to measure
the electronic structure using photo absorption spectroscopy. Fig.4 (b) and
(c) shows the calculated absorption intensity and the band structure change
when the bottom layer is fixed and top layer is shifted laterally from A-A
stacking to A-B’ stacking. Two absorption peaks are obtained when the
excitation energy is below 5 eV. Peak I at 2.75 eV which corresponds to the
vertical transition from VB to CBM does not change with the top layer lateral
displacement. This is because of this transition, the initial state in VB
changes in the same trend with CBM as shown in Fig.4 (c). Peak II, which
corresponds to the vertical excitation from VBM to CB, changes from 3.65 eV to
4.11 eV, showing significant sensitivity of top layer displacement because VBM
remains constant and the final state in CB varies distinctly. (See the details
of transitions I and II in supplementary materials.) Through the vertical
displacement of different stacking configurations, we found the energy of peak
II ($E_{PeakII}$) also follows the simple $d_{z}^{-2}$ trend as shown in Fig.4
(d). We propose that peak I can be used as a benchmark while the displacement
could be deduced from the position of peak II. Especially, regardless of
different stacking configurations, the vertical displacement could be easily
obtained through the simple $E_{PeakII}\propto d_{z}^{-2}$ relationship. In
addition, the BLBP electronic structure is also sensitive to the interlayer
twist angle. (see the supplementary materials) Therefore, this schematic
sensor can be upgraded to a four dimensional (three dimensions plus twist
angle) displacement sensor in nano-scale.

Figure 5: (color online).The change of $E_{gap}$ with $d_{z}$ for the vertical
interlayer displacement of A-A, A-B and A-B’ stacking configurations for BLBP,
bi-layer MoS2 and black phosphrus. For BLBP and MoS2, the fitting of
$E_{gap}=E_{0}-\beta d_{z}^{-2}$ is also shown.

We propose this principle displacement sensor could be realized using
different vertically stacked 2D bi-layer materials if they satisfy the
following two requirements: (i) The atomic layer is corrugated so that the
interlayer distance $d_{z}$ changes significantly with lateral displacement
and twisting. (ii) The CB or VB has significant delocalized interlayer
probability density distribution. Therefore the band structure of CB/VB
changes remarkably with different $d_{z}$. We notice that analogous
semiconductors as bi-layer MoS2 and black phosphorus are found to be sensitive
to the stacking configurations, which support our deduction 25 ; 48 ; 49 . To
convince this, we have calculated the $E_{gap}$ change versus $d_{z}$ of bi-
layer MoS2 and black phosphorus when the top layer is vertically displaced as
shown in Fig.5. Appreciable $E_{gap}$ change can be observed for these three
bi-layer materials. We found that $E_{gap}$ of bi-layer MoS2 and BLBP, which
are graphene-like 2D materials, follow the unified $d_{z}^{-2}$ law despite
different stacking configuration. For black phosphorus, since atomic structure
is different and the corrugation is very large (Fig. S4), the unified
$d_{z}^{-2}$ law for different stacking configurations fails. Thus we conclude
that both BLBP and bi-layer MoS2 can be easily used to measure the $d_{z}$
change since $E_{PeakII}\propto d_{z}^{-2}$, however, BLBP has larger $\beta$
value (shown in Table S1) meaning the $E_{gap}$ tunability of BLBP is larger.
Therefore we propose BLBP is a prominent candidate for displacement sensor for
weak force system with high resolution.

In summary, using first principles calculations, we found the electronic
structure of BLBP coupled through VdW interaction is sensitive to the
interlayer displacement. Despite of different stacking configurations, the
variation of band gap and interlayer distance $d_{z}$ follow a simple
relationship as: $E_{PeakII}\propto d_{z}^{-2}$. This simple $d_{z}^{-2}$ law
also works for other graphene-like corrugated bi-layer materials, such as
MoS2. Based on the tenability of electronic structure, we propose BLBP, as a
representative of weakly coupled 2D semiconductor family, can be used as nano-
scale displacement sensor.

This work is supported by National Key Basic Research Program (2011CB921404),
by NSFC (21421063, 11322434), by CAS (XDB01020300), and by USTCSCC, SC-CAS,
Tianjin and Shanghai Supercomputer Centers.

## References

  * (1) R. G. Knobel and A. N. Cleland, Nature 424, 291 (2003).
  * (2) C. Stampfer, A. Jungen, R. Linderman, D. Obergfell, S. Roth and C. Hierold, Nano Lett. 6, 1449 (2006).
  * (3) W. Suh, O. Solgarrd and S. H. Fan, J. Appl. Phys. 98, 033102 (2005).
  * (4) K. C. Schwab and M. L. Roukes, Phys. Today. July, 36 (2005).
  * (5) P. Hariharan, Optical Interferometry $2^{nd}$ edn (San Diego: Academic Press) (2003).
  * (6) L. K. Baxter, Capacitive Sensors: Design and Applications; IEEE Press: New York, NY, USA, (1997).
  * (7) F. Ebrahimi, Advances in Piezoelectric Transducers, Croatia, (2011).
  * (8) C. S. Smith, Phys. Rev. 94, 42 (1954).
  * (9) C. Stampfer, T. Helbling, D. Obergfell, B. Schöberle, M. K. Tripp, A. Jungen, S. Roth, V. M. Bright and C. Hierold, Nano Lett. 6, 233 (2006).
  * (10) S. N. Cha, J. E. Jang, Y. Choi, G. A. J. Amaratunga, D. CJ. Kang, D. G. Hasko, J. E. Jung and J. M. Kim, Appl. Phys. Lett. 86, 083105 (2005).
  * (11) M. Nishio, S. Sawaya, S. Akita and Y. Nakayama, Appl. Phys. Lett. 86, 133111 (2005).
  * (12) H. G. Craighead, Science 290, 1532 (2000).
  * (13) R. Lefèvre, M. F. Goffman, V. Derycke, C. Miko, L. Forró, J. P. Bourgoin and P. Hesto, Phys. Rev. Lett. 95, 185504 (2005).
  * (14) G. H. Li, A. Luican, J. M. B. Lopes dos Santos, A. H. Castro Neto, A. Reina, J. Kong and E. Y. Andrei, Nat. Phys. 6, 109 (2010).
  * (15) M. Yankowitz, J. M. Xue, D. Cormode, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, P. Jacquod and B. J. LeRoy, Nat. Phys. 8, 382 (2012).
  * (16) K. Kim, S. Coh, L. Z. Tan, W. Regan, J. M. Yuk, E. Chatterjee, M. F. Crommie, M. L. Cohen, S. G. Louie and A. Zettl, Phys. Rev. Lett. 108, 246103 (2012).
  * (17) R. W. Havener, H. L. Zhuang, L. Brown, R. G. Hennig and J. Park, Nano Lett. 12, 3162 (2012).
  * (18) G. T. de Laissardiere, D. Mayou and L. Magaud, Nano Lett. 10, 804 (2010).
  * (19) A. Luican, G. H. Li, A. Reina, J. Kong, R. R. Nair, K. S. Novoselov, A. K. Geim and E. Y. Andrei, Phys. Rev. Lett. 106, 126802 (2011).
  * (20) D. R. Hofstadter, Phys. Rev. B 14, 2239 (1976).
  * (21) R. Bistritzer and A. H. MacDonald, Proc. Natl Acad. Sci. USA 108, 12233 (2011).
  * (22) C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, T. Taniguchi, K. Watanabe, K. L. Shepard, J. Hone and P. Kim, Nature 497, 598 (2013).
  * (23) B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero and R. C. Ashoori, Science 340, 1427 (2013).
  * (24) L. A. Ponomarenko, R. V. Gorbachev, G. L. Yu, D. C. Elias, R. Jalil, A. A. Patel, A. Mishchenko, A. S. Mayorov, C. R. Woods, J. R. Wallband, M. Mucha-Kruczynski, B. A. Piot, M. Potemski, I. V. Grigorieva, K. S. Novoselov, F. Guinea, V. I. Fal ko and A. K. Geim, Nature 497, 594 (2013).
  * (25) K. H. Liu, L. M. Zhang, T. Cao, C. H. Jin, D. Qiu, Q. Zhou, A. Zettl, P. D. Yang, S. G. Louie and F. Wang, Nat. Comm. 5, 4966 (2014).
  * (26) Z. Zhu and D. Tománek, Phys. Rev. Lett. 112, 176802 (2014).
  * (27) G. Kresse and J. Furthmüller, J. Comput. Mater. Sci. 6, 15 (1996).
  * (28) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * (29) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * (30) S. Grimme, J. Comput. Chem. 27, 1787 (2006).
  * (31) S. Grimme, C. Muck-Lichtenfeld and J. Antony, J. Phys. Chem. C 111, 11199 (2007).
  * (32) J. Antony and S. Grimme, Phys. Chem. Chem. Phys. 10, 2722 (2008).
  * (33) T. Bucko, J. Hafner, S. Lebegueand and J. G. Angyan, J. Phys. Chem. A 114, 11814 (2010).
  * (34) N. Kharche and S. K. Nayak, Nano Lett. 11, 5274 (2011).
  * (35) J. Slawińska, P. Dabrowski and I. Zasada, Phys. Rev. B 83, 245429 (2011).
  * (36) R. Kagimura, M. S. C. Mazzoni and H. Chacham, Phys. Rev. B 85, 125415 (2012).
  * (37) L. Yuan, Z. Y. Li, J. L. Yang and J. G. Hou, Phys. Chem. Chem. Phys. 85, 8179 (2012).
  * (38) Y. Ma, Y. Dai, M. Guo and B. Huang, Phys. Rev. B 85, 235448 (2012).
  * (39) L. Yuan, Z. Y. Li and J. L. Yang, Phys. Chem. Chem. Phys. 15, 497 (2013).
  * (40) W. Hu, X. J. Wu, Z. Y. Li and J. L. Yang, Phys. Chem. Chem. Phys. 15, 5753 (2013).
  * (41) L. Chen, L. Wang, Z. Shuai and D. Beljonne, J. Phys. Chem. Lett. 4, 2158 (2013).
  * (42) W. Hu, X. J. Wu, Z. Y. Li and J. L. Yang, Nanoscale 5, 9026 (2013).
  * (43) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * (44) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * (45) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * (46) J. Heyd, G. E.Scuseria and M.Ernzerhof, J. Chem. Phys. 124, 219906 (2006).
  * (47) W. A. Harrison, Electronic Structure and Properties of Solids (Freeman, San Francisco, 1980).
  * (48) Q. H. Liu, L. Z. Li, Y. F. Li, Z. X. Gao, Z. F. Chen and J. Lu, J. Phys. Chem. C 116, 21556 (2012).
  * (49) J. Dai and X. C. Zeng, J. Phys. Chem. Lett. 5, 1289 (2014).