text:
\chapter{Snowmass 2013 Computing Frontier Storage and Data Management}
\label{chap:mag}

\begin{center}\begin{boldmath}

\input authorlist.tex

\end{boldmath}\end{center}

\section{Input from the Science Frontiers}
\label{sec:cpfi5-input}
\subsection{Energy Frontier}
Experiments at energy frontier hadron colliders already generate over a petabyte per
second of data at the detector device level.  Triggering and real-time event filtering
is used to reduce this by six orders of magnitude for a final rate to persistent storage
of around one gigabyte per second in the case of LHC experiments at the start of Run 2.
The cost of storing and analyzing data is a major fraction of the operational cost of LHC experiments, and thus the rate to persistent storage is determined by physics requirements and, perhaps even more, by what is affordable.

Science at energy frontier lepton colliders is unlikely to be constrained by data
management and storage issues.

The current practice of ATLAS and CMS is to treat the all the data written to persistent
storage equally through the production phases of reconstruction, re-reconstruction and
worldwide distribution of a complete set of data ready for physics analysis.  Flagging
a large fraction of the persistent data for storage on tape only, with no further
reconstruction or distribution unless a physics case arose, would cut costs or allow
the rate to persistent storage to be raised.

The vast majority of LHC data storage and data access relates to derived or simulated
data stored and accessed using ROOT\cite{Antcheva:2009zz} persistency\cite{Canal:2011zz}.  Efficient and agile data access is
already a major issue that is expected to rise in importance by the time of LHC Run 3.
Beyond HEP, efficient and agile data access will underpin future successes in data-intensive
scientific and commercial applications.  It is hard, but not impossible, to believe that HEP
will continue to be best served by a persistency solution that appears to be confined to HEP.

Distributed workflow and data management is proving to be a major success of the LHC experiments, as
illustrated in Figures \ref{fig:atlas-jobs},\ref{fig:atlas-ddm},\ref{fig:cms-aaa}, and is an area where the US brings much of the intellectual
leadership.  But this success is costly. For example the US-ATLAS M\&O-funded effort that
contributes to developing and operating the ATLAS distributed workflow and data management
system amounts to around \$4.5M per year.  This is in addition to the effort required to operate
the US Tier 1 and Tier 2 facilities.  It is difficult to imagine that the LHC-focused effort
could continue at this level for two decades, and even more difficult to imagine that other HEP
or wider science activities could take advantage of the LHC experiments current achievements
given their operational cost and complexity.
\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{ATLAS-jobs.pdf}
\caption{ATLAS jobs executing on the WLCG}
\label{fig:atlas-jobs}
\end{figure}
\begin{figure}[h!]
\centering
\includegraphics[width=0.7\textwidth]{ATLAS-DDM.pdf}
\caption{ATLAS distributed data management: Top level source-destination matrix showing success rates and data transfer rates}
\label{fig:atlas-ddm}
\end{figure}
\begin{figure}[h!]
\centering
\includegraphics[width=0.7\textwidth]{CMS-AAA.pdf}
\caption{Map showing the CMS ``Any data Anywhere Anytime'' distributed direct access system based on xrootd\cite{Bauerdick:2012st} software}
\label{fig:cms-aaa}
\end{figure}

With respect to data/physics preservation and open availability, the LHC experiments are actively
developing appropriate policies. Nascent efforts aimed at data preservation for internal
re-use and some moves towards public access have begun at the LHC.  The Tevatron experiments
have larger and more focused activities at this time; BaBar and the HERA experiments have
devloped substantial data preservation infrastructures.

\subsection{Intensity Frontier}
Lepton colliders in intensity frontier ``factory mode'' also run up against the cost of storage,
but the physics of lepton collisions is relatively clean and recording all events relevant to the
targeted physics has proved possible in the past and is a realistic expectation for the future.
The Belle II TDR estimates a data rate to persistent storage of 0.4 to 1.8 gigabytes/s, which is
comparable to LHC Run-2 rates but without the need to discard data with significant physics content.

Most of the many other intensity frontier experiments do not individually challenge storage
capabilities, but there is a recognition that data management (and workflow management) is
often inefficient and burdensome. Most experiments find it hard to escape from the comfort
and constriction of limiting all their data-intensive work to a single site -- normally Fermilab.
The statement ``all international efforts would benefit from an ATLAS-like model'' was made and
should probably be interpreted as a need for ATLAS-like data management functions at a much
lower cost and complexity than the current ATLAS system.

With respect to data preservation and open availability, the intensity frontier community
does not have a unified plan.  BaBar has implemented its Long-Term Data Analysis system in
the form of a stable, secure computing and storage cluster to support active analysis by
BaBar members for as long as possible. As a whole, the intensity frontier community recognizes
that the issue of data preservation and open availability exists and is pressing. It was widely
acknowledged that additional person-power is required to mount an extensive effort in this area.

\subsection{Cosmic Frontier}
The cosmic frontier presents several faces, each presenting its own challenges for data
management and storage: terrestrial sky survey telescopes, terrestrial radio telescopes,
HEP-detector-in-space telescopes, and large-scale simulations.

\subsubsection{Sky Surveys}
The Sloan Digital Sky Survey pioneered the use of innovative database technology to make
its data maximally useful to scientists.  This approach continues with LSST, notably the
development of a multi-petabyte scalable object catalog database that is capable of rapid
response to complex queries.  The data management needs of the sky surveys -- handling image
catalogs and object catalogs -- appear very different from those of experimental HEP, but
nevertheless, the baseline LSST object catalog employs HEP's xrootd technology in the key
role of providing a switchyard between MYSQL front ends and thousands of MYSQL backend servers.
LSST's 3.2 gigapixel camera will produce 15 terabytes per night, building up to over 100
petabytes of images and 20 petabytes of catalog database during the first ten years.

Although the basic data-access technology to make LSST science achievable has already been
demonstrated, it is certain that a vigorous LSST science community will want to attempt many
scientific studies that will be poorly served without major additional developments.
Not all LSST science will be possible using only the object catalog database. In particular,
studies such as those for dark energy effects, of particular interest to the HEP community,
are likely to require reprocessing of the LSST image data on HEP analysis facilities.
The model for funding and executing these studies is not yet clear.

The Dark Energy Survey (DES) can be considered a precursor to LSST, taking data with a
0.6 gigapixel camera for five years from 2012 culminating in a petabyte dataset.

CTA has a rather specialized real-time data challenge where some 30 gigabytes/s of data
must be gathered and processed in real time from about 100 telescopes spread over a square kilometer.

\subsubsection{Terrestrial Radio Telescopes}
Arrays of radio telescopes can present a data-volume challenge comparable with that
posed by energy frontier hadron collider experiments.  The most extreme example now being
planned is the European-led Square Kilometre Array (SKA) project that expects to complete
its Phase I system in 2020. SKA will feed petabytes/s to correlators that will synthesize
images in real time, producing a reduced persistent dataset on the scale of 300 to 1500
petabytes per year.  These volumes can only be realized if considerable evolution of
computing and storage costs happens by the time SKA data flows.  Although SKA currently
has no US involvement, it presents a concretely planned example of the technologies and
data-related challenges that will certainly be faced by US scientists involved in projects
in the same timeframe.

Today's example of the SKA concept is the Murchison Wide-Field array where a raw
15.8 gigabytes/s is processed to a produce a stored 400 MB/s.

\subsubsection{"HEP Detectors" in Space}
Examples include the Fermi Gamma Space Telescope (FGST) and the Alpha Magnetic
Spectrometer (AMS-02).  These detectors have front-end data rates far lower than
LHC experiments and would not be constrained by terrestrial storage and data analysis
capabilities.  The choke points determining their trigger rates to persistent storage
are the limited bandwidth of the downlinks that bring data back to Earth.  The
necessary conservatism applied by NASA and other space agencies to placing new
technology in space seems guaranteed to keep downlink bandwidths well below rates
that would make storage and data distribution a challenge in the future.
Nevertheless, these detectors are built and operated by large collaborations and
thus require functional distributed data management.

\subsubsection{Simulations}
Simulation provides our only way to perform ``experimental cosmology'' since only one
universe is observable.  Simulation also plays a vital role in understanding all
aspects of astrophysics, such as supernovae, for which only very limited observation
data can be collected for each occurrence. Finally, simulation is needed for the
design of observational programs and for their detailed technical elements.

Already today, post-processing of simulation data presents a major data-intensive
computing challenge, requiring data management, large-scale databases and tools for
data analytics.  Some of today's pain relates to the much more ready availability of
national resources for computation than those for data management and analysis: ``we
can easily generate many petabytes from simulations and have [almost] no place to
store them and analyze them''.

There is some expectation that compute-intensive simulation will be co-located with the
data-intensive facilities for analysis of the simulation, but powerful, easy to use
analysis tools will still need to be developed.

\subsubsection{Cosmic Frontier Data Preservation and Access}

The images, and tabular object catalogs of sky surveys and other image-based astronomy
are readily intelligible by other scientists and even the general public.

From the
experimental HEP perspective, data preservation and open availability is relatively
simple to achieve for image-based astronomy once policies have been decided.

However, the simple approach to open access will become more challenging for data-intensive image-based activities such as LSST.  Even though the data schema may be intelligible to all, analysis of the data can require massive resources, effectively limiting open access to the images or even to the multi-petabyte object catalog.  Current plans for LSST public access to
data include providing limited computing resources, with data access limited to those individuals or organizations
whose proposals for data analysis have been approved.

The raw persistent data from ``HEP-Detector-like'' devices contrasts markedly with that from the
image-based telescopes.  Like data from almost any HEP experiment, they are
intelligible only to a few experts until substantial reconstruction and analysis
has been performed.  They present the same data preservation and open availability
challenges as HEP experiments, and may be subject to the higher availability
expectations typical of image-based astronomy.

\subsection{Lattice Field Theory (LQCD)}
Like other simulation-based sciences, LQCD uses massive simulations on national
supercomputer facilities, followed by intense analysis of the resultant data.  The LQCD configuration generation step is performed on massively parallel supercomputers
because these are best adapted to the problem, whereas the subsequent, and by no means less compute-intensive, analysis step is performed on HEP-funded
throughput-optimized systems because these are the most appropriate for this step.

LQCD has significant, but not problematic data volumes that must be managed and transmitted
between the two steps, but does not face major data-related challenges.

\subsection{Perturbative QCD}
Data-related challenges are expected to remain minor in relation to those of other branches of HEP.

\subsection{Accelerator Science}
In broad summary, accelerator science is not a driver in data (or networking), but would
certainly welcome access to the easy-to-use data management and analysis tools that are
the goal of a wide range of HEP experiments.

Accelerator science has a long held dream of being able to perform predictive simulations
in close-to-real time so that feedback can be provided to physicists in the control room as
they strive to optimize accelerator performance.  A likely scenario involves running a massively
parallel simulation for a relatively short time on a remote Leadership Class Facility, followed
by the rapid transfer of tens or hundreds of terabytes of simulated data to local facilities
for rapid analysis.  This scenario is becoming achievable, but will stretch the limits of data
transmission bandwidths and of rapid data analysis.

