Our objective is to formally verify the correctness of the hundreds of expression
optimization rules used within the GraalVM compiler.
When defining the semantics of a programming language,
expressions naturally form abstract syntax trees, or, terms.
However, in order to facilitate sharing of common subexpressions,
modern compilers represent expressions as term graphs.
Defining the semantics of term graphs is more complicated
than defining the semantics of their equivalent term representations.
More significantly,
defining optimizations directly on term graphs and proving semantics preservation is
considerably more complicated than on the equivalent term representations.
On terms, optimizations can be expressed as conditional term rewriting rules,
and proofs that the rewrites are semantics preserving are relatively straightforward.
In this paper, we explore an approach to using term rewrites to verify term graph transformations
of optimizations within the GraalVM compiler.
This approach significantly reduces the overall verification effort and allows for simpler encoding of optimization rules.

§ INTRODUCTION

This paper focuses on verifying the optimization transformations of expressions within the GraalVM compiler <cit.>.
Optimizations for expressions represented as terms
can be specified by conditional term rewriting rules.
However, in order to handle common sub-expressions,
the GraalVM compiler utilizes a term graph intermediate representation (IR) <cit.>
that allows common sub-expressions to be shared within the graph.
Our goal is to formally verify these expression optimizations in Isabelle/HOL <cit.>.
Our approach is to verify the corresponding term rewriting rules on the (abstract) term representation
and then give a data refinement that represents terms as term graphs
and term rewriting rules as term graph transformations.

* A framework to express and prove soundness and termination properties of optimizations expressed as term rewriting rules.

* A domain specific language implemented in Isabelle to express term rewriting rules using familiar Java syntax
and automatically generate soundness and termination proof obligations for each rule.

* Proofs for 45 optimization rules implemented in the compiler,
including sophisticated bit-twiddling optimizations using static type information.

* A general proof that term optimizations can be lifted to the term graph optimizations used by the GraalVM compiler.

* Integration of proofs using a semi-automated translation from GraalVM compiler comments.

No sharingtree1

Some sharinggraph1

Maximal sharinggraph2

Three term graphs of the expression $(x*x)+(x*x)$ with various degrees of subgraph sharinggraph-sharing

(a) Abstract syntax tree for the term (x*x)+(x*x).
(b) Term graph of (x*x)+(x*x) such that each of the two multiplication (*) nodes share a single x node for both their arguments.
(c) Term graph of (x*x)+(x*x) with a single multiplication node acting as both inputs to the addition (+) node and a single x node for both arguments of the multiplication node.

§.§ Side-Effects and Sharing in the GraalVM IR

In the GraalVM IR,
expressions are represented by
a combined control-flow and data-flow (`sea of nodes') graph <cit.>.
All components of an expression that may have side effects, such as method invocation nodes,
are part of the control-flow graph. This ensures that their evaluation takes place in the correct order.
The side-effect-free components of expressions are represented by (data-flow) term sub-graphs
each of which is a directed acyclic graph.
Internal nodes of a term graph represent operators, while
leaf nodes represent either constants, references to a parameter of the current method invocation,
or control-flow nodes with a pre-calculated value,
such as a method invocation node.

For example, the expression $3 * f(x) + 1$ is represented in the IR by
a control-flow path that first visits a method invocation node, $n$, that represents invoking $f(x)$,
plus a term graph that represents the expression $3 * n + 1$.
When the node $n$ is encountered during the control-flow execution of the graph,
$f(x)$ is invoked (possibly causing side-effects),
and the value of its result is saved in the method's state, $\methodstate$,
that is a mapping from node identifiers to values.
Any node in the subsequent control-flow traversal may evaluate the term graph corresponding to
$3 * n + 1$, in which the evaluation of the leaf node $n$ simply returns its saved value, $\methodstate(n)$.
Note that, the method invocation node $n$ plays a role in both the control-flow execution of the program
and the (side-effect-free) data-flow evaluation of any expression that includes node $n$.
As well as method invocation nodes,
nodes corresponding to accesses to fields of objects, allocating new object instances,

the equivalent of Single Static Assignment (SSA) $\phi$-nodes <cit.> act as both control-flow nodes and
as leaf nodes in term graphs representing side-effect-free expressions.

The use of term graphs to represent side-effect-free expressions in the GraalVM IR
has the advantage, compared to SSA basic blocks, of not over-specifying the sequencing of evaluation of local values.
In addition, term graphs allow one to explicitly represent sharing of common sub-expressions.
For example, the expression $(x*x) + (x*x)$ can be represented by (amongst others) any of the three term graphs
in graph-sharing.
Each figure represents a different degree of sharing of the common sub-expressions $x$ and $x*x$.
Because these term graphs all represent the same side-effect-free expression,
they are semantically equivalent,
that is, they all evaluate to the same value for any given value of $x$.
However, from the point of view of generating optimized code,
graph2 is preferred because it maximizes sharing of common sub-expressions,
and hence it is the representation typically used in the GraalVM IR.
The use of a term graph with maximal sharing
replaces the use of (global) value-numbering schemes used in
conventional SSA basic block representations, such as LLVM <cit.>.
It also allows for (global) code motion optimizations,
for example, (sub-)expressions within loops,
whose values do not change from one iteration to the next,
may be moved outside the loop and hence evaluated only once,
rather than on each iteration.

To handle an expression that may generate a runtime exception,
such as an array index out of range or a null pointer dereference,
an explicit check for the exception is inserted within the control flow graph
before the expression evaluation is reached.
Thus in optimizing expression sub-graphs, one can assume that runtime exceptions are not possible,
for example, that indices of arrays are in range and dereferenced pointers are non-null.

§.§ Expression Semantics

While the semantics of term graphs can be defined directly on the graph structure
in doing so one needs to be concerned with well-formedness issues,
such as that the term graph contains no cycles in its evaluation and
that the term graph is closed, that is,
any node identifiers referenced in the term graph are actually present within the graph.

From an expression evaluation semantics viewpoint
(but not from the viewpoint of performance of generated code)
all of the term graphs in Figures <ref>–<ref> are equivalent.
Unlike Figures <ref>–<ref>,
tree1 also corresponds to a term or abstract syntax tree representation of the expression.
In Isabelle/HOL, the term representation of an expression can be defined by a conventional (inductive) datatype,
thus allowing the semantics of terms to be defined by a conventional inductive definition over the structure of the datatype.
Any term graph can be seen as a representation of a unique term,
for example, Figures <ref>–<ref> all represent the term $(x*x)+(x*x)$,
and hence one can inductively define a relation to extract the unique term corresponding to a term graph,
thus allowing the semantics of a term graph to be defined in terms of the semantics of the corresponding term.

§.§ Representing Expression Optimizations

For side-effect-free expressions, the optimizations used in the GraalVM compiler
are primarily based on relatively simple algebraic laws, for example, the following.
One can view these laws in a  number of ways:

on values
as properties of values, for example, if $x$ and $y$ are fixed width integers
(of type $'a~word$ denoting a fixed bit width machine word from the Isabelle $Word$ library for any bit-width given by the type variable, $'a$)

on terms
as rewriting rules representing optimizations

on term graphs
as term-graph rewriting rules,

Transformation of term graph representing x-x (with node x shared) to the term graph with a single constant 0 node,
and transformation of term graph for x-(x-y) (with node x shared) to the term graph with a single y node.

The latter corresponds to the implementation in GraalVM.
An expression optimization matches a term graph with root node $r$ against the left side pattern
and, if it matches, replaces node $r$ by a new term graph for the right side,
which may contain newly created nodes (e.g. the node containing 0 in the first rule)
and/or reuse existing nodes (e.g.  node $y$ in the second rule).
Proving the correctness of such an update directly on the graph representation is non-trivial because
(a) the node $r$ being updated may be referenced from multiple places within the overall graph $g$ for the method,
(b) a graph update has the potential to introduce cycles (although none of the optimizations updates do).

The approach that we have taken is to define each canonicalization optimization
[Optimizations of this form are commonly referred to as “canonicalization rules”,
even though they do not necessarily reduce semantically equivalent expressions to the same term.]
as a term rewriting rule
using a syntactic extension defined within Isabelle,
for example,

RedundantSubtract RedundantSubtract

A rewriting rule, $\Rewrite{e_1}{e_2}{}$, is considered valid if term $e_1$ is refined by term $e_2$,
written $e_1 \refsto e_2$,
meaning that $e_2$ is well formed in any context that $e_1$ is,
and $e_1$ and $e_2$ give the same value in any such context.
For example, the proof obligation for RedundantSubtract is $x - (x - y) \refsto y$,
which is straightforward to prove based on the corresponding property for fixed width integer values redundant-sub.

The bulk of the verification work for expression optimizations is to show
that each term rewriting rule is valid with respect to the semantics of terms,
where such proofs use the corresponding property on values.
Proofs at this more abstract level are considerably simpler than
attempting the same proof directly on the term graph representation.
To promote these rules to the term graph representation,
we make use of a theorem that states that if a term graph at node $r$ representing a term $e_1$
is replaced by a term graph representing term $e_2$, where $e_1 \refsto e_2$,
then the original graph is refined by the updated graph.

term-sem gives the semantics of side-effect-free expressions on their term representations.
The semantics needs to take into account the fixed bit-width representations of integers,
for which we make extensive use of the $Word$ library in Isabelle/HOL <cit.>.
tree-trans gives expression optimizations as term rewriting rules.
graph-sem defines the term graph representation for expressions,
relates that to the term representation,
and then defines the evaluation semantics of the term graph representation via the corresponding term semantics.
application overviews the results of applying our approach to canonicalization optimizations
within the GraalVM,
and related overviews related work.

§ EXPRESSIONS AS TERMS

§.§ Representing Terms

From this point on, expression refers to a side-effect-free expression, the focus of this paper.
In program language semantics, an expression is represented as an abstract syntax tree or term.
For example, the expression $(x*x) + (x*x)$ is represented in our Isabelle/HOL <cit.> theory by the term,

which corresponds to the tree in tree1.
Terms have the advantage that they can be represented using an inductive data type,
which is simpler to work with than a term graph structure.
Terms are represented by the following Isabelle/HOL
data type, in which $IRUnaryOp$ and $IRBinaryOp$ represent
an operator (e.g. $BinNegate, BinAdd, \ldots$).
To increase the scope for optimization, nodes in the GraalVM IR have an associated $Stamp$
that contains typing and static analysis information about the node.
For integer values the stamp includes its bit-width and statically determined lower and upper bounds on its value.
Parameter and leaf expressions in our $IRExpr$ data type include a stamp;
the stamps for the other forms of expression can be calculated.
A leaf expression includes the node $ID$ of the corresponding control-flow graph node (e.g. a method call node).

Terms are defined as an Isabelle/HOL data type.
[The full Isabelle/HOL theories
[are available as an attached artifact]can be found at <https://github.com/uqcyber/veriopt-releases/tree/cpp2023>.]

datatype IRExpr 0pt  UnaryExpr IRUnaryOp IRExpr  0pt BinaryExpr IRBinaryOp IRExpr IRExpr  0pt ConditionalExpr IRExpr IRExpr IRExpr  0pt ParameterExpr nat Stamp  0pt LeafExpr ID Stamp  0pt ConstantExpr Value

When expressions are evaluated they can return many different types of values, which we model
using the $Value$ datatype.
Integers are based on the 64-bit machine word representation using the $Word$ library of Isabelle,
which provides common arithmetic and bit-manipulation operators on fixed-width words.
Integer values additionally carry the bit width $b$, where $0 < b \leq 64$, to support integer values smaller than 64 bits.
Only the lower $b$ bits of the 64-bit word are significant; the rest are 0.

datatype Value 0pt  UndefVal  0pt IntVal nat 0pt64 word0pt  0pt …

Using this $Value$ type, we define evaluation functions for unary and binary operators respectively.
For each operator, if the input values are correctly typed (for example, addition requires operands to have equal bit widths)
then the result is calculated using fixed-width operators from the $Word$ library, otherwise
$UndefVal$ is returned.

Using Isabelle/HOL's $Word$ library to handle fixed bit-width integers,
rather than (unbounded) mathematical integers, is essential
because it allows us to verify optimizations on the fixed-width values used in the compiler,
including the corner cases where fixed-width integers have subtly different semantics to mathematical integers.
It also allows the verification of optimizations like
converting a multiplication by a power of two to a left-shift operation.

§.§ Term Semantics

A term is evaluated with respect to the context of the method invocation in which it occurs.
That context includes a list, $\parameters$, of the values of the parameters to the method,
and the pre-calculated values of leaf expression components that may have had side effects, such as method invocation nodes.
The pre-calculated values are stored in a mapping, $\methodstate$,
from (control flow) node identifiers to values;
we call this map the method state.

A parameter expression is well formed if its index is within (the domain of) $\parameters$.
A leaf expression is well formed if its node identifier is in (the domain of) $\methodstate$.
Unary, binary, and conditional nodes are well formed with respect to a context $[\methodstate,\parameters]$
if their sub-expressions are well formed with respect to that context.
The relation $\MT{\methodstate}{\parameters}{e}{v}$ represents that the term $e$ is
well formed in context $[\methodstate,\parameters]$
and evaluates to value $v$ in that context.
We use a relation in Isabelle/HOL, rather than a partial function,
because it matches a conventional (big-step) operational semantics for expressions.

Note that the evaluation relation is partial because, for example,
one cannot evaluate a parameter with index $i$
if $i$ is not a valid index into the list of parameters $p$.
Evaluation assumes the term being evaluated is type correct
because the evaluation functions (e.g. bin-eval) give the result $UndefVal$
if the arguments are incompatible with the operation being applied.
The notation, $v \in s$ says that a $Value$, $v$,
is contained within the range of possible values determined by the static type information ($Stamp$), $s$.

[term semantics]
The semantics of a term is given in the form of a set of inference rules,
one for each syntactic form of expression.

Because evaluation is modeled as a relation, rather than a partial function,
we need to explicitly show that evaluation is deterministic, as given by the following lemma, which has a two-line inductive proof.

[term evaluation deterministic]
For any given context of method state $m$ and parameters $p$,
the value of a term $e$ that is well formed in that context is unique.

§ DOMAIN SPECIFIC LANGUAGE

Initial termopt1

After $\Rewrite{(x - y) + y}{x}{}$opt2

After $\Rewrite{x - (x - y)}{y}{}$opt3

Applying two rewriting rules to the tree for expression $((x-y)+y)-(x-y)$ tree-rewrites
Three abstract syntax trees (a) for ((x-y)+y)-(x-y),
(b) for x-(x-y) after rewriting (x-y)+y to y, and
(c) for y after rewriting x-(x-y) to y.

To support the large number of expression optimizations in the GraalVM compiler,
we have devised a domain specific language (DSL)
that allows one to succinctly express each optimization as a term rewriting rule
and group together a set of rules into a phase.
Our terminology for term rewriting is taken  from Baader and Nipkow <cit.>.
To show that the set of rewriting rules in a phase is terminating, a measure function is associated with each phase.
Each optimization rule is required to:

* preserve the semantics of the expression being transformed (see preservation),

* reduce the size of the expression according to the measure function (see termination).

To this end, the DSL implementation automatically generates these two proof obligations for every rule,
and to assist in the proof of these obligations,
it automatically applies a proof tactic to simplify, and in some cases discharge, the proof obligations.

A sample optimization rule expressed as a term rewriting rule within Isabelle/HOL
and the proof obligations generated from it for (1) termination and (2) refinement.
Figure <ref> gives an optimization rule called InverseLeftSub that rewrites (x-y)+y to x.
The figure also gives the generated proof obligations:
(1) the termination condition that the “size” of the rewritten term is less than size of the original, and
(2) that the rewritten term, x, is a refinement of the original term, (x-y)+y.

An optimization rule is expressed as a term rewriting rule, $\Rewrite{e_1}{e_2}{}$,
with the expression patterns written with Java-like infix notation.
optimization-dsl-example gives the definition of one such rule,
along with its two generated proof obligations.
The DSL is implemented within Isabelle/HOL using its syntax translation facilities <cit.>.
It makes use of a shallow embedding in Isabelle/HOL that allows one to leverage its built-in pattern matching,
enabling subexpressions in the pattern expression to be used in the transformed expression.

For a rule, $\Rewrite{e_1}{e_2}{}$,
upon a match of the pattern, $e_1$, with a term, $e$,
the pattern variables of $e_1$ will be matched with sub-terms of $e$.
Term $e$ is then replaced with, $e'$, the transformed expression consisting of
$e_2$ instantiated with the values of the pattern variables from matching pattern $e_1$ with $e$.
tree-rewrites shows the application of rewriting rule
$\Rewrite{(x - y) + y}{x}{}$ to the left sub-tree of opt1 to give opt2 and
then rule $\Rewrite{x - (x - y)}{y}{}$ to opt2 to give opt3,
for the example expression $((x-y)+y)-(x-y)$.

§.§ Semantics Preservation

Given that term evaluation is deterministic, refinement between terms can be defined as follows.
[term refinement]
A term $e_1$ is refined by another term $e_2$, written $e_1 \refsto e_2$,
if and only if
for all contexts $[\methodstate,\parameters]$ for which $e_1$ is well formed,
$e_2$ is well formed and evaluates to the same value as $e_1$.

Note that $e_2$ must be defined in all contexts in which $e_1$ is defined
but it may be defined in more contexts than $e_1$,
e.g. the refinement $x-x \refsto 0$ holds but $x-x$ is only defined in contexts in which $x$ is defined,
whereas the constant $0$ is defined in any context.
Refinement only considers whether $e_2$ is a semantically valid replacement for $e_1$,
not whether $e_2$ outperforms $e_1$ at run time.

[refinement monotone]
Refinement of any sub-term of a term gives a refinement of the whole term.

Rule for normalizing constant expressions to the right hand side of addition operators, along with its generated proof obligations.conditional-rewrite
Figure <ref> gives a conditional rewriting rule for swapping a constant left argument of an addition node to the right,
provided the right argument is not also a constant.
The figure also shows the generated proof obligations:
(1) the termination condition that the “size” of the rewritten term is less than the original, if the condition is satisfied, and
(2) that the rewritten term is a refinement of the original, if the condition holds.

For an unconditional term rewriting rule, $\Rewrite{e_1}{e_2}{}$, we generate the refinement obligation, $e_1 \refsto e_2$.
Conditional rewriting rules of the form, $\Rewrite{e_1}{e_2}{cond}$, restrict the rewriting to
be applicable only to terms satisfying $cond$, and hence
the conditional term rewriting rule is valid if,
\(
cond \implies (e_1 \refsto e_2) .
\)
For example, rule AddCommuteConstantRight in conditional-rewrite normalizes expressions
involving constants by moving constants to the right
and has a constraint that the second argument of the plus is not already a constant.
The abbreviation $const~c$ in the optimization syntax stands for $ConstantExpr~c$.
The free pattern variables in $cond$ are instantiated with the values of the pattern variables from the match of $e_1$ with $e$.

§.§ Optimization Phases

A set of term rewriting rules is called a term rewriting system or an optimization phase.
The rules in the set can be repeatedly applied to optimize a term, and we are
interested in various properties of the optimization phase as a whole.

Within the Isabelle/HOL representation, sets of term rewriting rules are organized into phases,
each associated with the kind of the root node of its pattern, (e.g. $AddNode$, $AbsNode$, $MulNode$, ...).

Illustrative subset of conditional canonicalization rulesconditional-rules

conditional-rules gives a sample optimization phase for conditional expressions.
For the last optimization rule, LessCond,
the condition involves the stamps of expressions $u$ and $v$;
the stamps encode the results of static analyses of the expression, including range bounds,
that give upper and lower bounds on the possible values of $u$ and $v$,
allowing this optimization to be applied if the static analysis confirms that $u$ must be less than $v$
($stamp\isacharunderscore under~a~b$ is true when the upper bound of $a$ is less than the lower bound of $b$).

To optimize an expression it is sufficient to (recursively) first optimize each of the arguments of the root node
and then apply the optimization phase for the root node (i.e. a post-fix traversal).
Note that refinement monotone ensures that optimizing the arguments gives a refinement of the whole expression.

§.§ Termination

A term $e$ is in normal form (or is irreducible) if no rewrite rules are applicable to $e$ <cit.>.
A set of rewriting rules is terminating if there does not exist a term $e$
such that there is an infinite sequence of valid applications of the rewriting rules starting with $e$.
When devising a set of rewriting rules that will be repeatedly applied to a term,
it is important that the set of rules is terminating.
For example, a rewriting rule like

\begin{align}
\Rewrite{x + y}{y + x}{}  \label{opt-plus-comm}
\end{align}

(while semantics preserving) can lead to an infinite sequence of rewrites,
e.g. $x + y \RewriteOp y + x \RewriteOp x + y \RewriteOp \cdots$.
This is undesirable because optimization phases must terminate,
and hence this rule cannot be included in a set of optimizations.
In contrast, the rule in conditional-rewrite ensures termination because of its side condition.

We restrict each phase so that its set of rewriting rules must be terminating.
Termination for a set of rewrite rules can be shown by devising a well-founded measure on terms
(a variant function)
that strictly decreases on each application of a rewriting rule.
Phases are parameterized by a measure function, $trm$, of type, $IRExpr \fun nat$.
For each rewriting rule, $\Rewrite{e_1}{e_2}{cond}$, within a phase,
one must show that the size of the term strictly decreases if its condition holds,
i.e. $cond \implies (trm~e_1 > trm~e_2)$.
This is the first proof obligation generated by the optimization command.

The measure below is satisfactory for the canonicalizations which have been encoded in this paper;
the measure is more complicated than usual in order to handle the asymmetric commuting rule in conditional-rewrite.

Note that to handle the rule in conditional-rewrite,
a lower weight is given to binary expressions with a constant on the right,

This ensures that
$\measure(const~v + y) = \measure(const~v) + \measure(y) + 2 = 1 + \measure(y) + 2 >
\measure(y) + 2 = \measure(y + const~v)$.

Phases must be parameterized by a measure as there does not exist a general measure that shows termination across all phases.
For example, multiplication by two to the power of a constant is canonicalized to left shift by the constant,
allowing more efficient generated code.
A later reassociation phase groups constants of associative operations together so that they may be constant folded.
The reassociation phase undoes the multiplication canonicalization as left shift operations are not associative and therefore limit potential reassociation.
The combined application of rules in these phases will not terminate,
instead termination relies on the sequential application of phases.

§ EXPRESSIONS AS TERM GRAPHS

The optimization phase of the GraalVM compiler represents a side-effect-free expression as a term graph.
The main reason for using a term graph representation is to explicitly share common subexpressions
in order to generate code for a common subexpression only once.
From the point of view of the semantics of an expression, its term representation
and all its term graph forms are equivalent
(i.e. evaluation in any context $[\methodstate,\parameters]$ gives the same value)
but from the point of view of generating optimized code,
the graph with maximal sharing is preferred.
Using term graphs also improves the performance of the optimizer itself
because an optimization transformation is applied only once to a shared sub-graph,
rather than to every occurrence of the common sub-expression in a term representation

§.§ IR Graph Representation

Graph nodes of the GraalVM IR are represented by the type $IRNode$,
which uses a separate node constructor for each operator corresponding to the compiler implementation, a subset of this datatype is given below where $ID$ is a numeric node identifier —
[similar to Webb et al. <cit.> who formalize the GraalVM IR Graph semantics]see the Isabelle/HOL theories[<https://github.com/uqcyber/veriopt-releases/tree/cpp2023>] or <cit.> for details.

datatype IRNode 0pt  AbsNode ID  0pt AddNode ID ID  …  0pt MulNode ID ID  0pt NegateNode ID…

An IR graph, $g$, for a method is represented in Isabelle/HOL by a finite partial function from
node identifiers ($ID$) to pairs of $IRNode$ and $Stamp$;
the latter containing static type information about the node.
Finiteness is required to allow one to use induction over the nodes of the graph.
[graph representation]
The graph for a method is represented by a finite mapping from node identifiers ($ID$) to nodes and their stamps.

For a graph to be well formed any node identifiers referenced within nodes
must be in the graph (in $dom~g$) and
the graph must not lead to cycles in evaluation of a term sub-graph.
Termination of the semantic evaluation function on the term representation is straightforward,
but not so for graphs because they may contain cycles.
[In the GraalVM IR, the graph for a term associated with a $\phi$-node may
contain a cycle back to the $\phi$-node
but during evaluation of that term, the reference to the $\phi$-node
is treated as a leaf node, thus avoiding a cycle in evaluation.]
Both these constraints are implicit in the term representation
and hence rather than explicitly constraining an $IRGraph$ to satisfy these constraints,
it is sufficient to insist that each term sub-graph represents some unique term.

§.§ Relating Terms and Term Graphs

For a given term $e$
there may be multiple term graph representations with different degrees of sharing of common subexpressions.
However, for a given a term graph, there is a unique corresponding term
and hence the relation between the term graph and term representations can be defined as
a function from term graphs to terms.
The notation $\Kind{g}{n}$ stands for the $IRNode$ with id $n$ in graph $g$,
or $NoNode$ if $n$ is not in $g$,
and the notation $stamp~g~n$ extracts the stamp of the node with id $n$ within $g$.

[term graph to term]
The ternary relation $\GT{g}{n}{e}$ denotes
that the term sub-graph of $g$ at $n$ represents the term $e$,
and is defined inductively by a set of inference rules, from which we give an illustrative subset below.
For rule (<ref>), a node $\Kind{g}{n}$ is pre-evaluated if it is of a kind that can be used as a leaf node of a dataflow expression,
such as a method call node.

[deterministic representation]
If the term sub-graph of $g$ at $n$ corresponds to a term $e$,
that term is unique.

If the term sub-graph of $g$ at $n$ corresponds to a term $e$,
that sub-graph must be closed
(because for a reference to node id $n$ not in graph, $\Kind{g}{n} = NoNode$,
which does not correspond to any term),
and not lead to a cycle in evaluation
(because the corresponding term cannot).
Furthermore, the sub-graph is well formed with respect to a context, if its corresponding term is.

[term to term sub-graph]
We use the notation $\Unrep{g}{e}{g'}{n}$ to represent an operation that updates graph $g$ to give graph $g'$,
in which node $n$ represents the term $e$.
If $e$ is already represented by some node $n$ in $g$, $g'$ is the same as $g$,
otherwise nodes are added to $g$ to give $g'$ with node $n$ being the root node of the sub-graph representing $e$.
Reusing an existing node that represents $e$, if there is one, ensures the graph preserves maximal sharing.

This relation is defined inductively by a set of inference rules, a subset of these rules are given below.
The function find-matching optionally returns any node that matches both the node type and stamp in the graph.
Rule <ref> illustrates the case where the inserted node can be represented by an existing node,
the other inference rules for these cases are similar and therefore have been omitted.
Rules <ref>-<ref> define inserting new nodes for constants, parameters, unary, binary, and leaf expressions respectively.

[term graph reconstruction]
For any insertion of a term graph representing a term, $e$, into graph, $g$,
the resulting root node identifier, $n$, in the updated graph, $g'$,
faithfully represents the original term $e$.

§.§ Term Graph Semantics

The semantics of a term graph representation of an expression is given simply
by the semantics for its corresponding term.
[term graph semantics]
Within a graph, $g$,
the semantics of a term sub-graph with root node $n$ that represents a term $e$
that is well formed for the context pair $[\methodstate,\parameters]$
is given by the semantics of the corresponding term $e$.

[term graph semantics deterministic]
If the semantics of a term graph is defined, it is unique.

The relation $\GT{g}{n}{e}$ states that the term graph at $n$ in $g$ (exactly) represents the term $e$,
but any term graph that represents a term $e'$ that refines $e$ is also acceptable,
and hence we define the relation $\GTR{g}{n}{e}$ that allows the term at $n$ to represent any refinement of $e$
(including $e$ because refinement is reflexive).
[term graph refines term]
For graph $g$, node identifier $n$ and expression tree $e$,

[term graph evaluation]
If the sub-graph of $g$ at node $n$ represents an expression that refines $e$,
evaluating the graph gives the same value as evaluating $e$ in the same context.

[term graph refinement]
A graph $g_1$ is refined by a graph $g_2$ if,
all node identifiers in $g_1$ are also in $g_2$, and
if a node $n$ in $g_1$ represents a term $e_1$,
it is refined by the corresponding term represented at $n$ in $g_2$.

A method graph $g$ also contains control-flow nodes,
which may reference term graph node ids representing expressions used by the control flow node.
[]Our earlier work gives the semantics of control-flow nodes, including method invocation <cit.>.
In all cases, the semantics of a control-flow node depends solely on
the value of a term graph in a particular context $[m,p]$,
and not on the form of the term graph.
Hence the semantics of a control-flow node is preserved if a sub-graph representing a term $e$ it uses
is replaced by a sub-graph representing a refinement of $e$.
Detailed treatment of the semantics of control-flow nodes is given in <cit.>.

Initial graph, $((x - y) + y) - (x -y)$opt1g

Graph $x - (x - y)$ after
applying $\Rewrite{(x - y) + y}{x}{}$opt2g

Graph $y$ after
applying $\Rewrite{x - (x - y)}{y}{}$opt3g

Maximal sharing graph for expression $((x-y)+y)-(x-y)$ with two successive rewrites.
The blue subgraph in (a) is rewritten to node 4 in (b).
The colored solid/dashed/dotted incoming arcs to nodes 1, 2 and 3 remain unchanged throughout (despite moving their positions in the figures).
The dashed (purple) nodes in (b) and (c) are reference (indirection) nodes.graph-rewrites
Three term graphs at various stages of transformation:
(a) term graph for ((x - y) + y) - (x -y) with sharing of the subgraph for (x-y) and of the subgraph for just y,
(b) after replacing the subgraph for (x - y) + y by x to give the term graph for x - (x - y)
but with a reference (indirection) node to the node for x as the the left argument of the root node, and
(c) after replacing the subgraph for x-(x-y) by y the original root node is now a reference node to the node for y.

§.§ Term Graph Rewriting

graph-rewrites gives the equivalent of the term rewriting example given in tree-rewrites
for the term graph representation.
In opt1g, node 3, representing the term $x - y$,
is shared as the left argument of node 2 and the right argument of node 1.
In addition,
the wider graph containing opt1g may have multiple references to any of the nodes in opt1g.
These possible references are indicated by
the solid (red) arrow for node 1,
the dotted (green) arrow for node 2,
and the dashed (purple) arrow for node 3;
nodes 4 and 5 may also have additional references
but these have been omitted from the figure to avoid clutter.
When a term graph with root node $r$ is rewritten,
it is replaced by either a new node or some existing node $n$.
This could be handled by finding every reference to $r$ within the method graph
and replacing each with a reference to $n$
(as is done in the GraalVM implementation)
but with our partial mapping representation of the graph it is simpler to
replace the node at $r$ by a special indirection node, $RefNode~n$,
that redirects any reference to $r$ to node $n$.
Indirection nodes are shown as dashed (purple) circles within graph-rewrites.
The first rewrite from opt1g to opt2g
is applied to the term graph at node 2.
It rewrites $(x-y)+y$ to $x$ by updating node 2 to be an indirection node to node 4 ($x$).
The second rewrite from opt2g to opt3g is applied to node 1.
It rewrites $x-(x-y)$ to $y$ by replacing node 1 with an indirection node to node 5 ($y$).
Note that nodes other than $y$ remain in opt3g
because they may be referenced elsewhere in the larger graph.

An optimization rule, $\Rewrite{e_1}{e_2}{cond}$, succeeds at a node $n$ within $g_1$,
if the term graph at $n$ represents a (ground) term, $e_1'$, that matches $e_1$,
and $cond$ holds for that match, otherwise it fails.
If it succeeds, a term graph is constructed that represents a term $e_2'$
that is an instantiation of $e_2$ corresponding to the match.
The new term graph replaces the node $n$ in $g_1$ to give the new graph $g_2$.
All nodes within $g_1$ other than $n$ remain unchanged within $g_2$.
The new term graph for $e_2'$ may create new nodes (not in $g_1$) and/or
reuse existing nodes within $g_1$.
If the optimization rule is valid, (i.e. $cond \implies e_1 \refsto e_2$),
for a successful application of the rule that rewrites $e_1'$ to $e_2'$, it follows that $e_1' \refsto e_2'$.
The following theorem promotes term refinement to a term graph refinement between $g_1$ and $g_2$,
where for a set of node identifiers $nids$, $nids \ndres g$ stands for the graph $g$
with all nodes in the set $nids$ removed.
[term graph semantics preservation]
Refining a term graph at some node $n$ within a graph, gives a refinement of the graph.

Note that, except at $n$, the initial graph $g_1$ remains a sub-graph of the resulting graph $g_2$,
even if nodes reachable from $n$ within $g_1$ are no longer reachable from $n$ within $g_2$.
That is because such nodes may have been shared by some other expression within $g_1$
with a different root node, and hence may still be reachable.
If they are not used as part of some other expression in $g_1$,
they will become unreachable from the start node of $g_2$ but that does not affect the semantics.
Because the node being updated may be (indirectly) referenced from anywhere in the graph,
the proof of Theorem <ref> requires an induction over all possible data-flow node types,
making it quite lengthy.
[In prior work on verifying GraalVM optimizations ]In our earlier approach <cit.>,
such a (lengthy) proof was required for each optimization rule.
Our approach []within this paper requires only the one proof to show that any
optimization that is a refinement preserves the graph's semantics.

The optimization of expressions within the GraalVM compiler takes place in two ways:

when a term graph is constructed as part of building the initial IR representation of a method,

when a term graph rewriting rule is applied to an existing sub-graph to optimize it.

In both cases a new term graph is constructed and then it is either
(<ref>) added to the graph, or
(<ref>) replaces an (unoptimized) term graph.
[term-graph construction]
If a term graph representing $e_1$ is required to be added to a graph,
a term graph representing a refinement $e_2$ of $e_1$ may be added instead.

[refined term insertion]
If a term, $e_2$, refining, $e_1$, is inserted into a graph, $g_1$, resulting in graph, $g_2$,
then the node $n'$ in $g_2$ represents $e_1$ and $g_2$ is a refinement of $g_1$.

§ APPLICATION TO GRAALVM OPTIMIZATIONS

We are applying this approach to the verification of the hundreds of expression optimization
rules used within the GraalVM compiler.
As such, it is beneficial to automatically discharge as many proofs as is practical.
We utilize the Eisbach <cit.> tactic language to develop tactics that are automatically applied to any encoded optimizations.
These tactics commonly discharge optimization termination obligations and can occasionally eliminate semantic preservation obligations.
The tactics are continually evolving as common proof patterns are identified and included in the rule set.

As an example we consider the suite of optimization rules for
conditional expressions in conditional-rules.
In the GraalVM compiler, such a set of canonicalization rules is expressed as Java code that
matches and then rewrites the IR graph to apply the rules in a defined sequence.
[Prior GraalVM verification research]Our previous approach <cit.> involved manually translating this Java code into Isabelle IR
graph transformation rules, and then proving directly that each rule preserves the IR graph semantics.
This involved a significant amount of tedious reasoning about the graph properties and semantics.
For example, verifying the set of rules for conditional expressions took a total of
58 lines of Isabelle proof commands.

With the more abstract term-based approach described in this paper, optimization rules can be specified
directly in Isabelle/HOL as term rewrite rules, as shown in conditional-rules.
Each of these rules generates two proof obligations:

* a correctness (refinement) obligation — with automated tactics disabled,
this takes two or three proof steps for each rule,
with an average of 2.4 lines/rule.  For the rules in conditional-rules,
this is just 12 lines of proof,
which is nearly five times smaller than the 58 lines required with the approach in <cit.>.
When automated tactics are enabled, all obligations without side-conditions are eliminated and
obligations with side-conditions are reduced to one or two proof steps.
This results in a total of 3 lines of proof for all 5 optimizations.

* a termination obligation — with automated tactics disabled,
this is easily proven in a single line by unfolding the
measure function then applying the simp or auto tactic.
When automated tactics are enabled, all termination obligations are eliminated.

The specifications of each optimization rule are more concise using the new approach —
typically one line for each rule (though in conditional-rules we have
added a line break before the side-conditions, for clarity), compared to an average of three lines per rule
for the term-graph approach.  As one would expect, this is also more concise than the original Java
code that implements the optimizations, which takes an average of five lines of Java code for each of these
conditional expression optimizations.

Another significant advantage is that our []new rewrite rules are more general.
The []previous graph-rewrite optimization rules [ ]<cit.> required several restrictive assumptions
(well-formed graphs, well-formed stamps, and well formed values) that are no longer required,
because those well-formed properties can be derived from the semantics of term rewrites.
Additionally, [<cit.>]previously it was assumed that the optimized term would maintain the expression evaluation relation,
i.e. the proofs showed that if the unoptimized and optimized terms evaluated to a value then that value was the same.
Now the definition of term refinement ensures that the optimized term can be used in all situations where the unoptimized term was used.
And of course, our new version also proves termination, which [was not addressed in <cit.>]our previous graph-rewrite approach <cit.> did not address.

In summary, our term-rewriting approach:

* enables more concise specification of optimization rules,

* reduces the proof burden significantly,
while proving a stronger property, and

* enables termination to be proved as well.

§.§ Workflow Integration

Ideally, we want our verification of expression optimizations to be integrated into the development workflow of the GraalVM compiler
so that new and modified optimizations continue to be verified.
We do not want the compiler developers to have to learn Isabelle or become proof experts,
but we do want a close integration between the DSL rewrite rules and the Java code that implements those rewrites.

Our approach for verifying optimizations from the compiler is to annotate the GraalVM Java source code with inline comments that utilize our optimization DSL syntax.
We also include a keyword () to identify these comment lines, and a name for each rewrite rule.
For example, the annotated GraalVM compiler code for most of the optimization rules in conditional-rules are shown in conditional-code.
Each return statement corresponding to a rewrite is preceded by a comment describing the rule.

if (condition instanceof LogicNegationNode)
LogicNegationNode negated = (LogicNegationNode) condition;
// veriopt: NegateCond: ((!c) ? t : f) |-> (c ? f : t)
return ConditionalNode.create(negated.getValue(), falseValue, trueValue, view);

if (condition instanceof LogicConstantNode)
LogicConstantNode c = (LogicConstantNode) condition;
if (c.getValue())
// veriopt: TrueCond: (true ? t : f) |-> t
return trueValue;
// veriopt: TrueCond: (false ? t : f) |-> f
return falseValue;

if (trueValue == falseValue)
// veriopt: BranchEqual: (c ? x : x) |-> x
return trueValue;

GraalVM compiler code to optimize a .
The code performs most optimization rules from conditional-rules,
the final rule is omitted for brevity.
Each optimization rule is annotated with an inline comment which explains the rule and can be used to generate proof obligations.

Fortunately the GraalVM compiler already includes similar unstructured annotations to informally describe some optimizations, so this is an extension of existing practice.
Our more precise annotations use Java-like syntax and are easily understood by developers so they are useful documentation of what transformation the Java code is performing, as well as being more concise than the corresponding Java code.
A canonicalization method in the GraalVM compiler typically handles multiple optimizations,
each corresponding to a rewriting rule.
Each optimization corresponds to a non-null return statement
[Null returns indicate that no optimization rule for that node matches.]
that returns the rewritten term graph for the corresponding rewrite rule,
and hence our annotation comments are placed before each such return statement.
This is a useful style guideline for the optimization phases of the compiler.

We apply an extraction tool to the Java source code files to extract these special comment lines, and group them into Isabelle phases.
An appropriate termination measure is manually created or re-used for each phase.
Isabelle then generates two proof obligations for each rule to (1) ensure that the rule preserves the program semantics and (2) that the set of rules is terminating.
Tactics are applied to these obligations to attempt to automatically discharge them or else simplify the obligation.
Any remaining proof obligations must be interactively discharged by an Isabelle proof expert.

The extraction tool also records statistics about the number of optimizations in each Java source file.
This semi-automated process allows us to measure progress on verifying optimization rules.
It would be useful in the future to extend the tool to check the style guideline mentioned above,
and possibly add some fully-automated testing support that searches for counter-examples
(similar to the Isabelle `nitpick' command)
to the rewrite rules before they are sent to Isabelle for full verification.
Such tooling (independent of Isabelle) will be important for these annotation comments to be accepted into the standard development lifecycle of the compiler.

§.§ Validation

Differential testing work as part of this project <cit.>
provides confidence that the IR semantics and the specification of optimizations correspond to that of the existing GraalVM compiler.
Java programs included in the compilers existing unit tests are translated into Isabelle and executed with our executable semantics.
The output of the programs executed in Isabelle are compared with the output computed through compilation and execution with GraalVM.
Likewise, we perform optimizations of programs using GraalVM and our Isabelle specification and validate that the same optimized program is produced by both approaches.

§ RELATED WORK

Representing expressions as term graphs has a long history
[Doug McIlroy communicated privately that the earliest overt exploitation of the model that he happens to know is
Vyssotsky's use of it for Fortran data-flow diagnostics in 1962.]
The GraalVM compiler makes use of this representation as part of a larger graph combining
control-flow and data-flow as a “sea of nodes” <cit.>.

Representing program transformations as conditional term rewriting rules also has a long history
They have been used in language development environments for many different purposes, including
optimization transformations.  For example, the Spoofax Language Workbench for developing
domain-specific languages (DSLs) provides the Stratego transformation language,
which extends a pure rewriting approach by introducing traversal combinators and programmable
rewriting strategies <cit.>.  Stratego has been used to define many kinds
of transformations, including optimizations and code generation.

Rewriting of term graphs <cit.>
has also been studied extensively, with the
TERMGRAPH conference[See <http://www.termgraph.org.uk>] running biannually since 2002.
The main motivation for term graph rewriting is to improve the efficiency of term rewriting
by avoiding the duplication of sub-trees.

Other research that has addressed the formal correctness of optimization rules includes the
Alive/LLVM project <cit.>, which uses an SMT solver to prove the correctness of peephole
optimization rules in the LLVM compiler and has found at least 64 bugs in LLVM and 13 in the SMT solver (Z3).[See <https://github.com/AliveToolkit/alive2> for a list of bugs found and fixed.]
However, the peephole optimizations are defined on sequences of assembler instructions rather than graphs, so the complexities of term graph rewriting are not required.

Wider scoped compiler verification projects, such as CompCert <cit.> and CakeML <cit.>, incorporate proofs of semantic preserving optimizations.
CompCert currently utilizes 10 intermediate languages,
including 2 control-flow graph based languages.
However, none of these intermediate languages explicitly represent shared sub-expressions,
instead a common sub-expression elimination phase transforms shared expressions into move instructions.
Expression optimizations similar to those presented in this paper,
are performed as peephole optimizations on sequences of assembler instructions as in Alive <cit.>,
avoiding the term graph complexities.
Likewise the CakeML intermediate representation does not explicitly represent shared sub-expressions.

The VellVM project <cit.> uses the Coq theorem prover to formally verify the correctness of
several LLVM optimizations on an SSA intermediate representation.
Mansky and Gunter <cit.>
and Strecker <cit.> have designed generic Isabelle/HOL frameworks for proving correctness of optimizations which we anticipate to be useful for handling control-flow graph based optimizations.
The former focuses on SSA-based optimizations and have specified and proved the correctness of an SSA construction algorithm
in their framework.
Both of those projects and VellVM performed verification directly on the
graph (or SSA) form of the program and were primarily concerned with control-flow optimizations,
whereas our focus in this paper is making it easier to verify the correctness of the multitude of expression optimization rules.

§ CONCLUSIONS

Our primary concern is formal proof of the correctness of each transformation rule, based
on an underlying semantics of the IR being transformed.  These proofs are more difficult when transformations
are expressed directly as term graph rewrite rules.
Our contribution is to combine the semantic simplicity of term rewriting
(leading to simplified proof obligations) with the performance-efficiency of term graph rewriting.
This gives an effective two-level approach to verifying optimizations on term graphs in Isabelle/HOL:

* the first layer of verification proves the semantic correctness of each optimization rule as a simple term rewrite rule;

* the second layer of verification is generic, and shows that any correct tree rewrite rule can be lifted into a term graph rewrite rule that preserves the semantics of the program being optimized.

This two-layered approach is simpler than directly defining the semantics of term graphs
and directly verifying term graph transformations <cit.>
which leads to complex proofs for each rule.
Using this approach we have built a framework that allows one to more easily verify expression
canonicalization optimizations for the GraalVM IR.

At this stage,
we have proved 45 optimization rules for the 18 currently supported expression nodes.
These rules include reasonably sophisticated bit-twiddling optimizations that rely on the static type information inferred by GraalVM.
However, they only represent a small set of the approximately 120 side-effect free compiler nodes.
The focus of the project to this stage has been on developing the proof infrastructure
to a level where relatively inexperienced Isabelle/HOL users now can verify optimization rules
and we plan to have honours student projects working on extending the set of optimization rules handled.

Our approach is currently limited to expression optimizations that are side-effect free.
This places optimizations involving division and remainder expressions out of reach.
[The GraalVM IR has recently been updated to include data flow versions of divide and remainder nodes,
which we do plan to handle.]
Additionally, optimizations that perform control-flow transformations are not suitable for this approach,
currently these optimizations are performed on the graph-based representation.
We plan to utilize techniques developed by Mansky et al. <cit.>,
Strecker <cit.>, and Appel <cit.> to improve our approach to this category of optimizations.

In future work, we plan to formalize a range of strategy operators <cit.>
for combining optimization rules in various ways, including traversing subexpressions, and
lift our correctness and termination results to the strategy level.  We also plan to model
the generation of rewrite implementation code from rules and strategies, so that we can generate
Java code that could be integrated into the GraalVM compiler, with a clear connection back to the
verified optimization rules.

Mark Utting's position and Brae Webb's scholarship are both funded in part by a gift from Oracle Labs.
Thanks especially to Cristina Cifuentes, Paddy Krishnan, Andrew Craik and Gergö Barany from Oracle Labs for their helpful feedback,
and to the Oracle GraalVM compiler team for answering questions.
Thanks also to Kristian Thomassen for his work on Isabelle/HOL proofs of canonicalization rules,
Kausthubram Rajesh for research on a domain specific language for rewriting,
and Jazmin McWilliam for work on verifying canonicalization optimizations from GraalVM using the framework of this paper.


# Vacuum Radiation Pressure Fluctuations on Atoms

L. H. Ford<EMAIL_ADDRESS>Institute of Cosmology, Department of
Physics and Astronomy, Tufts University, Medford, Massachusetts 02155, USA

###### Abstract

Recent work has shown that the stress tensor components, such as energy
density or pressure, of a quantum field can be subject to large vacuum
fluctuations. The energy density or pressure must be averaged in time before
the fluctuations can be finite, and the probability of a large fluctuation
depends upon the details of the averaging and can be much larger than that
predicted by a Gaussian distribution. This paper explores vacuum radiation
pressure fluctuation on Rydberg atoms and their possible observable effects.
The excitation and de-excitation of a Rydberg atom provide an explicit model
for the time averaging of the radiation pressure, as the atomic polarizability
becomes time dependent, first increasing and then decreasing again by several
orders of magnitude. This switched polarizability can induce large vacuum
pressure fluctuations, which can in turn temporarily transfer linear momentum
to the atom and cause a recoil which might be observable.

##  I Introduction

Recent work has shown that components of the stress tensor for a quantized
field can exhibit large vacuum fluctuations FewsterFordRoman2010 ; FFR2012 ;
FF2015 ; SFF18 ; FF20 ; WSF21 . The relevant quantum operator must be averaged
in time or in space and time, and the resulting probability distribution is
sensitive to the details of this averaging, but typically decreases as an
exponential of a small fractional power. Hence the probability of a large
vacuum fluctuation can be much larger than would be predicted by the Gaussian
distribution which governs random processes. This reflects the fact that
vacuum fluctuations exhibit strong correlations. The averaging in time, or
space and time, may be viewed as arising from the details of the measurement
of the stress tensor.

Quantum stress tensor fluctuations can lead to passive fluctuations of the
gravitational field, and are hence a type of quantum gravity effect which
might play a role in the early universe FMNWW10 ; WHFN11 . Analogous effects
might be observable in condensed matter system, such a quantum density
fluctuations in a fluid WF20 .

In this paper, a specific model for the measurement of vacuum radiation
pressure on an atom will be discussed. It involves a Rydberg atom which is
excited to a high energy level by a short laser pulse, and later de-excited by
another pulse. During this time, the polarizability of the atom can increase
by many orders of magnitude before decreasing again. This time dependent
polarizability acts to average the radiation pressure operator in time HF17
resulting in potentially observable velocity fluctuations of the atom.

The outline of this paper is as follows: We first review previous results on
the fluctuations of linear quantum field operators is Sec. II.1, and then
discuss the significance of measurements in finite intervals in Sec. II.2. The
resulting probability distributions for the stress tensor operators are
reviewed in Sec. II.3. In Sec. II.4, we introduce an averaging function
depending upon two time scales which will be used later in the paper. The
dependence of the probability distribution upon the details of the averaging
is discussed in Sec. II.5, and in Sec. II.6 we review the case of radiation
pressure fluctuations on atoms. In Sec. III, our model for switching the
polarizability of an atom is presented, and estimates for both the expected
atomic speed fluctuations and the probability of large fluctuations are
presented. The results of the paper are summarized in Sec. IV.

Lorentz-Heaviside units in which $\hbar=c=1$ are used throughout the paper,
except as otherwise noted.

##  II Vacuum Fluctuations of Quantum Field Operators

###  II.1 The Case of Linear Field Operators

It is well known that the vacuum fluctuations of a linear field operator, such
as the electric field, satisfy a Gaussian probabilitiy distribution. However
the operator must first be averaged in either time or space with a test or
sampling function. The use of test functions has long been used in rigorous
approaches to quantum field theory as a formal tool to obtain well-defined
operators PCT . However, we may also view these functions as describing the
physical effects of a measuring apparatus. Let $f(t)$ denote a temporal
sampling function which vanishes as $t\rightarrow\pm\infty$ and is normalized
so that

$\int_{-\infty}^{\infty}f(t)\,dt=1\,.$ (1)

Now let $E({\bf x},t)$ be one cartesian component of the free electric field
operator, and define its time average by

${\bar{E}}=\int_{-\infty}^{\infty}f(t)\,E({\bf x},t)\,dt\,.$ (2)

Now ${\bar{E}}$ will undergo Gaussian fluctuations around a mean value of
$\langle 0|{\bar{E}}|0\rangle=0$ with a finite variance of

$\sigma=\langle 0|{\bar{E}}^{2}|0\rangle\,.$ (3)

If $\tau$ is the characteristic duration of the sampling, then
$\sigma\propto\tau^{-4}$, with the constant of proportionality dependent upon
the detailed functional form of $f(t)$. As $\tau\rightarrow 0$, the variance
grows, $\sigma\rightarrow\infty$. The mathematical reason for this is that
smaller values of $\tau$ cause higher frequency modes to give the dominant
contribution to ${\bar{E}}$. The physical reason is that measurements on
shorter time scales probe larger vacuum fluctuations. Here we have
concentrated on time averaging, but averaging in space alone or in both space
and time produce similar results for linear field operators.

The vacuum fluctuations of the electric field can have observable effects.
Welton Welton  showed that the dominant contribution to the Lamb shift may be
estimated by a simple argument involving the effects of electric field
fluctuations on energy levels of the hydrogen atom. It was argued in Ref. HF15
that vacuum electric field fluctuations will slightly increase the probability
of barrier penetration by a charged particle, and that a simple estimate using
the time averaged electric field operator agrees with a one-loop perturbation
theory calculation FZ99 .

###  II.2 Measurements in a Finite Interval

It is often convenient to select simple functions, such as Lorentzians or
Gaussians for the sampling functions. However, these have tails which strictly
correspond to a measurement which began in the infinite past and extends to
the infinite future. A better choice is a function which is strictly zero
outside of a finite interval, a function with compact support. Such a function
is necessarily non-analytic, but can be taken to have all of its derivative
finite, a $C^{\infty}$ function. Some examples of compactly supported
functions are constructed and discussed in Refs. FF2015 ; FF20 . Here we
consider such a function of time, $f(t)$, with duration of the order of
$\tau$, the sampling time.

A key feature of compactly supported functions is the rate at which their
Fourier transforms decay. Let the Fourier transform of $f(t)$ be

$\hat{f}(\omega)=\int_{-\infty}^{\infty}dt\,{\rm e}^{-i\omega t}\,f(t)\,.$ (4)

If $f(t)$ is a compactly supported $C^{\infty}$ function, then
$\hat{f}(\omega)$ will fall more rapidly that any power of $\omega$, but more
slowly than an exponential function as $|\omega|\rightarrow\infty$. Here we
take both $f(t)$ and $\hat{f}(\omega)$ to be real, even functions, and
consider a class of compactly supported functions for which

$\hat{f}(\omega)\sim\gamma\,{\rm
e}^{-\beta|\omega\tau|^{\eta}},\qquad|\omega|\to\infty$ (5)

for some constants $\gamma$, $\beta>0$, and $0<\eta<1$. (Note that the
parameter $\eta$ is denoted by $\alpha$ in Refs. FF2015 ; FF20 , but in this
paper we reserve $\alpha$ to denote atomic polarizability.) In all cases,
$\hat{f}(\omega)$ decreases more slowly than exponentially for large
$|\omega|$. The value of $\eta=1/2$ is of physical interest. A simple
electrical circuit is given in Ref. FF2015  in which the current increases as
a function of time with $\eta=1/2$ after a switch is closed.

The value of $\eta$ is linked to the rate of switch-on and switch-off of
$f(t)$; the more sudden the switching, the smaller will be the value of
$\eta$, and the more important will be the contribution of high frequency
modes. This contribution will play an important role in the probability of
large vacuum fluctuations for quadratic operators.

###  II.3 Probability Distributions for Stress Tensor Operators

The case of operators which are quadratic in the fields, such the energy
density and other components of the stress tensor is more complicated, and has
been treated in several recent papers FewsterFordRoman2010 ; FFR2012 ; FF2015
; SFF18 ; FF20 . Here the operator must be averaged in time in order to have
finite fluctuations. Additional averaging in space is optional and is
discussed in Ref. FF20 . Let $\bar{T}$ be the time average of a normal ordered
quadratic operator $T(t)$, defined as in Eq. (2). This averaged quadratic
operator may be expanded in normal modes as

$\bar{T}=\sum_{i\,j}(A_{ij}\,a^{\dagger}_{i}\,a_{j}+B_{ij}\,a_{i}\,a_{j}+B^{*}_{ij}\,a^{\dagger}_{i}\,a^{\dagger}_{j})\,.$
(6)

Here $a^{\dagger}_{i}$ and $a_{i}$ are the creation and annihilation operators
for mode $i$ of a bosonic field. Quantization in a finite volume is assumed,
so the modes are discrete. The $A_{ij}$ and $B_{ij}$ are symmetric matrices
which are proportional to $\hat{f}$. Specifically, we expect

$A_{ij}\propto(\omega_{i}\,\omega_{j})^{p/2-1}\;\hat{f}(\omega_{i}-\omega_{j})\,,\qquad
B_{ij}\propto(\omega_{i}\,\omega_{j})^{p/2-1}\;\hat{f}(\omega_{i}+\omega_{j})\,,$
(7)

where $p$ is an odd integer which depends upon the dimensions of the operator
$\bar{T}$. For example, for $T=\varphi^{2}$, the square of a massless scalar
field , $p=1$, For $T=\dot{\varphi}^{2}$, or a typical component of the stress
tensor for a massless field, $p=3$. In Sec. II.6, we will consider operators
with higher time derivatives, and hence larger values of $p$.

One approach to finding the probability distribution for the vacuum
fluctuations of $\bar{T}$ is to study the behavior of the moments,

$\mu_{n}=\langle\bar{T}^{n}\rangle\,,$ (8)

where the expectation value is in the vacuum state. Because $\bar{T}$ is
normal ordered, the first moment vanishes, $\mu_{1}=0$. The probability
distribution has equal area on either side of $\bar{T}=0$, so a given
measurement is equally likely to return a negative value as a positive value,
even for quantities such as energy density which are non-negative in classical
physics. This does not mean that the distribution is symmetric. For quantities
such as the energy density, the odd moments are non-zero and the distribution
is skewed. The width of the distribution is described by the second moment, or
variance,

$\mu_{2}=2\sum_{i\,j}B_{ij}\,B^{*}_{ji}\,,$ (9)

The fact that $\omega_{i}+\omega_{j}\geq 0$, combined with the decrease of
$\hat{f}$ for increasing values of its argument, guarantees that the variance
is finite so long as a time average over a nonzero interval has been
performed.

We can also understand now why time averaging is essential for quadratic
operators. Space averaging of $T$ with a spatial sampling function $g({\bf
x})$ will introduce a factor of its spatial Fourier transform, $\hat{g}({\bf
k})$, into each of $A_{ij}$ and $B_{ij}$, if we use a plane wave mode basis.
In particular, now $B_{ij}\propto\hat{g}({\bf k_{i}}+{\bf k_{j}})$. However,
this factor is not suppressed for pairs of modes where ${\bf k_{i}}=-{\bf
k_{j}}$. If there is no time averaging, and hence no factors of
$\hat{f}(\omega_{i}+\omega_{j})$, then $\mu_{2}$ receives a divergent
contribution from such modes.

The dimensionless measure of the value of the averaged operator $\bar{T}$ is

$x=\tau^{p+1}\,\bar{T}\,.$ (10)

We seek the probability distribution $P(x)$, and especially its asymptotic
form for $x\gg 1$. This form is determined by the behavior of the moments
$\mu_{n}$ for large $n$. The general moment $\mu_{n}$ may be expressed as a
mode sum of an $n$-th degree polynomial in the $A_{ij}$ and $B_{ij}$, and
hence in the Fourier transform, $\hat{f}$. If the asymptotic form in Eq. (5)
holds, and there is no spatial averaging then FF2015

$\mu_{n}\propto(p\,n/\eta)!$ (11)

for large $n$. This rapid growth of the moments leads to a slow decrease in
$P(x)$ for large $x$. Specifically,

$P(x)\sim c_{0}\,x^{b}\,{\rm e}^{-ax^{c}}$ (12)

where the constants $c_{0}$, $b$, $a$, and $c$ depend upon the choice of
sampling function. In particular, with only time averaging, we have

$c=\frac{\eta}{p}$ (13)

and

$b=\frac{2-\eta}{p}-(\eta+1)\,.$ (14)

The most crucial of these constants is $c$, and its relatively small value
describes the slow rate of decrease of $P(x)$. For example, a stress tensor
component such as energy density, with $p=3$, sampled with a temporal sampling
function with $\eta=1/2$, will have $c=1/6$, implying a relatively large
probability for large fluctuations.

The vacuum is a state with zero total energy, but fluctuations in the energy
contained within finite regions of space or time are possible. We can view
these fluctuations as a temporary loan which must be returned within a finite
time, possibly through anticorrelated fluctuations. The timescale for the
return of energy to the vacuum is inversely related to the magnitude of the
energy. This is enforced by the probability distribution, Eq. (12). For a
fixed magnitude of the averaged energy density, the dimensionless variable $x$
increases with increasing $\tau$, causing $P(x)$ to decrease. Thus the
probability of finding a given value of the averaged energy density $\bar{T}$
decreases with increasing averaging time.

The case of averaging both in time and in space is treated in Ref. FF20 , and
here we summarize some of the key results. Now we have a temporal sampling
function of width $\tau$ and a spatial function $g({\bf x})$ of characteristic
radius $\ell$. Here we are interested in the case where $\ell\ll\tau$, so we
expect a regime where only the temporal sampling is significant. This arises
if $x$, still defined by Eq. (10), is not too large. If $1\ll x\lesssim
x_{*}$, where

$x_{*}=\left(\frac{2}{a}\right)^{p/\eta}\;\left(\frac{\tau}{\ell}\right)^{p}\,,$
(15)

then $P(x)$ is approximately given by Eq. (12) with $c$ given by Eq. (13).
However, for $x\gtrsim x_{*}$, $P(x)$ makes a transition to a regime with the
same functional form, but with $c$ given by

$c=\eta\,.$ (16)

Now $P(x)$ falls more rapidly than in the worldline approximation, where
spatial averaging is ignored, so the effect of spatial averaging is to
somewhat decrease the probability of large fluctuations. However, with
$\eta<1$, the decay is still slower than exponential.

The results in Refs FF2015 ; FF20 , which have just been summarized, were
obtained by a study of the behavior of the moments. There is an alternative
approach to finding $P(x)$ which involves diagonalization of the operator
$\bar{T}$ by a Bogolubov transformation. This allows the determination of the
eigenstates and eigenvalues of this operator. The probability of finding a
given eigenvalue in a measurement is the squared overlap of the corresponding
eigenstate with the vacuum, which yields $P(x)$. In practice, this procedure
needs to be performed numerically in a system with a finite number of modes.
This was done in Ref. SFF18  for the worldline case and in Ref. WSF21  for the
spacetime averaged case, and the results are in good agreement with those of
the moments approach.

###  II.4 Sampling Functions with Two Time Scales

In Sec. II.2, we discussed sampling functions of time which depend upon a
single time scale, $\tau$. Now we consider functions which depend upon two
different scales, $\tau$, which is the characteristic rise or decay time, and
a potentially much longer duration, $t_{0}$. Let $F(t)$ be a function which
varies from values of zero to one on a timescale of the order of $\tau$, so
$F(t)=0$ for $t<0$ and $F(t)\rightarrow 1$ for $t\gtrsim\tau$. An explicit
example is

$F(t)={\rm e}^{-\tau/t}\qquad t>0\,.$ (17)

We can define a two-time scale sampling function by

$f_{2}(t)=\frac{C_{2}}{t_{0}}\,F(t+t_{0}/2)\,F(t_{0}/2-t)\,.$ (18)

This function switches on at $t=-t_{0}/2$ and off again at $t=t_{0}/2$, and
hence has a total duration of $t_{0}\geq 2\tau$. The rise and fall time is
characterized by $\tau$. If $F(t)$ has the form given in Eq. (17), then $\tau$
is the time required for $F(t)$ to increase from zero to $1/{\rm e}$ of its
final value, so the total switch-on or switch-off time for $f_{2}(t)$ will be
a few times $\tau$. An example of such a function is illustrated in Fig. 1.
The switch-on region of the same function is shown in more detail in Fig. 2.

Figure 1: A two-time scale sampling function $f_{2}(t)$, with a duration of
$t_{0}$ and switch-on and switch-off times of the order of $\tau$. Here
$t_{0}=50\,\tau$, and $f_{2}(t)$ was constructed using Eqs. (17) and (18).
Figure 2: Here the switch-on of the function $f_{2}(t)$ with $t_{0}=50\,\tau$
is illustrated in detail.

We can define a single time scale sampling functions, such as those used in
Sec. II.2, by setting $t_{0}=2\,\tau$ and writing

$f_{1}(t)=\frac{C_{1}}{2\tau}\,F(t+\tau)\,F(t-\tau)\,.$ (19)

Here $C_{1}$ and $C_{2}$ are constants defined so that these functions are
normalized as in Eq. (1). Note that if $t_{0}\gg 2\tau$, then $C_{2}\approx
1$.

The Fourier transforms of these temporal functions, $\hat{F}(\omega)$,
$\hat{f_{1}}(\omega)$, and $\hat{f_{2}}(\omega)$, are defined as in Eq. (4),
and are related to one another. In particular, Eq. (18) and the convolution
theorem imply

$\hat{f_{2}}(\omega)=\frac{C_{2}}{2\pi
t_{0}}\,\int^{\infty}_{-\infty}d\omega_{1}\,\hat{F}(\omega_{1})\,\hat{F}(\omega_{1}+\omega)\,{\rm
e}^{i\omega_{1}t_{0}}\,.$ (20)

Note that the exponential factor will suppress contributions to the integral
from values of $\omega_{1}\gg 1/t_{0}$. Hence, if $\omega\gg 1/t_{0}$, we have
$\hat{F}(\omega_{1}+\omega)\approx\hat{F}(\omega)$ and

$\hat{f_{2}}(\omega)\approx\frac{C_{2}}{2\pi
t_{0}}\,\hat{F}(\omega)\,\int^{\infty}_{-\infty}d\omega_{1}\,\hat{F}(\omega_{1})\,{\rm
e}^{i\omega_{1}t_{0}}=\frac{C_{2}}{t_{0}}\,F(t_{0})\,\hat{F}(\omega)\,.$ (21)

If we set $t_{0}=2\tau$, then the above result gives the asymptotic form of
$\hat{f_{1}}(\omega)$ for large $\omega$. We see that $\hat{F}(\omega)$,
$\hat{f_{1}}(\omega)$, and $\hat{f_{2}}(\omega)$ all have the same asymptotic
functional forms, given by Eq. (5) with the same values of $\beta$ and $\eta$,
but different values of $\gamma$.

Because $f_{1}(t)$ and $f_{2}(t)$ are normalized sampling functions which
satisfy Eq. (1), Eq. (4) implies that $\hat{f_{1}}(0)=\hat{f_{2}}(0)=1$.
However, $F(t)$ is not normalizable, so we expect
$\hat{F}(\omega)\rightarrow\infty$ as $\omega\rightarrow 0$. Because $F(t)$ is
a dimensionless function, $\hat{F}(\omega)$ has dimensions of time and is
hence proportional to $\tau$. For $\omega\gg 1/\tau$, we expect

$\hat{F}(\omega)\approx\gamma_{0}\,\tau\,{\rm
e}^{-\beta|\omega\tau|^{\eta}}\,,$ (22)

where $\gamma_{0}$ is a constant of order one. If $t_{0}\gg\tau$, where
$C_{2}\approx 1$ and $F(t_{0})\approx 1$, Eq. (21) becomes

$\hat{f_{2}}(\omega)\approx\gamma_{0}\,\frac{\tau}{t_{0}}\,{\rm
e}^{-\beta|\omega\tau|^{\eta}}\,.$ (23)

###  II.5 The Parameters $\eta$ and $\beta$

The rate of decay of the Fourier transform of a sampling function for large
$\omega$ is given by Eq. (5), which effectively defines the parameters $\eta$
and $\beta$. Both of these parameters play important roles in our discussion,
as the probability distribution describing large fluctuations is sensitive to
both. The greater sensitivity is to the value of $\eta$, which appears in the
power of $\omega$ in the exponential in $\hat{f}(\omega)$, and in the power of
$x$ in $P(x)$, so smaller values of $\eta$ are linked to more slowly
decreasing $\hat{f}(\omega)$ and to a greater probability of large
fluctuations. The same is true to a lesser degree of the value of $\beta$,
which appears as a coefficient inside the exponential functions in both
$\hat{f}(\omega)$ and $P(x)$.

The value of $\eta$ is associated with the functional form of $f(t)$ near the
switch-on time, which we take here to be $t=0$. The class of compactly
supported functions, whose Fourier transforms have the asymptotic form in Eq.
(5), have the switch-on behavior

$f(t)\sim t^{-\mu}\,{\rm e}^{-w\,t^{-\nu}}\,,\quad t\rightarrow 0^{+}\,,$ (24)

where $\nu=\eta/(1-\eta)$, $w$, and $\mu$ are constants determined by $\eta$.
[See Eqs. (50)-(54) in Ref. FF2015 ]. In particular, $\eta=1/2$ is associated
with $\nu=1$, as illustrated in Eq. (17). As noted above, the functions
$F(t)$, $f_{1}(t)$, and $f_{2}(t)$, related by Eq. (18), or similar relations,
all have the same values for $\eta$ and $\beta$. Recall that the electrical
model in Ref. FF2015  creates a switch-on described by a function with
$\eta=1/2$. The value of $\beta$ can vary in this model, but is determined by
the parameters of the specific circuit.

We can better understand the role of the parameter $\beta$ from the study of a
class of functions discussed by Johnson Johnson  of the form

$f_{J}(t)=C\,{\rm e}^{-\beta^{2}\,[1-(t/\tau)^{2}]^{1-a}}\,,$ (25)

where $b>0$, $a>1$, $-\tau<t<\tau$, and $C$ is a normalization constant. Note
that the switching behavior is similar to that in Eq. (24), with $\nu=a-1$,
but without the factor of $t^{-\mu}$. Johnson evaluates the asymptotic form
for the Fourier transform $\hat{f_{j}}(\omega)$ using a saddle point
approximation, and finds for the case $a=2$ that

$\hat{f_{J}}(\omega)\propto\omega^{-3/4}\,{\rm
e}^{-\beta\,\sqrt{\tau\,\omega}}\,,\quad\omega\gg 1/\tau\,.$ (26)

This shows that $a=2$ corresponds to $\eta=1/2$, and in this case, the
parameter $\beta$ in Eq. (25) is the same as that in Eq. (5). Note that
$\hat{f_{J}}(\omega)$ contains a factor of $\omega^{-3/4}$ which does not
appear in Eq. (5), and $f_{J}(t)$ lacks the factor of $t^{-\mu}$ which appears
in Eq. (24), which show that $f_{J}(t)$ and $f(t)$ are somewhat different
functions. Nonetheless, they share the same exponential factor in their
Fourier transforms.

Here we assume that the $\beta$-dependence of $f_{J}(t)$ is a reasonable guide
to understanding that of $f(t)$. In Fig. 3, $f_{J}(t)$ is plotted for the case
$a=2$ for three values of $\beta$. We can see that smaller $\beta$ leads to
more rapid switch-on and switch-off, as one might expect from the increased
contribution of high frequency Fourier components as $\beta$ decreases.

Figure 3: The normalized function $f_{J}(t)$ is plotted for the case $a=2$, or
$\eta=1/2$, for three values of $\beta$

###  II.6 Vacuum Radiation Pressure Fluctuations on Atoms: Probability
Distribution

Vacuum radiation pressure fluctuations on a perfect mirror could in principle
be computed from those of the appropriate components of the electromagnetic
stress tensor. However, in the present paper we are concerned with atoms,
which scatter longer wavelengths of light by Rayleigh scattering. The vacuum
radiation pressure fluctuations on such polarizable particles were treated in
Sec. V of Ref. HF17 , and the results will be summarized here. The Rayleigh
cross section for the scattering of light with angular frequency $\omega$ from
a particle with a static polarizability $\alpha$ is

$\sigma_{R}=\frac{\alpha^{2}}{6\,\pi}\,\omega^{4}\,.$ (27)

The momentum flux in the incident beam is given by
$\mathbf{E}\times\mathbf{B}$, where $\mathbf{E}$ and $\mathbf{B}$ are the
quantized electric and magnetic field operators. The net momentum carried by
the scattered radiation vanishes, so the the force in the $z$-direction on the
particle may be written as

$F^{z}=\sigma_{R}\,(\mathbf{E}\times\mathbf{B})^{z}=\frac{\alpha^{2}}{6\,\pi}\;R^{z}\,,$
(28)

where

$R^{z}=(\mathbf{\ddot{E}}\times\mathbf{\ddot{B}})^{z}$ (29)

may be viewed as the radiation pressure operator for polarizable particles.
Note that the factor of $\omega^{4}$ in the Rayleigh cross section has
resulted in two time derivatives on each of the electromagnetic fields. As a
consequence, $R^{z}$ is an operator with $p=7$.

We may average $F^{z}$ in time by letting $\alpha=\alpha(t)$ and integrating
along the worldline of the atom. A specific model for the origin of this time
dependence will be presented in the next section. We take the time dependence
to be proportional to the two-scale sampling function, $f_{2}(t)$, discussed
in Sec. II.4 with $t_{0}\gg\tau$. Here we set
$\alpha^{2}(t)\propto\alpha_{0}^{2}\,f_{2}(t)$, where $\alpha_{0}^{2}$ is the
time average of $\alpha^{2}(t)$ The averaged force now becomes
$\bar{F^{z}}=\alpha_{0}^{2}\,\bar{R^{z}}/(6\pi)$, where

$\bar{R^{z}}=\int f_{2}(t)\,R^{z}(t)\,dt\,,$ (30)

The probability distribution for large fluctuations of $\bar{R^{z}}$, and
hence of $\bar{F^{z}}$, is given by Eqs. (12), (13), and (14), with $p=7$ and
$\eta$ determined by $\alpha^{2}(t)$.

Let $x=\tau^{8}\,\bar{R^{z}}$ be the dimensionless measure of the quantum
radiation pressure fluctuations. For $1\ll x<x_{*}$, its probability
distribution has the asymptotic form in Eq. (12), with $c=\eta/7$ and
$b=-(8\eta+5)/7$ from Eqs. (13) and (14). The remaining constants in Eq. (12)
may be found from Eqs. (100) and (101) in Ref. FF2015  with
$\alpha\rightarrow\eta$, $B_{0}=4$, and $B=1/(6\pi^{2})$. (See Eq. (35) in
Ref. HF17 .) In addition we now set $\gamma=\gamma_{0}\,\tau/t_{0}$ and
$f(0)\rightarrow\tau\,f_{2}(0)\approx\tau/t_{0}$. The last factor of $\tau$
arises because Ref. FF2015  uses units where $\tau=1$, so $f(0)$ would be
$\tau\,f(0)$ in general units. Thus we find

$a=2\beta\,\left(\frac{\tau}{3\pi\,t_{0}}\right)^{-\eta/7}\,,$ (31)

and

$c_{0}=\frac{\gamma_{0}^{2}}{128\,\pi^{2}\,(\beta\,\eta)^{8}}\;\left(\frac{\tau}{3\pi\,t_{0}}\right)^{2(4\eta-7)/4}\,.$
(32)

If we set $\eta=1/2$ then the constants in Eq. (12) become

$c=\frac{1}{14}\,,\qquad b=-\frac{9}{7}\,,\qquad
a=2\beta\,\left(\frac{\tau}{3\pi\,t_{0}}\right)^{-\frac{1}{14}}\,,$ (33)

and

$c_{0}=\frac{2\,\gamma_{0}^{2}}{\pi^{2}\,\beta^{8}}\;\left(\frac{\tau}{3\pi\,t_{0}}\right)^{\frac{2}{7}}\,.$
(34)

Now the asymptotic probability density has the form

$P(x)=c_{0}\,x^{\frac{9}{7}}\,{\rm e}^{-a\,x^{\frac{1}{14}}}\,.$ (35)

If $p=7$ and $\eta=1/2$, the upper limit for the validity of the worldline
approximation, given by Eq. (15) and the value of $a$ in Eq. (33), becomes

$x_{*}=\beta^{-14}\,\left(\frac{\tau}{3\pi\,t_{0}}\right)\,\left(\frac{\tau}{\ell}\right)^{7}\,.$
(36)

It will be useful to have a measure of the actual probability of a fluctuation
of the order of magnitude of a given value of $x$. For this purpose, we define
the quantity

${\cal P}(x)=\int_{\frac{1}{2}x}^{x}P(y)\,dy\,,$ (37)

which is probability of finding an outcome in the interval $\frac{1}{2}x\leq
y\leq x$. Use the form of $P$ in Eq. (35) and define $u=y^{\frac{1}{14}}$ to
write

${\cal P}(x)=14c_{0}\int_{u_{1}}^{u_{2}}u^{-5}\,{\rm e}^{-a\,u}\,du\approx
14\,c_{0}\,{u_{2}}^{-5}\,{\rm e}^{-a\,u}\,\Delta u\,,$ (38)

where $u_{1}=(\frac{1}{2}x)^{\frac{1}{14}}$, $u_{2}=x^{\frac{1}{14}}$ , and
$\Delta u=u_{2}-u_{1}\approx 0.0483\,u_{2}$. This leads to

${\cal P}(x)\approx 0.676\,x\,P(x)=0.676\,c_{0}\,x^{-\frac{2}{7}}\,{\rm
e}^{-a\,x^{\frac{1}{14}}}\,.$ (39)

Recall that in Sec. II.3 we discussed how the probability distribution
enforces repayment of an energy loan from the vacuum. Similarly, a radiation
pressure fluctuation on an atom amounts to a temporary loan of linear
momentum. Suppose that we hold $\tau$ fixed but increase $t_{0}$. In this
case, $x=\tau^{8}\,\bar{R^{z}}$ does not change. However, if $\eta=1/2$, then
from Eqs. (31) and (34), we have $a\propto t_{0}^{\frac{1}{14}}$ and
$c_{0}\propto t_{0}^{-\frac{2}{7}}$. Both of these effects act to decrease
${\cal P}(x)$ as $t_{0}$ increases.

##  III Switched Polarizability and Velocity Fluctuations of Atoms

In this section, we present a model for the origin and effects of a time-
dependent polarizability of an atom.

###  III.1 Switching with Rydberg Atoms

Rydberg atoms are hydrogen-like atoms in highly excited states, and have been
extensively investigated in recent years. For a review, see for example,
Gallagher Gallagher . If $n$ is the principal quantum number of the excited
state, several physical properties of the Rydberg atom scale as powers of $n$.
These include the mean radius, $r_{0}\propto n^{2}$, the polarizability,
$\alpha\propto n^{7}$, and the radiative lifetime, $\tau_{rad}\propto n^{3}$.
See Table 2 in Ref. Gallagher . This means that the polarizability, can
increase many orders of magnitude from its value in the ground state when the
atom is excited, maintain this value for a time as long as $\tau_{rad}$, and
then return to a much smaller value when the atom is de-excited. The
excitation can be produced by a short laser pulse, which at optical
frequencies can have a duration of the order of $1{\rm fs}$. This pulse is
assumed to be described by a compactly supported function of time, and a
frequency spectrum of the form of Eq. (5). The de-excitation can be induced by
a second such pulse, which causes the polarizability to decrease rapidly
either by stimulated emission or by ionization. In this case, the
polarizability $\alpha(t)$ may have a time-dependence of the form illustrated
in Fig. 1.

We may take

$\alpha^{2}(t)=\alpha_{0}^{2}\,t_{0}\,f_{2}(t)\,,$ (40)

where $\alpha_{0}$ is the polarizability in the excited state and $f_{2}(t)$
is a two-scale sampling function of the type discussed in Sec. II.4. The
factor of $t_{0}$ arises because $f_{2}(t)$ satisfies Eq. (1). Recall that
$f_{2}(t)$, and hence $\alpha^{2}(t)$, increases and later decreases on a time
scale $\tau$, the duration of the exciting and de-exciting pulses. Here we
take $t_{0}\gg\tau$ to be the approximate temporal duration of the excited
state. Because $\alpha^{2}(t)$ increases and then decreases by many orders of
magnitude, we may take its initial and final values to be zero.

This time dependence leads to vacuum radiation pressure fluctuations, with the
averaged radiation pressure in the $z$-direction being given by $\bar{R^{z}}$
defined in Eq. (30). Note that here the $z$-direction is randomly selected and
need not be correlated with the direction of either laser pulse. If the
excited state is approximately spherically symmetric, then we expect vacuum
radiation pressure fluctuations to be equally likely in any direction. In this
case, $\langle\bar{R^{z}}\rangle=0$, and all other odd moments will also
vanish. However, the variance, $\langle(\bar{R^{z}})^{2}\rangle$ and other
even moments are nonzero. This leads to a symmetric probability probability
distribution, $P(x)=P(-x)$, as fluctuations of either sign and in any
direction are equally probable. The asymptotic form of $P(x)$ for large $x$
can be taken to be given by Eq. (12), where the parameters such as $\beta$ and
$\eta$ are determined by the Fourier transforms of the pulses which cause the
switching of the polarizability. In particular, if these pulses have
$\eta=1/2$, then $P(x)$ takes the form in Eq. (35) and falls very slowly.

###  III.2 Velocity Fluctuations: Numerical estimates

Now we wish to estimate the magnitude and probability of the velocity
fluctuations of the atom due to vacuum radiation pressure fluctuations. For
simplicity, we assume that the fluctuations may be describe in the worldline
approximation and with $\eta=1/2$, so we may use Eq. (35) with $x\lesssim
x_{*}$, with $x_{*}$ given by Eq. (36). We take the probability of a
fluctuation with $x\approx x_{*}$ to be ${\cal P}(x_{*})$, given by Eq. (39).
Such a fluctuation produces a radiation pressure of order $x_{*}/\tau^{8}$ and
a force which magnitude is of the order of

$\bar{F}\approx\frac{\alpha_{0}^{2}}{6\pi\,\tau^{8}}\,x_{*}$ (41)

and a duration of the order of $t_{0}$. The characteristic recoil speed will
be of order

$\bar{v}\approx\frac{\bar{F}\,t_{0}}{m}\approx\frac{\alpha_{0}^{2}\,t_{0}}{6\pi\,m\,\tau^{8}}\,x_{*}$
(42)

where $m$ is the mass of the atom.

We may use Eq. (36) and set $\ell\approx r_{0}$, the characteristic size of
the Rydberg atom, to write

$\bar{v}\approx\frac{\alpha_{0}^{2}}{18\,\pi^{2}\,m\,\beta^{14}\,r_{0}^{7}}\,.$
(43)

There are two remarkable feature of this expression: first, the dependence
upon both $\tau$ and $t_{0}$ have cancelled. Second, the result is very
sensitive to the parameter $\beta$, which appears in the asymptotic form of
the Fourier transform of the sampling function, Eq. (5), and hence in the
frequency spectrum of the exciting and de-exciting laser pulses. The meaning
of this parameter was discussed in Sec. II.5. However, we should note that
$\beta$ makes its appearance in $\bar{v}$ through the $\beta$-dependence of
$x_{*}$, and hence linked with our restriction to using the worldline
approximation. We may express Eq. (43) as

$\bar{v}\approx\frac{2.1\times 10^{-8}}{\beta^{14}}\,\left(\frac{1{\rm
u}}{m}\right)\,\left(\frac{\alpha_{0}}{n^{7}\,a_{0}^{3}}\right)^{2}\,\left(\frac{n^{2}\,a_{0}}{r_{0}}\right)^{7}\,,$
(44)

where $a_{0}$ is the Bohr radius of hydrogen. Note our estimate for $\bar{v}$
is independent of $n$, and that we expect $\alpha_{0}\approx n^{7}\,a_{0}^{3}$
and $r_{0}\approx n^{2}\,a_{0}$, in which case the last two factors in the
above expression will be of order one. Recall that we are working in units
with the speed of light set to one, so $\bar{v}=2.1\times 10^{-8}=6.3\,{\rm
m/s}$. Note the very strong dependence upon the value of $\beta$. For example,
if $\beta=1/\sqrt{2}$, and the other factors in Eq. (43) are of order one, we
would obtain the estimate $\bar{v}\approx 2.7\times 10^{-6}\approx 800\,{\rm
m/s}$.

We can compare these estimates with the expected atomic speeds due either to
thermal motion, or the recoil from photon absorption or emission. At absolute
temperature $T$, the root-mean-square thermal speed is

$v_{T}=\left(\frac{3\,k_{B}\,T}{m}\right)^{\frac{1}{2}}\approx 5\times
10^{-10}\,\left(\frac{T}{1\mu{\rm K}}\right)^{\frac{1}{2}}\;\left(\frac{1{\rm
u}}{m}\right)^{\frac{1}{2}}\,.$ (45)

Thus at sufficiently low temperatures, the thermal motion can be smaller than
the effects of vacuum radiation pressure fluctuations. If an atom at rest
either absorbs or emits a photon with energy $E_{\gamma}$, the magnitude of
the recoil momentum of the atom will equal the photon’s momentum, resulting in
a recoil velocity of magnitude

$v_{R}\approx 10^{-9}\,\left(\frac{E_{\gamma}}{1{\rm
eV}}\right)\,\left(\frac{1{\rm u}}{m}\right)\,.$ (46)

This can also be smaller than the effects of vacuum radiation pressure
fluctuations, and is also correlated with the photon’s direction.

Let us now turn to estimates of the probability ${\cal P}(x_{*})$ for the
fluctuations in question. Equations (33), (34), (36), and (39) lead to

${\cal P}(x_{*})\approx
0.14\,\gamma_{0}^{2}\,\beta^{-4}\,\left(\frac{r_{0}}{\tau}\right)\,{\rm
e}^{-2\,\sqrt{\tau/r_{0}}}\,.$ (47)

This probability also increases as $\beta$ decreases, although not so strongly
as does $\bar{v}$. Recall that $r_{0}\approx a_{0}\,n^{2}$. If we set
$\tau=1\,{\rm fs}\approx 3\times 10^{-5}\,{\rm cm}$, then

$\frac{\tau}{r_{0}}\approx\frac{6000}{n^{2}}\approx\left(\frac{77}{n}\right)^{2}\,.$
(48)

Recall that for the worldline approximation to have any range of validity, we
need $n\lesssim 77$. In Fig. 4, the value of the expression for ${\cal
P}(x_{*})$ given by Eq. (47) is plotted as a function of $n$. If both $\beta$
and $\gamma_{0}$ are of order one, then we can have ${\cal P}(x_{*})$ of the
order of $1\%$ while $\tau>r_{0}$.

Figure 4: The estimate of the probability of a radiation pressure fluctuation
in a given measurement is plotted as a function of the principal quantum
number $n$. Here we have set $\tau=1{\rm fs}$.

We have presented estimates both for the mean speed, $\bar{v}$, and the
probability, ${\cal P}(x_{*})$, of a vacuum radiation pressure fluctuation
capable of producing this mean speed. Both of these estimates assumed the
worldline approximation, in which the radiation pressure operator is averaged
in time only, and are hence restricted to considering fluctuations for which
$x\lesssim x_{*}$. This restriction is probably responsible for the strong
dependence of ${\cal P}(x_{*})$, and especially of $\bar{v}$, on the parameter
$\beta$. This may be seen from the sensitivity of $x_{*}$ to the value of
$\beta$ exhibited in Eq. (36). Thus it is reasonable to expect that including
spatial averaging, which will allow consideration of larger fluctuations, may
reduce the dependence of the results upon $\beta$. This will be a topic for
future work.

##  IV Summary

This paper has dealt with vacuum radiation pressure fluctuations on Rydberg
atoms. Such fluctuations are governed by a probability distribution which
decreases very slowly and is very sensitive to the details of how the
radiation pressure is measured. The relevant quantum operator is quadratic in
second time derivatives of the electromagnetic field. Like other stress
tensor-type operators for the quantized electromagnetic field, it must be
averaged in time before its fluctuations become finite. This averaging
describes the measurement process of the radiation pressure. If the
measurement is to begin and end at finite times, then the temporal averaging
function must have compact support, meaning that it is a non-analytic, but
$C^{\infty}$, function which vanishes outside of a finite time interval. Here
a specific model for the measurement has been proposed, in which an atom is
excited from its ground state to a highly excited state, forming a Rydberg
atom with greatly increased polarizability. The excitation is assumed to be
caused by a short laser pulse of finite duration, after which the excited
state persists for a time long compared to the excitation time. A second short
pulse ends the excitation by either ionization or stimulated emission.

The effect of the resulting time-dependent atomic polarizability is to measure
the vacuum radiation pressure on the atom, which in turn causes a temporary
recoil of the atom that may be large enough to be observed. The atom in effect
borrows some linear momentum from the quantum vacuum for a finite time. The
probability distribution for large radiation pressure fluctuations decreases
remarkably slowly at a rate determined by the pulse frequency spectrum, such
as an exponential of the $\frac{1}{14}$ power of the pressure fluctuation. It
is this slow rate of decrease which allows the possibility of observable
recoil with non-negligible probability.

The treatment in this paper has assumed a worldline approximation in which the
radiation pressure operator is averaged in time, but not in space. This
approximation simplifies the analysis, but limits the magnitude of the
fluctuations which may be considered. An expanded treatment including spatial
averaging, using either the analytic methods of Refs. FF2015 ; FF20 , or the
numerical methods of Refs. SFF18 ; WSF21  will be a topic for future research.

###### Acknowledgements.

I would like to thank Haiyun Huang, Enrico Schiappacasse, and Peter Wu for
useful conversations. This work was supported in part by the National Science
Foundation under Grant PHY-1912545.

## References

  * (1) C.J. Fewster, L.H. Ford and T.A. Roman, Probability distributions of smeared quantum stress tensors, Phys. Rev. D 81, 121901 (2010), arXiv:1004.0179 [quant-ph].
  * (2) C. J. Fewster, L. H. Ford and T. A. Roman, Probability distributions for quantum stress tensors in four dimensions, Phys. Rev. D 85, 125038 (2012), arXiv:1204.3570 [quant-ph].
  * (3) C. J. Fewster and L. H. Ford, Probability Distributions for Quantum Stress Tensors Measured in a Finite Time Interval, Phys. Rev. D 92, 105008 (2015), arXiv:1508.02359 [hep-th].
  * (4) E. D. Schiappacasse, C. J. Fewster and L. H. Ford, Vacuum Quantum Stress Tensor Fluctuations: A Diagonalization Approach, Phys. Rev. D 97, 025013 (2018), arXiv:1711.09477 [hep-th].
  * (5) C. J. Fewster and L. H. Ford, Probability Distributions for Space and Time Averaged Quantum Stress Tensors, Phys. Rev. D 101, 025006 (2020), arXiv:1909.07295
  * (6) P. Wu, E. D. Schiappacasse and L. H. Ford, Space and Time Averaged Quantum Stress Tensor Fluctuations, manuscript in preparation.
  * (7) L. H. Ford, S. P. Miao, Kin-Wang Ng, R. P. Woodard, and Chun-Hsien Wu, Quantum stress tensor fluctuations of a conformal field and inflationary cosmology, Phys. Rev. D 82, 043501 (2010).
  * (8) Chun-Hsien Wu, Jen-Tsung Hsiang, L. H. Ford, and Kin-Wang Ng, Gravity waves from quantum stress tensor fluctuations in inflation, Phys. Rev. D 84, 103515 (2011).
  * (9) P. Wu and L. H. Ford, Large Zero Point Density Fluctuations in Fluids, Phys. Rev. Research 2, 032028 (2020), arXiv:2005.04266.
  * (10) H. Huang and L. H. Ford, Vacuum Radiation Pressure Fluctuations and Barrier Penetration, Phys. Rev. D 96, 016003 (2017); arXiv:1610.01252
  * (11) R.F. Streater and A.S. Wightman, PCT, Spin, Statistics, and All That, (Princeton, 1989).
  * (12) T.A. Welton, Some Observable Effects of the Quantum Mechanical Fluctuations of the Electromagnetic Field, Phys. Rev., 74 1157, (1948).
  * (13) H. Huang and L. H. Ford, Quantum Electric Field Fluctuations and Potential Scattering, Phys. Rev. D 91, 125005 (2015); arXiv:1503.02962.
  * (14) V.V. Flambaum and V.G. Zelevinsky, Radiation Corrections Increase Tunneling Probability, Phys. Rev. Lett 83, 3108 (1999).
  * (15) S.G. Johnson, Saddle point integration of $C_{\infty}$ “bump” functions, arXiv: 1508.04376.
  * (16) T. F. Gallagher, Rydberg Atoms, Rep. Prog. Phys. 51, 143 (1988).