\subsection{General Considerations for Large-Scale HEP}
Large and costly experiments or telescopes, and even some simulations, require
international collaboration.  Whatever could be done by one nation can be done
on a larger scale with more science reach as a collaborative project.  In such
international projects, the data storage and analysis can take advantage of funding
from many nations, access to large shared resources, opportunistic access to a wide
range of resources, and access to and development of distributed expertise.
The price to be paid is the complexity of a geographically distributed system.

Thus HEP, and indeed any data-intensive science operating at a comparable
international scale, needs to use a combination of wide area networks,
distributed storage, distributed computing, and the software technologies
to co-manage efficient worldwide work flow and data flow.  It seems clear
that future HEP does not need many different solutions to this challenge,
and may have much to gain by identifying commonalities with other
scientific or even commercial fields.

\section{Technology Outlook}
\label{sec:cpfi5-technology}
A simplistic prediction for the future evolution of technology is that it
would to continue to
evolve as it has in the past.  Figure \ref{fig:technology-evolution} shows a highly selective, but relevant,
view of technology evolution, charting how much you can get for \$1M in disk storage,
CPU power and long-haul network links that have been bought over three decades
for some particular experimental HEP activities.  The figure is  selective in that
it involves real purchases of hardware and services for a limited range of
experiments -- L3, BaBar, and parts of the LHC program. Some notable features of the
evolution are:
\begin{itemize}
\item
Over three decades the line ``doubling every 1.3 years'' is a good match to the average
CPU evolution and is not far from the average disk capacity evolution.
\item
The data do not exclude a marked slow down in evolution from about 2010.
\item
The network evolution shows a major discontinuity, corresponding to de-regulation
and the end of European PTT monopolies.
\item
Disk accesses per \$M -- almost totally dominated by the unchanging rotational
speed of disks -- has changed hardly at all in two decades.
\item
Comparison of technology evolution with the BaBar and ATLAS raw data rates to
persistent storage (just two examples of HEP's ``data frontier'') shows a similar
distance from the technology evolution and hence, perhaps a similar level of
technological challenge but separated by more than a decade in time.
\end{itemize}

\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{technology-evolution.pdf}
\caption{Historical trends in technology used for HEP computing.  The data reflect cpu and disk purchases for L3 at CERN, and BaBar and ATLAS at SLAC.  The network data are for the LEP3Net and USLHCNet links managed by Caltech}
\label{fig:technology-evolution}
\end{figure}

\subsection{Storage Futures}
\subsubsection{Tape}
The death of tape storage has been predicted for more than a decade, but today
its future seems to be more assured than at any time in the last decade.
Much of this change relates to the problems with disk technology described
below, but tape storage does have the intrinsic properties of negligible
power usage, a different set of failure modes and often lower failure rate
than disks, and persistently lower cost than disks.  The scientific and the
commercial world has and will continue to have, a need for archival or ``just
in case we need it'' storage with these properties.

Tape does fail, so HEP has traditionally taken advantage of its enforced distributed
approach to computing to make a copy or copies elsewhere in the world.  A cheaper
approach is possible with RAIT (redundant array of individual tapes), analogous
to RAID, provided that the data are stored away from areas prone to major natural
disasters.

Data volumes in science beyond HEP seem to be doubling every year according to NCSA
experience.  This makes sense as driven not just by the increasing data hunger of
individual sciences, but also by the increasing number of fields that are becoming
data intensive.  Hence the view from outside HEP is of a growing role for tape in an
optimized scientific data management hierarchy.

Estimating tape costs is challenging.  The smallest component is usually the tape
media, followed in ascending order by the acquisition of the drive-plus-robot system,
the maintenance of the drives and robots, and the highly skilled labor needed to
operate arcane tape-data-management systems such as HPSS.  It is reasonable to
expect that as tape consolidates its position in scientific computing, less
labor intensive approaches to its integration in the storage hierarchy will appear,
not least due to development efforts at NCSA and other major scientific computing
centers that see a clear need for this integration.

Tape technology availability has, for many decades, been limited by marketing considerations
more than by technology itself.  Tape is very far from its fundamental physical
limitations and suffers mainly from the existence of very few -- approximately
two today -- major drive manufacturers who can see no reason to compete with
themselves.  The technology could deliver twice the capacity per dollar every
18 months, but it seems reasonable to expect a doubling time of around three
years for the marketed products.

\subsubsection{Rotating Disks}
The 30-year run of exponential growth in capacity per dollar is almost certainly over.
The factors responsible are both market related and technology related.

The consumer market for rotating disks is now declining.  An increasing number of
consumer devices -- for example most HEP laptops -- are now being bought with flash
memory taking over the historical role of disk.  One consequence is that an appealing
and for-a-short-time successful HEP storage strategy -- buy the cheapest, slowest,
largest consumer-market disks and hide them behind a layer of faster enterprise disks
or solid-state-disks -- no longer works as well as it did.

In the enterprise-disk world, devices are loosely classified by their interfaces.
SATA disks with low rotational speed and maximum capacity are used for less demanding
"nearline" applications and much more expensive SAS disks are used for the more
demanding "online" applications.  The cost differences are not really related to the
SATA or SAS interfaces, but rather to properties such as rotational speed, quality of
mechanical engineering, device monitoring probes,  and automatic recovery
firmware.   The market for the most
expensive disk capacity, currently provided by 2.5 inch 15k rpm drives used
primarily for database applications, is already threatened by solid-state
storage with much lower latency.

The technological problems faced by magnetic disk are easy to understand. The area
of the bits written to current 4 terabyte drives is close to the limit of magnetic
stability.  The density could be raised if more atoms could be involved by writing
bits that used more of the magnetic material below the surface of the disk platter,
and/or by using a higher coercivity material that would allow smaller bit sizes.
Technological developments along these directions have been in the works for some
years, but according to industry experts, none will be advanced enough to bring the
next jump in disk capacity to market much before the end of the decade.  The prime
``use-a-larger-volume-of-magnetic material'' approach is shingled recording,
conceptually illustrated in Figure \ref{fig:shingled-recording}. In this approach a specially shaped write
head lays down tracks with a triangular cross section, which become a set of
``shingles'' as successive tracks are written.  Each shingle goes deep into the
material, so can be stable even if it is much narrower than current tracks.
It is impossible to re-write sectors or single tracks with this approach.
To change any information on the disk, a large block of tracks must be erased
and re-written.  HEP would have no problem with such a device, but its success
in the wider marketplace could be problematic.  The ``high-coercivity'' approaches
include Heat-Assisted Magnetic Recording (HAMR) using a uniform surface with higher
coercivity and ``bit patterned recording'' using lithography to lay down a circular
pattern of bits that are close together but physically isolated from each other.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{shingled-recording-crop.pdf}
\caption{Disk track layout produced by shingled recording.  Tracks are written from left to right in the diagram by a head that produces a magnetized volume with a triangular cross section. Individual tracks cannot be modified or erased}
\label{fig:shingled-recording}
\end{figure}

The next ten years of disk evolution is thus highly uncertain.  It would be
prudent and probably correct to assume that the doubling time is now around four
years.  The impact of this will be considerable.  The logic that made it
cost-effective to replace storage after as few as three years (approximately
twice the doubling time) might now argue for keeping disks eight years or more,
if the disks were sufficiently reliable.  Disks may also last longer if most of
the access traffic can be fielded by a substantial solid-state storage cache as
part of the overall storage hierarchy. There will be strong market pressures
favoring, long lasting, low power, physically dense, and therefore very heavy
disk systems.

\subsubsection{Solid-State Storage}
The certainty of an increasing role for solid-state storage is suggested by
the twenty-year failure of rotating disk to provide more accesses per second
per dollar.  In terms of this metric, today's solid-state storage provides
vastly better value than rotating disk.  The solid state storage being used
in (exploratory) production at HEP computing centers costs about ten times as
much per unit capacity as rotating disk.  Reducing a rotating disk purchase
by 10\% and spending the money on solid-state cache is cost neutral and can
improve physics analysis throughput.

In this caching mode, solid-state storage is certain to play an increasing but not
dominant role in HEP, and perhaps an even larger role in the commercial world.
The woes of rotating disk described above are likely to allow solid-state --
specifically flash-memory-based -- storage to lower its cost relative to disk
and gain market share.  It seems plausible that the cost differential would
go down from around a factor ten now to around a factor 3 within ten years.
Industry experts do not expect current solid-state technologies to "kill" disk
in the foreseeable future, not least because it is almost inconceivable to
create the large number of multi-billion dollar ``fabs'' -- chip fabrication
facilities -- that would be needed to displace the world's many exabytes of disk.

\subsection{Storage Middleware and Data Management}
With the notable exception of the adoption by several HEP sites of the
semi-commercial HPSS (High Performance Storage System)\cite{HPSS:url} to manage tape data,
single-site storage middleware in HEP has been almost entirely home made.  In the
late 1990s, it seemed that it would be an achievement to limit HEP's invention
of tape data management systems to only one per lab.

HEP's middleware stacks have proved quite successful, for example the seamless
integration for BaBar at SLAC of petabyte (HPSS) tape storage and tens of terabytes
of disks distributed over 50 servers, but these successes have not translated into
wide HEP adoption, let alone interest by other sciences or industry.

In the wider scientific world, middleware such as Globus Online (GO)\cite{GlobusOnline:url}, which is
a data transfer mechanism that has nice retry and a graphical display, is
currently the norm.  These tools continue to grow and change, but are still
based on the gridftp protocol.  There are other tools for data transfer, but
it still has to be managed by the application.  Tools are also emerging that
that transfer data at the file system level.

Scalable file systems such as Lustre\cite{Lustre:url} and GPFS\cite{GPFS:url} connect
the file system with the archive or nearline environment so there is seamless
access, from the users prospective, to all files.  The data-managed file
system means that the inode remains in the file system with a stub of the first
blocks of the data, or maybe nothing of the file at all, and the I/O is captured
until the data returns to the local disk of the file system server.
The data can be recalled from any ftp archive as long as permissions and
security have been prearranged.   This builds basically an endless file system.
Connecting these large endless file systems is then next on the horizon.

In the next few years some test sites will be using GO-Storage
to connect these large managed file systems together. The meta-data will be
separate, in a replicated environment, and the data could be anywhere in the
world.  As data is retrieved from locations, it moves closer and locally
depending on the application requirements.   These methods all have small
latencies added and depending on the actual location of the data, the latency could
be great especially if it's on tape somewhere.   If these efforts yield systems and services
that are widely adopted by science they must become serious candidates for
meeting HEP needs in the future.

Even so, at the current time there are no commercial or widely used open-source
offerings meeting the HEP needs for worldwide data and workflow management.
The absence of commercial offerings is probably largely due to HEP's need to
integrate tens to hundreds of autonomous computer centers, having many different
funding sources, many different technologies, and a varying level of affiliation,
from loose to none at all, with any particular HEP project.  Multinational commercial
entities are generally orders of magnitude more coherent in their technologies and
management.  The wider world of collaborative science does face many of the HEP
challenges, but HEP is still at the bleeding/leading edge in terms of overall
complexity.

The commercial world of distributed data-intensive applications should not be
ignored -- individual technologies and {\it de-facto} standards relevant to HEP are
almost certain to appear.  The wider scientific world must also not be ignored.
HEP will be able to bring benefit, reap credit and attract non-HEP funding by
helping adapt some of its most successful distributed computing technologies
for other sciences and smaller HEP activities. Beyond this, HEP can benefit
from participating in the development of widely useful distributed computing
tools, and should studiously avoid the arrogance that could postpone the adoption
of tools developed outside HEP to serve the exponential growth of data-intensive science.

A major stimulus to obey the exhortations above must be the labor intensive aspects,
in both ongoing development and operations, of the distributed computing systems
used in HEP today.

\section{Data, Software and Physics Preservation}
It is widely recognized that "knowledge preservation" in all frontiers is becoming
increasingly important, or even vital, as the size of data samples and the time
scale of experiments increase.  Here, ``knowledge preservation'' includes the preservation
of data, the accompanying software, the capability to apply the software to the data
in a reproducible way, and a record of how the software was applied to obtain a given
result.  There are several motivations for effort in this area:
\begin{itemize}
\item
Enlightened self-interest: foresight and infrastructure is required to be able to
re-use data and reproduce physics results several years later, even within a single
experiment.  The LHC experiments are already seeing limitations in this area, even
after only three years of data-taking.  Solving this problem in a generic way will
enable most aspects of the issues behind knowledge preservation to be addressed in
a manner that will allow many different experiments across the different frontiers
to adopt a common solution.
\item
Outreach: limited access to some data and the accompanying documentation provides one
means of preserving the ensemble of data and software.  Resources such as Rivet\cite{Rivet:url},
RECAST\cite{RECAST}, and the extensions to HepDB provide a means of encapsulating the structure
of a given analysis in a manner that allows easy communication between experts and
relatively simple re-use.  More broad outreach efforts, such as those supported by
nearly all Energy Frontier experiments, have also developed simple yet powerful
``public-friendly'' interfaces to simplified experimental data and analysis tools.
\item
Mandate: stewardship of the public investment in fundamental research is very important.
(Have to write this carefully) While this is recognized by many of the experiments,
HEP as a field has been reluctant to make large datasets public due their inherent
complexity and the perceived difficulty in their interpretation by the non-expert.
However, there is increasing pressure to make some the data available to the public
in some form.  Doing this in a manner that is sustainable and does not require immense
resources poses many difficult problems which would have to be solved if this pressure
turns into a mandate from the funding agencies.
\end{itemize}

\subsection{Current Efforts in Knowledge Preservation}
Worldwide, and across many disciplines, the past couple of years have seen widespread
recognition that Data/Knowledge preservation is critical to the long-term exploitation
of any science that qualifies as ``Big Data''.  This broad interest has spawned a
dizzying array of national and international disciplinary and cross-disciplinary
efforts to attempt to address the related issues.  Up to now, many of these have
focused solely on providing a means of storage and indexing of datasets.  This is
insufficient for many of the particle and astrophysics experiments, who require a
means of recording processing workflows as well.  Indeed, particle experiments
probably have the most stringent and complex requirements for knowledge
preservation, driven primarily by the complexities of the datasets and the amount
of processing required to achieve a physics result.

The DPHEP effort\cite{DPHEP}, a Study Group for Data Preservation in High Energy Physics under the auspices of the International Committee on Future Accelerator (ICFA), has produced a study outlining the current state of data preservation within HEP, including an extensive overview of other disciplines which will not be reproduced here. They suggest a series of guidelines for HEP data preservation efforts, as well as a framework for global coordination.  Their conclusions include a recognition of the scientific potential for data re-use, especially the desirability to preserve full analysis capability.  They also emphasize the urgency required to begin and sustain global, coordinated data preservation efforts.

Among colliders whose data-taking period has ended, two currently have advanced efforts in data preservation and access.  The BaBar Long Term Data Access archival system\cite{BaBar} is comprised of several racks of compute nodes and data and database servers, all completely walled-off from the host laboratory.  The deployment of virtual machines ensures a stable, non-evolving software environment.  This system is designed to provide active analysis use through 2018, at which point only the legacy data would be preserved.  At DESY, a combined effort is underway to maintain sustained access to the HERA datasets of ZEUS, H1, and HERMES.  A cross-HERA virtual environment based on current grid infrastructure that will allow new production of simulated events already exists; each of the experiments has slightly different plans for the long-term archiving and access to their data and the corresponding documentation.  A general software/archival validation tool is also under development at DESY\cite{DESY}; this could become a general component in any long-term software and data preservation effort.
Planning at the Tevatron is underway\cite{Tevatron}. The Tevatron experiments finished data collection at the end of September, 2011. Unlike at the LHC, very little thought was given to data preservation and re-use during Tevatron running, which will make the preservation of the data especially challenging.

At the LHC, K. Cranmer's RECAST project is notable for its advanced state of deployment and its differences from traditional archival methods.  RECAST is a computational framework designed to encapsulate the "institutional knowledge" about a particular analysis in a way that allows outsiders to query the dataset for a certain physics result.  The analysis methods, the simulated events, and the results from data are stored in a manner that can be used to probe the presence or absence of new physical processes that were not considered at the time of the analysis.  So, the analysis algorithms and the dataset are fixed, but the results and constraints on new physics models can be explored.

The LHC experiments are currently considering data preservation and public access policies.  CMS and LHCb have approved policies on public access to data and on the prospects of releasing processed data to the public.

The Data and Software Preservation for Open Science (DASPOS) project, funded by the National Science Foundation, has established a collaboration tying together US physicists from the CMS and ATLAS experiments at the LHC and the Tevatron experiments, scientists from other disciplines, and experts in digital preservation, heterogeneous high-throughput storage systems, large-scale computing systems, and grid access and infrastructure. The DASPOS activities are connected with the DPHEP coordination effort, the experimental collaborations, and other related multi-disciplinary projects in Europe, Asia, and the US. Together, this group represents the US in international efforts, and acts as a coordinating point of contact, a partner in dialogue, and a technological consort.  The DASPOS project represents an initial exploration of the key technical problems that must be solved to provide appropriate data, software, and algorithmic preservation for HEP and will address the importance of preserving the contexts necessary to understand, trust, and re-use data. It will attempt to create a template architecture for
knowledge preservation in data-based science, focused on HEP, but representing broader consultation and applicability.  The
activities of the project comprise a series of workshops aimed at fact-finding and cross-disciplinary exploration coupled with a technical demonstration of a knowledge preservation infrastructure.

While DPHEP serves as a central forum for discussion of Data/Knowledge preservation in high energy physics, there are many other efforts in other scientific fields that are attempting to grapple with similar problems.  This is one area where central guidance and coordination of expertise across funding agencies and scientific entities would be very welcome.

\section{HEP Outlook and Recommendations}
\label{sec:cpfi5-hep-outlook}
\subsection{Impact of the Technology Outlook on HEP}
The coming decade and beyond will see experiments and observations at several
frontiers straining against the limits of data management and storage technology.
Storage technology is likely to evolve in capacity/cost more slowly than in the
last decades, making it ever more important that the role of storage is carefully
optimized with respect to other costly elements of the scientific programs.
Figure \ref{fig:storage-evolution} attempts to illustrate this challenge by comparing the data rates
expected from the most challenging activities with the disk technology evolution
predicted earlier in this report.

\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{storage-evolution-v2.pdf}
\caption{Raw data rates of HEP experimental/observational programs in comparison with likely capacity/unit price evolution}
\label{fig:storage-evolution}
\end{figure}

As outlined earlier, the rate of writing raw events to persistent storage for the
LHC experiments is substantially influenced by the cost of storage and can be
expected to evolve at a similar rate.
It is likely that other activities with petabyte/s rates for raw data from
front-end devices will be compelled to employ real-time data processing and
selective rejection in order to reach persistent rates that are affordable
in the context of each activity.

At a more detailed level, the roles of tape, rotating disk and solid-state
storage will evolve:
\begin{itemize}
\item
Tape will fall slightly in relative cost should and play a more important role
than, for example, in the LHC experiments today where it is arguably underused.
\item
Rotating Disk, will not improve in capacity/cost as rapidly as tape or solid-state
storage and for the next decade may take around four years to double capacity/cost.
As a result, it could be cost-effective to wait about eight years before replacing
disks that are not giving trouble, provided the space, power and cooling are available to make this possible.
\item
Solid-state storage will not kill rotating disk in the enterprise and data-intensive
science markets, but it will become relatively cheaper, perhaps by a factor of three,
and will merit full consideration in the overall optimization of HEP computing
environments, where it offers the promise of efficient sparse or random access
to HEP data.  To make effective use of the still small affordable quantities of
solid state storage, application-aware approaches to caching data on such storage
appear essential.
\end{itemize}

The relative rates of cost evolution of tape, rotating disk, solid-state storage,
networks and CPU are clearly uncertain.  However, it is certain that HEP computing
models will need to adapt to make the most efficient use of all these elements.
The best way to adapt would be to make the implementation of the models
intrinsically adjustable to serve a wide range of situations. The now old,
but never really implemented, concept of ``virtual data'' would go a long way to
making HEP computing dynamically adaptable to wide range of relative storage and
CPU costs.  In a virtual data system, all derived (or simulated) data products
begin existence solely as the rigorously complete recipes for creating them.
These virtual data products are then instantiated or replicated based on real
or algorithmically anticipated demand.  The physical instances are retained
based on anticipated future use and the relative cost of storage versus re-creation.
HEP has often appeared to be close to developing the rigorous provenance databases
needed to implement virtual data and a full virtual data implementation may be in
reach of the LHC experiments for Run 3.

\section{Findings and Recommendations}
\label{sec:cpfi5-findings}
{\bf CpF15 Finding 1a:} The largest HEP experiments have developed, and are
improving functional distributed data and workflow management systems meeting
their needs. These systems are expensive to develop and operate and are thus
rarely appropriate for smaller experiments.

{\bf CpFI5 Finding 1b:} HEP currently benefits from, but can also be
constrained by, the highly successful ROOT features supporting reading
and writing of persistent data. No other major scientific field uses
ROOT or appears interested in it. Major developments in persistency
technology will be required to take advantage of storage hardware on the
timescale of LHC Run 3.

{\bf CpFI5 Recommendation 1:} HEP should maintain and promote a vision of the future
in which fully functional and low-operational-cost distributed computing and
persistency management is supported by software that is widely used in
data-intensive science.  To this end, developments in industry and the wider
science community should be monitored actively, HEP should work with the wider
science and computer science community to export and adapt HEP technologies and
vice-versa. In distributed computing, HEP should organize itself to significantly
reduce the number of diverse approaches and provide the benefits of ideas and
software developed in the largest experiments to other activities where they are needed.

{\bf CpFI5 Finding 2a:} Rotating disk storage will suffer a marked slowdown in the
evolution of capacity/cost.  This may be the largest perturbation of HEP computing
models that must attempt to optimize the roles of tape, rotating disk, solid-state
storage, networking and CPU.

{\bf CpFI5:} Finding 2b: Many of the components required to support virtual data
already exist in the data and workflow management software of the largest experiments.
The rigorous provenance recording required to support the virtual data concept would
also benefit data preservation.

{\bf CpFI5 Recommendation 2:} Computing model implementations should be flexible
enough to adapt to a wide range of relative costs of the key elements of HEP
computing.  In preparing for Run 3, the LHC program should seriously consider
virtual data as a way to accommodate scenarios where storage for derived and
simulated data becomes relatively very costly.

{\bf CpFI5 Recommendation 3:} All experiments across all frontiers need
infrastructure that will allow scientists to store, catalogue, access, and
reprocess datasets years after the original physics results are produced.
The inherent similarity of the requirements across experiments and disciplines
call for a coordinated investment in common infrastructure to enable easy
access and adoption of best practices in knowledge preservation.  Solutions
should be developed that meet the needs of the particle and astrophysics
communities before widespread release of data to the public can be expected
or mandated.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
Quantum optomechanics has become an active area of research in which the quantized center-of-mass motion of a tiny mirror plays the central role.
One aim is to probe the existence of coherent superpositions of macroscopically different quantum states and study their decay mechanisms \cite{marshall2003towards,kleckner2008creating}. A quantum mirror may also be used for quantum information processing purposes, e.g., to
store optical information in the mechanical motion
of the mirror. By mapping this information back into a different optical mode, frequency conversion
of photonic quantum states may be achieved \cite{tian2010optical}.

So far, the quantized motion has been studied experimentally by means of laser light interacting with the moving mirror.
Here we study the possibilities that arise from having a non-negligible interaction between a single photon and a mirror (a mildly futuristic possibility, but one that starts being taken very seriously \cite{nunnenkamp2011single,pepper2012optomechanical,hammerer2012nonclassical,akram2013,sekatski2014macroscopic}).

In particular, we will describe theoretically how the {\em time-dependent  spectrum} of a single photon is modified by its interaction with a moving mirror. The time-dependent spectrum of a single photon is defined as follows:
Imagine that a single-photon wavepacket enters a (Lorentzian) filter cavity, described by a resonance frequency $\omega$ and a filter bandwidth $\Gamma$. We can record as a function of time when the photon exits the filter. If we represent the detuning between filter resonance and optomechanical cavity frequency $\omega_{c}$ by $\Delta=\omega-\omega_{c}$, then we can express the expected counting rate at time $t$ in terms of the continuous field annihilation and creation operators $\hat{a}(t)$ and  $\hat{a}^\dagger(t)$ of the photon wavepacket as
\begin{equation}\label{TDPS}
\begin{split}
&N(t;\Delta,\Gamma)=\Gamma^{2}\int_{0} ^t\int_{0} ^t e^{-(\Gamma-i\Delta)(t-t^{'})}e^{-(\Gamma+i\Delta)(t-t^{''})}\times\\
&\hspace{20mm}\langle\hat{a}^{\dagger}(t^{'})\hat{a}(t^{''})\rangle dt^{'}dt^{''}.
\end{split}
\end{equation}
This is the time-dependent spectrum (it depends both on $\omega$ and $t$). It has the same form as that introduced by Eberly and Wodkiewicz \cite{eberly1977time}. The main difference is that the classical field amplitudes in the Eberly and Wodkiewicz spectrum are here replaced by quantum annihilation and creation operators.

We can also integrate the time-dependent spectrum over time, and thus define
 \begin{equation}\label{SPS}
 N_S(t;\Delta,\Gamma)=
 \int_0^t N(t';\Delta,\Gamma)dt'.
 \end{equation}
 In the limit of $t\rightarrow\infty$ this quantity would equal the spectrum for a stationary process as obtained from the Wiener-Khinchine theorem \cite{mandel1995optical}, which is perhaps the more familiar quantity.

We are going to use two theoretical methods in the following. First, the whole process of detecting
a single photon emanating from a cavity is very well described by the quantum trajectory method \cite{molmer1993monte, dum1992monte, wiseman1996quantum, plenio1998quantum,carmichael2008statistical}, especially when combined with input-output theory \cite{gardiner2004quantum, ciuti2006input}.
As we have shown before \cite{mirza2013single}, the time-dependent single photon spectrum (as well as its infinite-time limit) can be straightforwardly calculated using the method developed in \cite{carmichael2008statistical}. Second, we find that a simple dressed-state picture
suffices to understand the locations and heights of the resonances that become visible in these spectra.

\section{Theoretical description}
\begin{figure*}[t]
\hspace{-20mm}\includegraphics[width=4.4in]{Fig1.pdf}
\captionsetup{
  format=plain,
  margin=1em,
  justification=raggedright,
  singlelinecheck=false
}
  \caption{An optomechanical cavity (OMC): A Fabry-Per\'ot cavity with a fixed mirror on the left and a movable mirror on the right. Initially the OMC is assumed to contain a single photon and zero phonons. The right mirror (modeled as a harmonic oscillator) is coupled to a finite temperature heat bath with mechanical decay rate $\gamma_{M}$, while the temperature of the bath is given in terms of the average thermal phonon number $\overline{M}$. The decay rate of the optical cavity field amplitude (due to leakage of photons through the left mirror) is denoted by $\kappa$. The output field is detected and its time-dependent spectrum is used to probe the motion of the movable mirror.}\label{Fig1}
\end{figure*}
As shown in FIG.~1 we consider an optomechanical cavity (OMC) system \cite{aspelmeyer2013cavity, kippenberg2008cavity, park2009resolved}. The OMC is a Fabry-P\'erot cavity with a perfectly reflecting movable right mirror, which is modeled as a quantum harmonic oscillator with frequency of oscillation $\omega_{M}$. The annihilation of quanta of mechanical vibrations (phonons) is described by the operator $\hat{b}$. Initially the OMC is assumed to have a single photon inside, and we assume  that only a single resonant optical mode is relevant. The frequency of that mode is $\omega_{c}$,  and annihilation of photons in that mode is described by the operator $\hat{a}$. We study the situation where the right mirror is so thin that even a single photon can affect its mechanical motion through radiation pressure \cite{jackson1998classical, girvin2009trend}. The rate describing the coupling strength between optical and mechanical degrees of freedom is denoted by $g_{M}$.

We shall start by writing down the standard Hamiltonian \cite{aspelmeyer2013cavity, milburn2011introduction, meystre2012short} of such an  OMC while leaving out the zero-point energies of the oscillators, so that
\begin{equation}\label{Hsys}
\hat{H}_{sys.} = \hbar\omega_{c}\hat{a}^{\dagger}\hat{a}+\hbar\omega_{M}\hat{b}^{\dagger}\hat{b}-\hbar g_{M}\hat{a}^{\dagger}\hat{a}(\hat{b}+\hat{b}^{\dagger}).
\end{equation}
The  left mirror is taken to be partially transmitting and hence  the photon wavepacket can leak out through that left mirror. We assume this happens at a rate $\kappa$ and that the photon wave packet enters an optical fiber (which is assumed to have a continuum of modes). The fiber is introduced just for the purpose of ensuring unidirectional propagation of the emitted photon wavepacket towards a detector. The transmission, reflection and emission spectra of a single photon wavepacket from either a one-sided or a two-sided Fabry-P\'erot cavity (with both mirrors fixed) have all been very well studied \cite{ujihara2010output, saleh1991fundamentals}, but here we address the  situation where one of the mirrors of the OMC is capable of harmonic motion.

Due to the presence of the continuum of modes in the fiber, our problem is essentially an open quantum system problem. There are several different approaches which can be applied in order to calculate the spectrum emitted by such a system \cite{carmichael1993open, breuer2002theory, weiss2012quantum}. Here we will make use of the Quantum Jump/Trajectory (QJT) approach \cite{molmer1993monte, dum1992monte, wiseman1996quantum, plenio1998quantum,carmichael2008statistical} combined with input-output theory \cite{gardiner2004quantum, ciuti2006input}.

There is one output detector in our system and we introduce the corresponding output operator $\hat{J}_{out}$ (which is basically the continuous-mode annihilation operator for the output field). This output operator  is related to the input operator $\hat{J}_{in}$  through the standard input-output relationship \cite{gardiner2004quantum} as:
\begin{equation}\label{JO}
\hat{J}_{{\rm out}}(t)=\hat{J}_{{\rm in}}(t)+\sqrt{\kappa}\hat{a}(t)
\end{equation}
where we have neglected the trivial fiber time delays between the OMC and the detector (following the standard cascaded quantum jump approach \cite{gardiner1993driving}). We also can often disregard the input operator in the above equation as it is not going to contribute to normally ordered observables (all our observables of interest are of that form). This is because if we denote by $\ket{\Psi}$ the initial state of the global system (cavity and fiber) we have $\hat{J}_{in}\ket{\Psi} = 0$ as we assume there is no (fiber) photon present to serve as input for the OMC.
\section{Single-phonon mechanical oscillations}
In this section, we consider a rather simple situation where we assume that initially there is no phonon present and the photon inside the OMC will generate at most a single phonon.
\subsection{Absence of mechanical losses and zero temperature}
We first consider the simplest case of no mechanical losses and zero temperature for the mechanical heat bath. That is,   we assume a zero phonon leakage rate ($\gamma_{M} = 0$), and we assume the average thermal phonon number $\overline{M}$ to be zero. Note that the first condition can always be physically realized for short times $t$ such that
\[
\frac{1}{\kappa}<t<\frac{1}{\gamma_{M}}.
\]
In the next subsection we will discuss the effects of mechanical losses and nonzero temperature on the behavior of the spectrum.

According to the QJT approach we have the following picture. In any given small time interval we have one of two
situations:

(I) Occurrence of a quantum jump: whenever our detector records a click, a quantum jump takes place, and we apply the output (annihilation) operator $\hat{J}_{out}$ (which in this context is also termed the ``jump operator'')
to the state of the system.

(II) No quantum jump takes place: while no detector clicks, the system evolves according to the following non-Unitary Schr\"odinger equation:
\begin{equation}\label{NUSE}
i\hbar\frac{d}{dt}\ket{\tilde{\psi}(t)}=\hat{H}_{NH}\ket{\tilde{\psi}(t)},
\end{equation}
where the ``Hamiltonian'' appearing in this equation is a non-Hermitian operator, which is the sum of two parts. The first part is Hermitian and is given by the system's Hamiltonian, the second part is anti-Hermitian and  is constructed from the jump operator. In total we have
\begin{equation}\label{NHH}
\begin{split}
&\hat{H}_{NH}=\hat{H}_{sys}-i\hbar\hat{J}^{\dagger}_{out}\hat{J}_{out}\\
&=\hbar\omega_{c}\hat{a}^{\dagger}\hat{a}+\hbar\omega_{M}\hat{b}^{\dagger}\hat{b}-\hbar g_{M}\hat{a}^{\dagger}\hat{a}(\hat{b}+\hat{b}^{\dagger})-i\hbar\frac{\kappa}{2}\hat{a}^{\dagger}\hat{a}.
\end{split}
\end{equation}
Restricting the phonon number to be at most 1 (one) we can also define a mechanical annihilation operator as $\hat{b}=\ket{0}_{b}\bra{1}$, in terms of the states $\ket{0}_{b}$ and $\ket{1}_{b}$, the zero- and one-phonon number states, respectively. This restriction also implies that $\hat{b}^{\dagger}\ket{m}=0$, $\forall$ $m\geq 1$. Notice that with this new notation the Hermitian nature of $\hat{H}_{sys}$ will not be disturbed.

The unnormalized state $\ket{\tilde{\psi}(t)}$ appearing in Eq. [\ref{NUSE}] is called the ``No-jump state'' in QJT and it is a superposition of all the different possibilities of finding the excitation in the system before it is being lost by the system and registered by the detector. For our setup it can be written as:
\begin{equation}\label{NJS}
\ket{\tilde{\psi}(t)} = c_{1}(t)\ket{10}+c_{2}(t)\ket{11},
\end{equation}
where we use the following  notational convention: the first place in the ket gives the number of photons in the OMC and the second place the number of phonons in the mechanical oscillator. The probability amplitudes appearing in the No-jump state can easily be worked out by using Eqs.[\ref{NJS}] and [\ref{NHH}] in Eq.~[\ref{NUSE}]. In Laplace space these amplitudes then turn out to be
\begin{subequations}
\begin{eqnarray}
C_{1}(s) = \Bigg[\frac{s+\frac{\kappa}{2}+i\omega_{M}}{(s+\frac{\kappa}{2})(s+\frac{\kappa}{2}+i\omega_{M})+g^{2}_{M}}\Bigg],\\
\label{C1s}
C_{2}(s) = \Bigg[\frac{ig_{M}}{(s+\frac{\kappa}{2})(s+\frac{\kappa}{2}+i\omega_{M})+g^{2}_{M}}\Bigg],\label{C2s}
\end{eqnarray}
\end{subequations}
where $C_{i}(s)$ is the Laplace transform of $c_{i}(t)$ with $i = 1, 2$.

For the spectrum calculations in QJT,  we first have to calculate the time-dependent spectrum and then by taking the $t\rightarrow\infty$ limit we can obtain the infinitely long time spectrum. We use the equations given in the Introduction, with the generic operator $\hat{a}$ replaced by $\hat{J}_{out}$. The latter is determined by Eq.~[\ref{JO}], and
expectation values involving $\hat{a}$ and $\hat{a}^\dagger$ are determined by the coefficients $C_i(s)$.
Our infinitely long time spectrum $P_{D}({\Delta,\infty})$, can be evaluated analytically, and the result is
\begin{equation}\label{ILTS}
\begin{aligned}
&\hspace{-5mm}P(\Delta,\infty)=\kappa\Gamma\Bigg(\Bigg|C_{1}(s=-i\Delta)\Bigg|^{2}+\Bigg|C_{2}(s=-i\Delta-i\omega_{M})\Bigg|^{2} \Bigg)\\
&\hspace{8mm}=\kappa\Gamma\Bigg[\Bigg|\frac{i(\omega_{M}-\Delta)+\kappa/2}{\lbrace i(\omega_{M}-\Delta)+\kappa/2 \rbrace \lbrace \kappa/2-i\Delta \rbrace+g^{2}_{M}}\Bigg|^{2}\\
&\hspace{10mm}+\Bigg|\frac{ig_{M}}{\lbrace \kappa/2-i(\omega_{M}+\Delta)\rbrace \lbrace \kappa/2-i\Delta \rbrace+g^{2}_{M}}\Bigg|^{2}\Bigg].
\end{aligned}
\end{equation}
This is the main result for this subsection. We plot this spectrum  in various different regimes:\\
\begin{figure*}
\hspace{7mm}\vspace{-5mm}\text{{\bf (a)}}\hspace{80mm}\vspace{-2mm}\text{{\bf (b)}}
\begin{center}
\begin{tabular}{cccc}
\subfloat{\includegraphics[width=3in,height=6in]{Fig2a.pdf}}&
\hspace{10mm}\subfloat{\includegraphics[width=3in,height=6in]{Fig2b.pdf}}\\
\end{tabular}
\captionsetup{
  format=plain,
  margin=1em,
  justification=raggedright,
  singlelinecheck=false
}
\caption{Single photon spectrum in the limit $\rightarrow\infty$ as emitted by the OMC, in the case that there is at most one phonon present.
The common parameters used in both parts (a) and (b) are the frequency filter parameters $\Delta/\omega_{M} = 0.5$ and $\Gamma/\omega_{M}=0.1$.
{\bf (a)} Comparison of the good cavity ($\kappa<\omega_{M}$) and bad cavity ($\kappa > \omega_{M}$) regimes. From top to bottom curves one sees that side bands in the spectrum become better and better resolved. In all curves $g_{M}/\omega_{M} = 1.25$.
 {\bf (b)} Varying the values of $g_{M}$ as compared to $\kappa$ (in units of $\omega_{M}$) so that a comparison between the strong and weak coupling regimes can be made. For all curves we have $\kappa/\omega_{M} = 0.25$.
In both parts (a) and (b) of the figure  the separation between the two middle peaks is determined by $g_{M}$ while the two side peaks are $\omega_{M}$ farther away from the middle peaks.}\label{Fig2}
\end{center}
\end{figure*}
{\bf Good and Bad cavity limits:} First we plot this spectrum in the good ($\kappa < \omega_{M}$) and bad ($\kappa > \omega_{M}$) cavity limits \cite{marquardt2007quantum, marquardt2009optomechanics}. FIG.~[2(a)] shows our results. For the parameters chosen (especially the value of $g_{M}/\omega_{M}$ ) we note that even in the bad cavity limit there are some resonant structures visible, even though all resonances are overlapping. Gradually going to the good cavity limit and finally approaching $\kappa = 0.10 \omega_{M}$ we see that all resonances are now separated and can easily be distinguished from one another.\\
{\bf Strong and weak coupling regimes:} Next we vary $\kappa$ and $g_{M}$ in units of $\omega_{M}$. In FIG.~[2(b)] we have plotted our results in both weak $(\kappa>g_{M})$ and strong $(\kappa<g_{M})$ coupling regimes \cite{groblacher2009observation, weiss2013strong}.
We note that starting from the weak coupling regime (top red curve in the figure) we have just one major peak centred at the optical resonance frequency and a very tiny shoulder (side-band) on the left side of the peak.  This left side shoulder (the ``red side band'' \cite{chan2011laser}) indicates a process in which the photon first produces a phonon and then leaves the cavity with a frequency smaller than $\omega_c$.
 For the red curve, the red side band is very small because the coupling between the single photon and the single phonon is kept small compared to the photon escape rate $\kappa$. In the next three curves (blue, pink and black) we enter into the strong coupling regime. With this change we start to observe all side bands clearly. Specifically, we now see the appearance of blue side bands as well (on the positive side of the $\Delta$ axis). These new side bands refer to the processes in which the photon first produces a phonon and then takes energy from away from that phonon and leaves the cavity with energy greater than $\hbar\omega_c$.

From this figure we can conclude that a fully resolved OMC spectrum (the black curves in FIG.~2) can be observed if we work in the strong coupling regime within the good cavity limit.
\subsubsection{Locations of the resonances}
The locations of the resonances can be worked out by setting the real part of the poles in the spectrum of Eq.~[\ref{ILTS}]
equal to zero. There are two terms in the spectrum and each term gives us two resonances (poles).  We thus find four resonances determined by
\begin{equation}
\Delta=\pm\frac{\omega_{M}}{2}\pm\sqrt{\frac{\omega^{2}_{M}+\kappa^{2}}{4}+g^{2}_{M}}.
\end{equation}
The locations of the resonances can also be worked out by performing a dressed-states analysis \cite{barnett2002methods} of this problem.

\begin{figure*}
\hspace*{-3cm}%
\centerline{

  \resizebox{15cm}{!}{
    \begin{tikzpicture}
     \draw[thick](-5,0) -- (-3,0) node at (-5,0)[left] {{$0$}}node at (-2,0)[left] {{$\ket{0,0}$}};
     \draw[thick](-5,2) -- (-3,2) node at (-5,2)[left] {{$\hbar \omega_{c}$}}node at (-2,2)[left] {{$\ket{1,0}$}};
     \draw[thick](-5,2.5) -- (-3,2.5) node at (-5,2.5)[left] {{$\hbar (\omega_{c}+\omega_{M})$}}node at (-2,2.5)[left] {{$\ket{1,1}$}};
     \draw[thick](-5,0.5) -- (-3,0.5)node at (-5,0.5)[left] {{$\hbar\omega_{M}$}}node at (-2,0.5)[left] {{$\ket{0,1}$}}node at (-4,3.75){{{$ \small \text {Bare/Uncoupled states}$}}}node at (5.5,4){{{$ \small \text{Dressed  states}$}}};
    \hspace{20mm}\draw[thick](1.8,0) -- (5,0) node at (6,0)[left]{$\small{\ket{0,0}}$}node at (1.4,0)[right]{$0$};
                 \draw[thick](1.8,0.5) -- (5,0.5)node at (6,0.5)[left]{$\small {\ket{0,1}}$}node at (1,0.5)[right]{${\rm \hbar\omega_{M}}$};
                 \draw[thick](1.8,1.5) -- (5,1.5)node at (8,1.5)[left]{$\small {C_{2}\ket{1,0}+C_{1}\ket{1,1}}$}node at (-2.5,1.5)[right]{${\rm \hbar(\omega_{c}+\frac{\omega_{M}}{2}-\frac{1}{2}\sqrt{4g_{M}^{2}+\omega_{M}^{2}})}$};
                 \draw[thick](1.8,3) -- (5,3)node at (8.2,3)[left]{$\small {-C_{1}\ket{1,0}+C_{2}\ket{1,1}}$}node at (-2.5,3)[right]{${\rm \hbar(\omega_{c}+\frac{\omega_{M}}{2}+\frac{1}{2}\sqrt{4g_{M}^{2}+\omega_{M}^{2}})}$};
                  \draw[thick][->](2.5,3) -- (2.5,0)node at (2.5,1.25)[left]{$\texttt{(I)}$};
                  \draw[thick][->](3.25,3) -- (3.25,0.5)node at (3.3,1.25)[left]{$\texttt{(II)}$};
                  \draw[thick][->](4.1,1.5) -- (4.1,0)node at (4.2,1.25)[left]{$\texttt{(III)}$};
                  \draw[thick][->](4.85,1.5) -- (4.85,0.5)node at (4.9,1.25)[left]{$\texttt{(IV)}$};
    \end{tikzpicture}
    }
    }
\captionsetup{
format=plain,
margin=1em,
justification=raggedright,
singlelinecheck=false
}
\caption{A energy level diagram explaining the appearance of dressed states as a result of coupling between bare optomechanical states. Different transitions are numbered corresponding to spectrum peaks in FIG.~4}
\end{figure*}
Diagonalizing the system Hamiltonian (Eq.[\ref{Hsys}]) using $\lbrace \ket{100},\ket{110},\ket{001},\ket{011}\rbrace$ as a basis, produces the set of eigenvalues:
\begin{equation}
\begin{aligned}
& \Bigg\lbrace 0,\hspace{3mm}\hbar\omega_{M},\hspace{3mm} \hbar\Bigg(\omega_{c}+\frac{\omega_{m}}{2}\pm\frac{1}{2}\sqrt{4g_{M}^{2}+\omega_{M}^{2}}\Bigg)\Bigg\rbrace.
\end{aligned}
\end{equation}
The corresponding eigenvectors $\ket{\lambda_{i}}$, $\forall$ $1\leq i\leq4$ (the so-called dressed states) take the form
\begin{subequations}
\begin{eqnarray}
\ket{\lambda_{1}}=\ket{0,0},\ket{\lambda_{2}}=\ket{0,1},\\
\ket{\lambda_{3}}=-N_{1}\ket{1,0}+N_{2}\ket{0,1},\\
\ket{\lambda_{4}}=N_{2}\ket{1,0}+N_{1}\ket{0,1},
\end{eqnarray}
\end{subequations}
with
\begin{equation}
N_{1}=\frac{2g_{M}}{\sqrt{4g_{M}^{2}+(\omega_{M}+\sqrt{4g_{M}^{2}+\omega_{M}^{2}})^{2}}},
\end{equation}
 and
 \begin{equation}
 N_{2}=\frac{(\omega_{M}+\sqrt{4g_{M}^{2}+\omega_{M}^{2}})}{\sqrt{4g_{M}^{2}+(\omega_{M}+\sqrt{4g_{M}^{2}+\omega_{M}^{2}})^{2}}}.
 \end{equation}
The energy level diagram showing the dressed states resulting from the optomechanical coupling between bare states is shown in FIG.~3. We note that because of the coupling between optical and mechanical degrees of freedom, the two bare states having a single photon combine to form two dressed states (upper two states in right part of the figure). In the dressed states picture there are  four transitions  possible among the different states. All of these four transitions have a different frequency. And these frequencies turn out to be exactly located at the peak locations in FIG.~2  worked out above by setting the real part of the poles zero.

There is another and rather simpler way of expressing the resonances. For that we can summarize the peak locations as $\Delta=\frac{g^{2}_{M}}{\omega_{M}}-m\omega_{M}$ for $m=0,1$ for both positive and negative axes of $\Delta$. This expression is just a compact form of writing the peak positions by looking at the eigenvalues of the system Hamiltonian which we stated above, and it is consistent with results reported in \cite{nunnenkamp2011single}.

\subsubsection{Asymmetry in the peak heights} We note that two of the peaks in the fully resolved spectrum are of equal height but the other two are asymmetric. Asymmetric and symmetric peaks are associated with the first and the second terms in the emission spectra of Eq.(\ref{ILTS}), respectively. We notice that mathematically this asymmetry can be attributed to the presence of the detuning parameter $\Delta$ in the numerator of first term, which under the exchange $\Delta\longleftrightarrow -\Delta$ breaks the symmetry in the heights.

 We can further explain this asymmetry by looking at the dressed state picture of the problem. Once the photon interacts with the mechanical motion it can either create a single phonon or no phonons at all. Corresponding to these two choices, the system can be either in one dressed state $\ket{\lambda_{3}}$ or in the other one, $\ket{\lambda_{4}}$. In general  we can express the state of the system as a superposition of both these dressed states. This superposition is in general imbalanced and hence the photonic transitions from such a superposed state to the two lower states ($\ket{\lambda_{1}},\ket{\lambda_{2}}$) will lead to unequal peak heights.

This explanation predicts, in fact, that in general all peaks should be of different heights. We confirmed this prediction by varying the system parameters: we found that indeed the symmetry in the two peaks of the spectrum as shown in FIG.~2 (bottom plots) is accidental and is not always present (as Figures 5 and 6 below will confirm). Moreover, we expect that, by starting with one phonon present in the initial state, these peaks should interchange their heights, as in that case the transitions between the one-phonon excited state to the one-phonon ground state will become the main process. The plot in FIG.~4 confirms this expectation.
\begin{figure}
    \includegraphics[width=3in,height=2in]{Fig4.pdf}
    \captionsetup{
  format=plain,
  margin=1em,
  justification=raggedright,
  singlelinecheck=false
}
 \caption{Effect of changing the initial phonon number in the OMC from zero to one on the long time spectrum. The parameters used are: $\kappa/\omega_{M}=0.25$, $g_{M}/\omega_{M} = 1.25 $ and $\Gamma/\omega_{M}=0.1$. Notice the interchange of peak heights with this change of initial phonon number.}
\end{figure}

Here we would also like to note that by setting $g_{M}$ and $\omega_{M}$ both equal to zero (i.e. neglecting the mechanical oscillations completely) we recover the usual Lorentzian spectrum emitted by a Fabry Per\'ot cavity (as one can expect and anticipate from the behaviour of top most red curve in FIG.~[2(b)]). We found that the general features of the spectrum (which we calculated using QJT) are consistent with the already present literature about the single photon optomechanical spectrum \cite{nunnenkamp2011single,rabl2011photon,ren2013single}. The spectra in these references were calculated mainly using Quantum Langevin Equations and input-output theory.

\subsubsection{Time-dependent spectrum} After finding the regime of parameters where one can detect the fully resolved single-photon optomechanical spectrum, next in FIG.~[5] we plot the time-dependent spectrum emitted by such a system.
\begin{figure}
\includegraphics[width=7.5cm,height=5cm]{Fig5.pdf}
\captionsetup{
  format=plain,
  margin=1em,
  justification=raggedright,
  singlelinecheck=false
}
\caption{Time-dependent spectrum emitted by OMC in the strong coupling regime with parameters $\kappa/\omega_{M}=0.25$, $g_{M}/\omega_{M} = 1.25$ and $\Gamma/\omega_{M}=0.1$}\label{Fig3}
\end{figure}
We note that with the passage of time the spectrum starts to grow, initially (until about $t\sim 6\omega^{-1}_{M}$) in the form of a broad curve shifted towards the $-\Delta$ axis. This situation corresponds to the times when the mechanical oscillator (which has zero phonons to begin with) has just started to vibrate and the photon leaked out before it could interact with the moving mirror. But at later times  we notice that side-bands start to emerge in the spectrum, indicating that now the photon has interacted with the mirror and has produced a phonon. Finally, after a long enough time ($t\sim 20\omega^{-1}_{M}$), the peaks become sharper and more pronounced, which is showing us directly how the mechanical oscillations have had an effect on the photon spectrum.
\subsection{Non-zero phonon leakage and the presence of mechanical thermal bath}
In this subsection we are going to include the losses from the mechanical oscillator by assuming that the mechanical oscillator is interacting with a finite temperature Markovian mechanical heat bath with average number of thermal phonons $\overline{M}$. This coupling opens the possibility of phonons escaping from the OMC with a rate $\gamma_{M}$. Assuming that initially the heat bath and the mechanical oscillator are in thermal equilibrium, we can specify the initial thermal state of the mechanical motion by giving the probability of finding $m$ phonons as
\begin{equation}\label{IBD}
p_{m}(t_{0})=\frac{\overline{M}^{m}}{(1+\overline{M})^{m+1}}.
\end{equation}
\begin{figure}
    \includegraphics[width=3in,height=2in]{Fig6.pdf}
    \captionsetup{
  format=plain,
  margin=1em,
  justification=raggedright,
  singlelinecheck=false
}
 \caption{Infinitely long time spectrum emitted by OMC including mechanical losses $\kappa/\omega_{M}=0.25$, $\gamma_{M}/\omega_{M}=0.1$, $\overline{M}=0.8$, $g_{M}/\omega_{M} = 1.25 $ and $\Gamma/\omega_{M}=0.1$.}
\end{figure}
The  spectrum in this case will be different from our previously calculated spectrum $P({\Delta,\infty})$ (Eq.[\ref{ILTS}]), and it will  be a weighted average
over spectra calculated from different initial
numbers of phonons.
The new parameters $\gamma_{M}$ and $\overline{M}$ will enter into the calculations when identifying the non-Hermitian Hamiltonian $\hat{H}_{dNH}$, which contains some additional terms now,
\begin{equation}\label{NHHd}
\hat{H}_{dNH}=\hat{H}_{NH}-i\hbar\frac{\gamma_{M}}{2}(\overline{M}+1)\hat{b}^{\dagger}\hat{b}-i\hbar\frac{\gamma_{M}}{2}\overline{M}\hat{b}\hat{b}^{\dagger},
\end{equation}
where $\hat{H}_{NH}$ is the same as before, as displayed in Eq.[\ref{NHH}]. Following then the same line of calculations developed in last subsection, we finally arrive at the following time-independent spectrum:
\begin{widetext}
\begin{equation}\label{ILTSD1}
\begin{split}
&P^{(1)}_{d}(\Delta,\infty)=\Bigg(\frac{1}{1+\overline{M}}\Bigg)\Bigg[\Bigg(\frac{\kappa\Gamma^{2}}{\Gamma+\overline{M}\gamma_{M}} \Bigg)\Bigg\lbrace \Bigg|\frac{i(\omega_{M}-\Delta)+\kappa/2+(\overline{M}+1)\frac{\gamma_{M}}{2}}{\lbrace i(\omega_{M}-\Delta)+\kappa/2+(\overline{M}+1)\frac{\gamma_{M}}{2}\rbrace \lbrace \kappa/2-i\Delta +\overline{M}\frac{\gamma_{M}}{2}\rbrace+g^{2}_{M}}\Bigg|^{2}\Bigg\rbrace  \\
&\hspace{17mm}+\Bigg(\frac{\kappa\Gamma^{2}}{\Gamma+(\overline{M}+1)\gamma_{M}} \Bigg)\Bigg\lbrace \Bigg|\frac{ig_{M}}{\lbrace \kappa/2-i(\omega_{M}+\Delta)+(\overline{M}+1)\frac{\gamma_{M}}{2}\rbrace \lbrace \kappa/2-i\Delta +\overline{M}\frac{\gamma_{M}}{2}\rbrace+g^{2}_{M}}\Bigg|^{2}  \Bigg\rbrace      \Bigg].
\end{split}
\end{equation}
\end{widetext}
This spectrum is shown in FIG.~[6]. For the sake of comparison we have plotted three different situations: (i) when there is no mechanical decay (red curve), (ii) when the thermal bath is at zero temperature but there is decay (black curve), and (iii) when both phonon decays and finite temperature effects are considered.

 We notice in case (ii) that there are still four peaks but the two symmetric peak heights are considerably reduced compared to other two asymmetric peaks. This fact can be explained by looking at the dressed state picture (FIG.~3). With the possibility of phonon decay, the transition between the states of one phonon to states with zero phonons will be possible. Hence the transitions (I) and (III) should be more probable now and hence peaks corresponding to these transitions (both are asymmetric) become higher than the other two symmetric peaks, which correspond to the situations in which the final state still contains a phonon. In case (iii)  all four peaks are still present and centered at the same positions, but they are now all considerably reduced in height compared to case (i).

\section{Two-phonon mechanical oscillations}
 The more strongly the photon interacts with the movable mirror the more phonons it can generate. The restriction used in the previous Section to just a single phonon won't remain a valid assumption. With this motivation in mind, we now allow the possibility of two phonons to get an idea of what aspects of the spectrum will change with the presence of additional phonons. In this Section we present results for the case when mechanical losses are included. In the QJT approach we still have the same system (as described by the non-Hermitian Hamiltonian, Eq.[\ref{NHHd}]), but the mechanical motion annihilation operator $\hat{b}$ should now be expressed  as
 \[\hat{b}=\ket{0}_{b}\bra{1}+\sqrt{2}\ket{1}_{b}\bra{2}\] with $\ket{m}_{b}$ is the phonon number state with $0\leq m\leq2$.
 The restriction to at most two phonons  implies that $\hat{b}^{\dagger}\ket{m}=0$, $\forall$ $m\geq 2$. The No-Jump state  must now account for additional possibilities of finding excitations in the system, and we write
 \begin{equation}
\begin{aligned}
&\ket{\tilde{\psi}(t)} = d_{1}(t)\ket{10}+d_{2}(t)\ket{11}+d_{3}(t)\ket{12}.
\end{aligned}
\end{equation}
Here we are using a different symbol $d$ for the amplitudes just to make the distinction with the previous Section more clear. We do use the same notational convention as before so that the second slot is reserved for displaying the number of phonons, and the first slot gives the number of photons. Assuming again that initially there was no phonon in the system, the infinitely long time spectrum now has the form
\begin{widetext}
\begin{equation}\label{ILTSD2}
\begin{split}
&P^{(2)}_{d}(\Delta,\infty)=\Bigg(\frac{1}{1+\overline{M}}\Bigg)\Bigg[\Bigg( \frac{\kappa\Gamma^{2}}{\Gamma+\overline{M}\gamma_{M}}\Bigg)\Bigg|D_{1}(s=-i\Delta)\Bigg|^{2}+\Bigg(\frac{\kappa\Gamma^{2}}{\Gamma+(3\overline{M}+1)\gamma_{M}}\Bigg)\Bigg|D_{2}(s=-i\Delta-i\omega_{M})\Bigg|^{2}  \\
&\hspace{20mm}+\Bigg(\frac{\kappa\Gamma^{2}}{\Gamma+(\overline{M}+1)\gamma_{M}}\Bigg)\Bigg|D_{2}(s=-i\Delta-2i\omega_{M})\Bigg|^{2}  \Bigg].
\end{split}
\end{equation}
\end{widetext}
The expressions of the amplitudes in Laplace space are rather involved and will not be shown here. Spectrum plots (both with mechanical losses and without losses) are shown in FIG.~[7]. We note that with the inclusion of one more phonon in the system, multiple additional side bands appear. Corresponding to the Laplace amplitude $D_{1}(s=-i\Delta)$ there are three peaks now and their location can be found by setting the real part of its pole equal to zero. The remaining peaks are then located at integer multiples of $\omega_{m}$ away from these three peaks where the integer here can either be 1 or 2.

\begin{figure}[h]
    \includegraphics[width=3in,height=2in]{Fig7.pdf}
    \captionsetup{
  format=plain,
  margin=1em,
  justification=raggedright,
  singlelinecheck=false
}
 \caption{Two-phonon infinitely long-time spectrum emitted by the OMC with mechanical losses included. Notice the appearance of many side bands when going to two phonons. Parameters are the same as in FIG.~5}
\end{figure}

In the lossless case, the analytic expression of the spectrum Eq.[\ref{ILTSD2}] predicts nine resonances, while we find eight peaks in the plot. A close inspection of the peak locations calculated from the spectrum indicates that the location of the highest peak (peak located at $\Delta = -1.25 \omega^{-1}_{M}$) occurs twice (both from $D_{1}(s=-i\Delta)$, $D_{3}(s=-i\Delta-2i\omega_{M})$). Hence the total number of peaks is one less than the naively expected number of resonances. This also explains  why this peak is the highest among all resonances.

When losses are included (red curve in FIG.~[7]) transitions among different phonon number states become more probable than before. This causes a redistribution in the peak heights, and the peaks  become wider as well. This broadening feature at finite bath temperatures also causes a dissolution of the smaller peaks.

Finally, we note that by following the same procedure of calculations our analysis of the  two-phonon scenario can be straightforwardly extended numerically to multi-phonon situations.
\section{Conclusions}
We demonstrated how a single photon interacts with a moving mirror, and how the single-photon time-dependent spectrum reveals the properties of that interaction. Resonances in the spectrum show how many phonons were created
and what the strength of the photon-mirror interaction is. A simple dressed-state picture
suffices to explain the positions and relative strengths of those resonances. The time-dependent spectrum shows how the resonances are built up in time by the photon interacting with the moving mirror and generating phonons.


text:
\section{Introduction}
Climate archives from the North Atlantic region show repeated
shifts in glacial climate, the Dansgaard-Oeschger (DO) events
\citep{Grootes1993, Andersen2006}. During Marine Isotope Stages 2 and 3 the
intervals between the events exhibit a tendency to coincide approximately with
multiples of 1470 years \citep{Rahmstorf2003}, as depicted in figure 1. The
statistical significance of this pattern and the responsible mechanism,
however, is still a matter of debate \citep{Ditlevsen2005}.
Several hypotheses were proposed to explain the timing of DO
events. One of these relates the events to two century-scale
solar cycles with periods close to 1470/7 (=210) and 1470/17
($\approx$87) years \citep{Braun2005}, the so-called De Vries/Suess
\citep{Wagner2001} and Gleissberg \citep{Peristykh2003}
cycles. Support for a leading solar role comes from deep-sea sediments,
which indicate that during the Holocene century-scale solar variability was
a main driver of multi-centennial scale climate changes in the North
Atlantic region \citep{Bond2001}.

\begin{figure}
\noindent
\includegraphics[width=20pc,height=15pc]{Figure1.eps}
 \caption{DO events as seen in two ice cores from Greenland. Top:
   NGRIP. Bottom: GISP2. Labels mark the events 1-10, the Aller{\o}d (A) and
   the end of the Younger Dryas (0). Dashed lines are spaced by
   1470 years.}
\end{figure}

Recently the phase relation between solar variability (deduced
from $^{10}$Be) and 14 DO events was analyzed \citep{Muscheler2006}. A
relation far from fixed was found and was interpreted as being in
contradiction to Braun et al.'s hypothesis \citep{Braun2005}. While in linear systems a
constant phase relation between the forcing and the response
is expected, such a relation does not necessarily exist in
non-linear systems. But climate records and ocean-atmosphere models
\citep{Ganopolski2001}, which are not yet suitable for statistical analyses on
DO events because of their large computational cost, suggest that the events
represent switches between two climate states, consequently implying an
intrinsically non-linear dynamical scenario. Thus, to interpret the reported
lack of a fixed phase relation between the DO events and solar proxies, it is
crucial to analyze their phase relation in simple models.

\section{A Simple Model of DO Events}

Here we investigate this phase relationship in a very simple model of DO
events. A comprehensive description of this model has already been published
before, including a detailed discussion of its geophysical motivation and its
applicability, as well as a
comparison with a much more detailed ocean-atmosphere model \citep{Braun2007}.
In the simple model, which was derived from the dynamics of the
events in that ocean-atmosphere model, DO events represent repeated
switches between two possible states of operation of a bistable,
excitable system with a threshold (Figure 2). These states
correspond to cold and warm periods of the North Atlantic region
during DO cycles. The switches are assumed to occur each time the
forcing function (f), which mimics the solar role in driving DO events,
crosses the threshold (T). Transitions between the two states are accompanied
by an overshooting of the threshold, after which the system relaxes back to
its respective equilibrium following a millennial-scale relaxation.

The rules for the transitions between both states are illustrated
in figure 2. It is assumed that the threshold function $T$ is
positive in the interstadial (``warm'') state and negative in the
stadial (``cold'') state. A switch from the stadial state to the
interstadial one is triggered when the forcing $f$ is smaller than
the threshold function, i.e. when $f(t)<T(t)$. The opposite switch
occurs when $f(t)>T(t)$. During the switches a discontinuity in
the threshold function is assumed, i.e. T overshoots and takes a
non-equilibrium value ($A_{0}$ during the shift into the stadial
state, $A_{1}$ during the opposite shift). Afterwards, T
approaches its new equilibrium value ($B_0$ in the stadial state,
$B_1$ in the interstadial state) following a millennial scale
relaxation process:
\begin{equation}
{\frac{dT}{dt} = -\frac{T-B_{s}}{\tau_{s}}.}
\end{equation}
Here, $\tau_0$ and $\tau_1$ represent the relaxation time in the
stadial state (s=0) and in the interstadial state (s=1), respectively.

Both the overshooting relaxation assumption and the transition
rules in our simple model are a first order approximation of the
dynamics of DO events in a coupled ocean-atmosphere model
\citep{Ganopolski2001}. In that model the events also represent
threshold-like switches  in a system with two possible states of
operation (corresponding to two fundamentally different modes of
deep water formation in the North Atlantic) and with an
overshooting in the stability of the system during these shifts
\citep{Ganopolski2001, Braun2007}. Analogous to the simple model,
switches from the stadial mode into the interstadial one are
triggered by sufficiently large negative forcing anomalies (i.e.
by a reduction in the surface freshwater flux to the North
Atlantic that exceeds a certain threshold value), whereas the
opposite shifts are triggered by sufficiently large positive
forcing anomalies (i.e. by an increase in the freshwater flux that
exceeds a certain threshold value). It has further been
demonstrated that the simple model is able to reproduce the timing
of DO events as simulated with the ocean-atmosphere model, as well
as the occurrence of non-linear resonance phenomena such as
stochastic resonance and ghost resonance, which were shown to be
properties exhibited by that model \citep{Ganopolski2002,
Braun2005, Braun2007}.

An obvious advantage of the conceptual model compared with the
ocean-atmosphere model is its low computational cost, which allows
for extensive statistical analyses on the timing of DO
events. All model-parameter values chosen here are the same as in two earlier
publications (A$_0=-27$ mSv, A$_1=27$ mSv, B$_0=-9.7$
mSv, B$_1=11.2$ mSv,
 $\tau_0=1200$ years, $\tau_1=800$ years; 1 mSv = 1
milli-Sverdrup
 = $10^{3}$ m$^{3}$/s) \citep{Braun2005, Braun2007}.

\begin{figure}
\noindent
\includegraphics[width=18pc,height=14pc]{Figure2.eps}
 \caption{Dynamics of DO events in the model. Top: Forcing f (grey)
   and threshold function T (black). Bottom: Model state s (s=0: ``cold''
   state, $s=1$: ``warm'' state). A switch from the cold to the warm
   state is triggered when $f<T$, which happens in this example at time $t''$.
   During the transition, interpreted as the start of a DO event, $T$
   overshoots and relaxes back towards its new equilibrium $B_{1}$ ($B_{1}>0$)
following a millennial time scale. The events are terminated by a switch back
to the cold state, which is triggered when
   $f>T$ (at time $t'$ in the figure). Again, T overshoots
   and approaches its new equilibrium $B_{0}$ ($B_{0}<0$).}
\end{figure}

\section{Testing Fixed Phase Relationships}

To test the assumption of a fixed phase relationship between
solar-forced DO events and solar variability, we here drive our model
by a simple input consisting of noise $\sigma \cdot n(t)$ and of two
sinusoidal cycles with equal amplitudes A:
\begin{equation}
{f(t) = -A \cdot ( \cos[\frac{2 \pi t}{T_1}] + \cos[\frac{2 \pi t}{T_2}]) +
  \sigma \cdot n(t).}
\end{equation}
 $\sigma$ is the standard deviation
of the noise and $n(t)$ the standard unit variance white noise, with a cutoff
frequency of 1/50 years (figure 3).
Following \citet{Braun2007} the cutoff is
used to account for the fact that the model shows an unrealistically large
sensitivity to decadal-scale or faster forcing.
In analogy to \citet{Braun2005} the periods of the two cycles are chosen to be $T_{1}$ =
1470/7 (=210) years and $T_{2}$ = 1470/17 ($\approx$86.5) years, i.e. close to
the leading spectral components of the solar De Vries and Gleissberg
cycles.

In the simulations shown in Figure 4 we use three
different signal-to-noise ratios (SNR): A = 8 mSv and $\sigma$ =
5.5 mSv ($SNR \approx 2.1$), A = 5 mSv and $\sigma$ = 8 mSv ($SNR
\approx 0.4$), A = 3 mSv and $\sigma$ = 9 mSv ($SNR \approx 0.1$).
For all of these, the waiting time distribution of the simulated
events is centered around a value of 1470 years, with several peaks of only
decadal-scale width (figure 4). The relative
position of these peaks is well understood in the context of the
ghost resonance theory \citep{Chialvo2002, Chialvo2003, Calvo2006,
  Braun2007}. The peaks result from constructive interference between the two
sinusoidal forcing cycles which produces particularly large magnitude
variations in the bi-sinusoidal forcing and -- when noise is added --
leads to favored transitions at the corresponding waiting times. Depending on
the relative amplitude values of the noise
and the periodic forcing this synchronization is more or less efficient,
as is seen for the different signal-to-noise ratios in Figure 4. Even for the
lowest ratio, however, the synchronization is
still notable. The waiting time distributions shown in figure 4 are thus almost
symmetrically centered around a preferred value of 1470 years because the
sinusoidal cycles enter in phase every 1470 years, creating forcing peaks of
particularly large magnitude. This 1470-year repeated
coincidence of the bi-sinusoidal forcing, however, does not show
up as a corresponding forcing frequency, since no sinusoidal cycle
with that period is present. Thus, when linear spectral analysis is performed on
the forcing, only the two century-scale sinusoidal cycles are detected as
outstanding components.

\begin{figure}
\noindent
\includegraphics[width=20pc,height=14pc]{Figure3.eps}
 \caption{Forcing and response. The forcing consists of two
  added sinusoidal cycles with equal amplitudes ($A=8$ mSv) and white noise
   ($\sigma=5.5$ mSv). (A) Total forcing f (grey) and threshold function T
   (black). (B) Forcing components
   (grey); from top to bottom: 210-year cycle, 86.5-year cycle, noise. Dashed
   lines indicate the onset of the simulated DO events. Despite the tendency
   of the three events to recur approximately every 1470
   years, only the first two events coincide with minima of the 210-year
   cycle. The third event, in contrast, occurs
   closely after a maximum of that cycle.}
\end{figure}

\begin{figure*}
\noindent
\includegraphics[width=39pc,height=18pc]{Figure4.eps}
 \caption{Waiting time distribution of the simulated events and their phase
   relation with the forcing. The amplitude of the two
   sinusoidal forcing cycles and the standard deviation of the noise
   are: A = 8 mSv, $\sigma$ = 5.5 mSv (top); A = 5 mSv, $\sigma$ = 8 mSv
   (middle); A = 3 mSv, $\sigma$ = 9 mSv (bottom). First column: normalized
   distribution of the spacing between successive DO events in the
   simulation. Second column: probability distribution of the
   phase relation between the onset of the events and the
   210-year cycle (zero corresponds to a start of the events at the
   minimum of that cycle, $\pm \pi$ to a start at the maximum). Third
   column: Rayleigh's R value \citep{Ditlevsen2007}, a measure for the
   phase correlation between the forcing cycles and the simulated events,
   as a function of the noise level $\sigma$ ($R = \frac{1}{N} | \sum_{n=1}^{N} \cos \frac{2
     \pi t_n}{T_1} + i \cdot \sin \frac{2 \pi t_n}{T_1} | $, where $t_n$
   denotes the timing of the events, $T_1$ = 210 years and $N \to
   \infty$). R is 1 if and only if a fixed phase relation exists between
   the 210-year forcing cycle and the events. Fourth column: Superposition of
   the model response (i.e. of the state variable s) and of the bi-sinusoidal
   forcing over a series of 14 simulated events, aligned by the onset of the
   events following \citet{Muscheler2006}. The
   unaveraged bi-sinusoidal forcing is normalized with maximum and minimum values of $\pm$1.}
\end{figure*}

Despite the robustness of the synchronization effect, none of the two
sinusoidal cycles in our forcing shows a fixed phase relationship with {\it{all}} of
the simulated DO events, due to the presence of noise and the existence of
a threshold. In our model, a fixed phase relation can only be present in the low
noise limit (i.e. either for $\sigma \to 0$ [with a supra-threshold
bi-sinusoidal forcing] or for the lowest noise level that still enables
repeated threshold crossings [with a sub-threshold bi-sinusoidal forcing],
thus corresponding to DO events with extremely long waiting times), compare
third column in
figure 4. Even for the largest of the three signal-to-noise ratios in our simulations, the events thus only show a
tendency to cluster around a preferred phase of the two sinusoidal
cycles. Outliers, however, can still occur in almost opposite phase, at least once over a sufficiently
large number of events in the simulation. For the highest signal-to-noise
ratio, for example, there is still a probability of about 35 percent to find
at least one out of 14 events in opposite phase (i.e. outside of the
interval $[-\pi/2,+\pi/2]$). And for the other two signal-to-noise ratios, the
corresponding probability is even much higher (i.e. 92 percent and 99.5
percent, respectively). Since a fixed phase relationship between the simulated
events and the forcing cycles does not even exist in our very simple model
system, it appears unrealistic to us to assume the existence of such a
relationship in the climate system. Thus, the reported lack of a fixed phase
relationship between DO events and solar variability (deduced from $^{10}$Be)
would also be expected with the proposed ghost resonance solar forcing. We note that
superimposed epoch analyses of 14 simulated events can produce forcing-response
relations similar to the one reported by \citet{Muscheler2006}: The onset of
the superimposed events (at $t=0$ in the fourth column in figure 4) typically
coincides with a minimum in the averaged bi-sinusoidal forcing which, however,
is not more pronounced than other minima and is highly damped as compared with
the unaveraged forcing. A considerable statistical spread exists in the
magnitude of this damping because the small number (14) of events is not yet
sufficient to infer reliable information concerning the average phase
fluctuation between the input and the output.

This lack of phase correlation between forcing and response is explained by
the threshold character of DO events: The simulated events are triggered when
the total forcing
(the sum of the two sinusoidal cycles and the noise) crosses the
threshold function. Some of the threshold-crossings are in the first place
caused by constructive interference of both cycles. These events
coincide with near-minima of the two forcing cycles. Other
threshold-crossings are, however, in the first place caused by constructive
interference of just one cycle and noise. These events thus
coincide with a near-minimum of only that cycle (compare figure 3), whereas a fixed
phase-relation with the second cycle does not necessarily exist. And at least
for low signal-to-noise ratios, some of the threshold crossings are in the first
place caused by the noise alone. These events thus do not show a fixed phase
relation with any of the two forcing cycles.

The inherently nonlinear noisy synchronization mechanism exhibited
by our model is not unique to DO events. In fact, it
has originally been proposed to explain the perception of the
pitch of complex sounds \citep{Chialvo2002, Chialvo2003} and, as a
general concept, has already been used to describe theoretically
and experimentally similar dynamics in other
excitable \citep{Calvo2006} or multi-stable systems with thresholds,
e.g. in lasers \citep{Chialvo2002, Chialvo2003, Buldu2003}.
Because of the fact that the leading output frequency is
absent in the input, this type of resonance is called ghost
stochastic resonance.

\section{Conclusions}

We here used a simple model of DO events, driven by a bi-sinusoidal forcing
plus noise, to show that a fixed phase relation between the forcing cycles
and {\it{all}} simulated events does not exist, apart from the unrealistic low
noise limit. As argued above, in this model the fluctuations in the
phases between the forcing and the response are related to the
process giving rise to the transition itself. Each event is
generated by a threshold crossing resulting from a cooperative
process between the two periodic driving forces (i.e., the
centennial-scale input cycles) and the stochastic fluctuations. In
this nonlinear scenario, as we showed explicitly in our simulations,
millennial-scale events with fixed input-output phase relations are impossible for any nonzero noise amplitude.

While one could disagree on the interpretation and the
statistical significance of the pattern described in Figure 1, our
results show that the reported lack of a fixed phase
relationship between 14 DO events and solar proxies is consistent
with the suggested solar role in synchronizing DO events. At the
same time our results have further implications for a second so-far
unexplained oscillation during Pleistocene climate, i.e. the
glacial-interglacial cycles, which also show strong indications for the
existence of threshold-like dynamics during glacial terminations
\citep{Paillard1998, Huybers2005}: Since the existence of a causal
relation between threshold-crossing events and their quasi-periodic
forcing does not necessarily imply the existence of a clear phase relation
over {\it{all}} events, the lack of such a relation between glacial
terminations on one hand and the orbital eccentricity and precession cycles on
the other hand is not sufficient to reject a leading role of these cycles
during terminations, in contrast to the interpretation proposed by
\citet{Huybers2005}. More insight in the cause of Pleistocene climate cycles
might thus be gained from more adequate statistical approaches, based e.g. on
Monte-Carlo simulations \citep{Ditlevsen2007} with simple models
\citep{Paillard1998, Braun2007} that mimic the nonlinear dynamics which seems
to be relevant during these oscillations.

\begin{acknowledgments}
H.B. thanks A. Ganopolski for instructive discussions on the dynamics of DO
events and A. Svensson for his hospitality at the Niels Bohr
Institute. H. B. was funded by Deutsche Forschungsgemeinschaft (project MA 821/33).
\end{acknowledgments}