# Vacuum torque, propulsive forces, and anomalous tangential forces: Effects
of nonreciprocal media out of thermal equilibrium

Kimball A. Milton<EMAIL_ADDRESS>H. L. Dodge Department of Physics and
Astronomy, University of Oklahoma, Norman, OK 73019, USA  Xin Guo
<EMAIL_ADDRESS>H. L. Dodge Department of Physics and Astronomy, University
of Oklahoma, Norman, OK 73019, USA  Gerard Kennedy<EMAIL_ADDRESS>School of Mathematical Sciences, University of Southampton, Southampton, SO17
1BJ, UK  Nima Pourtolami<EMAIL_ADDRESS>National Bank of Canada,
Montreal, Quebec H3B 4S9, Canada  Dylan M. DelCol<EMAIL_ADDRESS>H.
L. Dodge Department of Physics and Astronomy, University of Oklahoma, Norman,
OK 73019, USA

###### Abstract

From the generalized fluctuation-dissipation theorem, it is known that a body
at rest made of nonreciprocal material may experience a torque, even in
vacuum, if it is not in thermal equilibrium with its environment. However, it
does not experience self-propulsion in such circumstances, except in higher
order. Nevertheless, such a body may experience both a normal torque and a
lateral force when adjacent to an ordinary surface with transverse
translational symmetry. We explore how these phenomena arise, discuss what
terminal velocities might be achieved, and point out some of the limitations
of applying our results to observations, including the Lorenz-Lorentz
correction, and the cooling due to radiation. In spite of these limitations,
the effects discussed would seem to be observable.

##  I Introduction

There is a long history of theoretical predictions of quantum or Casimir
friction, where a particle or extended body that moves parallel to a surface
experiences a force opposing its motion. The subject seems to have originated
with Teodorovich teod  and Levitov lev . For a selected bibliography on this
subject before 2016, see Ref. rev . The friction is typically conceived to
arise because of dissipation in the surface. For a subset of papers on this
subject, the reader is referred to Refs. pendry ; volokitin ; dedkov ; barton
; intravaia ; davies . For a readable overview, see Ref. marty . Quantum
friction can also result if the body itself is made of dissipative material.
However, much earlier, it was recognized that, even in vacuum, a moving body
or an atom without dissipation will experience friction due to the surrounding
radiation field—this is the famous Einstein-Hopf effect EinsteinHopf .

In previous papers, we have considered quantum vacuum friction due to field
and dipole fluctuations guo ; guo2 . For low velocities, the condition that
the atom or particle not gain or lose energy, the Nonequilibrium Steady State
condition (NESS) reiche , implies equal temperatures of the body and the
environment, while relativistic velocities typically imply that the body be
substantially hotter than the environment. (For other earlier work on
nonequilibrium friction, see Refs. dedkov2 ; pieplow ; vp  for example.) The
forces we considered there are true frictions, in that they always oppose the
motion, and they vanish at zero velocity and at zero temperature.

Here, we consider forces and torques that arise in the vacuum, or near other
bodies, when the relative velocity is zero. This requires not only that the
system be out of equilibrium, so that the temperature of the body is different
from that of its environment, but also that the electrical properties that
characterize the material constituting the bodies be exotic, “nonreciprocal,”
at least in lowest order. Nonreciprocity seems not to be possible for an
isolated body; the typical way it can be achieved is through the introduction
of an external magnetic field, or some other appropriate external influence.
Thus, it is something of an oxymoron to discuss nonreciprocal vacuum torque or
friction.

Of course, there is much earlier work on such nonequilibrium phenomena,
involving heat transfer, torque, and nonreciprocal surface forces krugeremig ;
Ott ; maghrebi ; Khandekar ; khandekar0 . (Further references will be provided
as our discussion continues.)

The outline of this paper is as follows. In Sec. II we discuss how the
fluctuation-dissipation theorem is modified for a nonreciprocal
susceptibility. We then display, in Sec. III, a simple model of such a
nonreciprocal material as a simplification of that given in Ref. guofan . The
corresponding quantum vacuum torque, first found in Refs. fan ; strekha , is
derived in Sec. IV. In Sec. V we rederive the modified torque if the body is
slowly rotating, which was first worked out in Ref. guofan . If the body is
hotter or not too much colder than the environment, the ordinary quantum
frictional torque acts as a drag, and the body acquires a terminal angular
velocity which should be readily observable, provided this temperature
difference can be maintained. In Sec. VI, the effect on the torque of an
underlying plate, be it a dielectric slab or a perfect conductor, is
investigated. Then we turn to the force, which can, of course, be inferred
from the torque. The quantum vacuum force is shown to vanish, in the weak-
susceptibility approximation that we use [Sec. VII], while, if an underlying
surface is present, there is a component of the force parallel to the surface,
as shown in examples of an imperfectly and a perfectly conducting plate [Sec.
VIII]. Again, if the nonequilibrium temperature difference could be
maintained, a substantial terminal velocity could be achieved. In Sec. IX we
calculate the time it would take for a body at rest to reach thermal
equilibrium with the environment, unless some mechanism were supplied to keep
it hotter or colder than the background. Possible suppression effects due to
the Lorenz-Lorentz (Clausius-Mossotti) correction are discussed in Sec. X,
although the resulting torques and forces should still be experimentally
measurable. Conclusions round out the paper.

Throughout, we use Heaviside-Lorentz (HL) electromagnetic units. We also set
$\hbar=c=1$, except when numerical values are given. It should be noted that
many authors (including some of us) often use Gaussian (G) units, for which
the polarizability differs by a factor of $4\pi$.

##  II Generalized fluctuation-dissipation theorem

Let $\mathbf{x}(t)$ be some dynamical variable (such as an electric dipole
moment). In terms of its frequency Fourier transform, the fluctuation-
dissipation theorem (FDT) tells us the expectation value of the symmetrized
quadratic product of the frequency-transformed variables:

$\langle
S\mathbf{x}(\omega)\mathbf{x}(\nu)\rangle=2\pi\delta(\omega+\nu)\Im\bm{\chi}(\omega)\coth\frac{\beta\omega}{2},$
(1)

where $\beta=1/T$ is the inverse temperature of the system, and $\bm{\chi}$ is
the generalized susceptibility; alternatively, the fluctuating quantity might
be the electric field, driven by the electric polarization, and the
susceptibility would be the retarded electric Green’s function.

Typically, we regard the susceptibility tensor as diagonal, or at least
symmetric. More generally, the “imaginary part” that occurs in the FDT is the
anti-Hermitian part,

$\Im\bm{\chi}=\frac{1}{2i}(\bm{\chi}-\bm{\chi}^{\dagger}),$ (2)

that is,

$(\Im\bm{\chi})_{ij}(\omega)=\frac{1}{2i}[\chi_{ij}(\omega)-\chi_{ji}^{*}(\omega)]=\frac{1}{2i}[\chi_{ij}(\omega)-\chi_{ji}(-\omega)],$
(3)

which uses the fact that $\chi_{ij}(\omega)$ is the Fourier transform of a
real response function. Unusual properties emerge from this if $\bm{\chi}$ is
not symmetric: This means that $\Im\bm{\chi}$ has both real and imaginary
parts in the conventional sense. The real part is

$\operatorname{Re}(\Im\bm{\chi})_{ij}(\omega)=\frac{1}{2}[\operatorname{Im}\chi_{ij}(\omega)+\operatorname{Im}\chi_{ji}(\omega)],$
(4a) which is symmetric in the indices but odd in $\omega$. These are the
components that give rise to the quantum friction force and torque. The
imaginary part of $\Im\bm{\chi}$ is
$\operatorname{Im}(\Im\bm{\chi})_{ij}(\omega)=-\frac{1}{2}[\operatorname{Re}\chi_{ij}(\omega)-\operatorname{Re}\chi_{ji}(\omega)],$
(4b) which is antisymmetric in the indices but even in $\omega$.

This property, as we shall see, leads to unusual phenomena for nonreciprocal
bodies, spontaneous quantum torque and quantum propulsion. The term
“nonreciprocal” seems to have a variety of meanings in the literature; in this
paper we will take it to mean that Eq. (4b) is nonzero.

For an ordinary material, $\chi_{ij}(\omega)$ is symmetric in the indices,
which means that the anti-Hermitian part coincides with the usual imaginary
part:

$\chi_{ij}(\omega)=\chi_{ji}(\omega)\Rightarrow(\Im\bm{\chi})_{ij}(\omega)=\operatorname{Im}\chi_{ij}(\omega).$
(5)

Where a susceptibility depends on continuous coordinates as well, such as the
Green’s dyadic that describes the electric field, reciprocity means invariance
under interchange of discrete indices and continuous coordinates:

$\Gamma_{ij}(\mathbf{r,r^{\prime}};\omega)=\Gamma_{ji}(\mathbf{r^{\prime},r};\omega)\Rightarrow(\Im\bm{\Gamma})_{ij}(\mathbf{r,r^{\prime}};\omega)=\operatorname{Im}\Gamma_{ij}(\mathbf{r,r^{\prime}};\omega).$
(6)

It is easy to check that this is satisfied by the Green’s dyadic for a
dielectric half-space, for example, which has off-diagonal elements, symmetric
in the tensor indices. When the latter is expressed as a two-dimensional
Fourier transform, as is convenient when the environment consists of a
dielectric slab perpendicular to the $z$ axis,

$\Gamma_{ij}(\mathbf{r,r^{\prime}};\omega)=\int\frac{(d\mathbf{k}_{\perp})}{(2\pi)^{2}}e^{i\mathbf{k_{\perp}\cdot(r-r^{\prime})}_{\perp}}g_{ij}(z,z^{\prime};\mathbf{k}_{\perp},\omega),$
(7)

the FDT is expressed in terms of

$(\Im\mathbf{g})_{ij}(z,z^{\prime};\mathbf{k}_{\perp},\omega)=\frac{1}{2i}\left[g_{ij}(z,z^{\prime};\mathbf{k}_{\perp},\omega)-g_{ji}(z^{\prime},z;-\mathbf{k}_{\perp},-\omega)\right].$
(8)

##  III Model for nonreciprocal material

In order to create a nonreciprocal response, one needs an external influence.
Such is supplied by a magnetic field. Let us suppose a oscillator with damping
$\eta$ is driven by both an electric and a magnetic field:

$m\frac{d^{2}\mathbf{r}}{dt^{2}}+m\eta\frac{d\mathbf{r}}{dt}+m\omega_{0}^{2}\mathbf{r}=e\left(\mathbf{E}+\frac{d\mathbf{r}}{dt}\times\mathbf{B}\right).$
(9)

If the magnetic field lies in the $z$ direction, this immediately yields an
electric susceptibility that is nonsymmetric and nonreciprocal:

$\bm{\chi}=\omega_{p}^{2}\left(\begin{array}[]{ccc}\frac{\omega_{0}^{2}-\omega^{2}-i\omega\eta}{(\omega_{0}^{2}-\omega^{2}-i\omega\eta)^{2}-\omega^{2}\omega_{c}^{2}}&\frac{-i\omega\omega_{c}}{(\omega_{0}^{2}-\omega^{2}-i\omega\eta)^{2}-\omega^{2}\omega_{c}^{2}}&0\\\
\frac{i\omega\omega_{c}}{(\omega_{0}^{2}-\omega^{2}-i\omega\eta)^{2}-\omega^{2}\omega_{c}^{2}}&\frac{\omega_{0}^{2}-\omega^{2}-i\omega\eta}{(\omega_{0}^{2}-\omega^{2}-i\omega\eta)^{2}-\omega^{2}\omega_{c}^{2}}&0\\\
0&0&\frac{1}{\omega_{0}^{2}-\omega^{2}-i\omega\eta}\end{array}\right),$ (10)

in terms of the plasma frequency $\omega_{p}^{2}=ne^{2}/m$, $n$ being the
density of charges, and the cyclotron frequency $\omega_{c}=eB/m$. For a
metal, we would set the restoring force to zero, so $\omega_{0}=0$, and we
exactly recover the form given by Guo and Fan guofan . In particular, this
provides us with a model for the anti-Hermitian part of the susceptibility,

$\chi_{xy}=-\chi_{yx}=-i\frac{\omega_{p}^{2}\omega_{c}/\omega}{(\omega+i\eta)^{2}-\omega_{c}^{2}}.$
(11)

Numerically, for the charge and mass of the electron,

$\omega_{c}=\frac{eB}{m}=m\frac{B}{B_{c}},\quad
B_{c}=\frac{m^{2}}{e}=4.41\times 10^{9}\,\mbox{T},$ (12)

so for a magnetic field of strength 1 T, $\omega_{c}\sim 10^{-4}$ eV, far
smaller than the damping parameter for gold, for example, $\eta\approx
3.5\times 10^{-2}$ eV. Thus, to a good approximation, we can use for a metal,

$\chi_{xy}-\chi_{yx}\approx-2i\frac{\omega_{c}\omega_{p}^{2}}{\omega}\frac{1}{(\omega+i\eta)^{2}}.$
(13)

##  IV Quantum vacuum torque

The vacuum torque on a stationary body was discussed in general terms in Ref.
fan , and subsequently in Ref. strekha . Somewhat earlier, Ref. maghrebi
showed that a topologically insulating film in a magnetic field, out of
thermal equilibrium, experiences a torque, which seems to be an instance of
this same phenomenon. A torque on a body at rest requires that it be composed
of nonreciprocal material, which is characterized by having a real part of the
susceptibility which is nonsymmetric, and that the temperature of the body,
$T^{\prime}$, be different from that of the environment, $T$.

The torque on an arbitrary body, described by electric susceptibility,
$\bm{\chi}(\mathbf{r};\omega)$, is, in terms of the polarization,
$\mathbf{P}$,

$\bm{\tau}(t)=\int(d\mathbf{r})\,\mathbf{r}\times\left[-\bm{\nabla}\cdot\mathbf{P}(\mathbf{r},t)\mathbf{E}(\mathbf{r},t)+\frac{\partial\mathbf{P}(\mathbf{r},t)}{\partial
t}\times\mathbf{B}(\mathbf{r},t)\right].$ (14)

Writing this in terms of Fourier transforms, and eliminating $\mathbf{B}$ in
favor of $\mathbf{E}$, we have

$\displaystyle\bm{\tau}$ $\displaystyle=$
$\displaystyle\int(d\mathbf{r})\,\mathbf{r}\times\int\frac{d\omega}{2\pi}\frac{d\nu}{2\pi}e^{-i(\omega+\nu)t}\left\\{-\bm{\nabla}\cdot\mathbf{P(r};\omega)\mathbf{E(r};\nu)-\frac{\omega}{\nu}\mathbf{P(r};\omega)\times[\bm{\nabla}\times\mathbf{E(r};\nu)]\right\\}$
$\displaystyle=$
$\displaystyle\int(d\mathbf{r})\frac{d\omega}{2\pi}\frac{d\nu}{2\pi}e^{-i(\omega+\nu)t}\mathbf{r}\times\bigg{[}\frac{\omega}{\nu}\left\\{\bm{\nabla}\cdot[\mathbf{P(r;\omega)E(r;\nu)}]-\mathbf{P(r;\omega)}\cdot(\bm{\nabla})\cdot\mathbf{E(r;\nu)}\right\\}-\left(1+\frac{\omega}{\nu}\right)[\bm{\nabla}\cdot\mathbf{P(r;\omega)]E(r;\nu)}\bigg{]},$

where the notation in the first term in the last equality means that the
$\bm{\nabla}$ is dotted only with $\mathbf{P}$, but acts on both variables,
while in the second term $\bm{\nabla}$ is the vector crossed with
$\mathbf{r}$—the parentheses are intended to insulate it from the two vectors
surrounding it. Here, the source of the electric field is the electric
polarization,

$\mathbf{E}(\mathbf{r};\omega)=\int(d\mathbf{r^{\prime}})\,\bm{\Gamma}(\mathbf{r,r^{\prime}};\omega)\cdot\mathbf{P}(\mathbf{r}^{\prime};\omega),$
(16a) where $\bm{\Gamma}$ is the retarded electromagnetic Green’s dyadic,
while the polarization is linearly related to the electric field,
$\mathbf{P}(\mathbf{r};\omega)=\int(d\mathbf{r})\,\delta(\mathbf{r-r^{\prime}})\bm{\chi}(\mathbf{r};\omega)\cdot\mathbf{E(r^{\prime};\omega)}=\bm{\chi}(\mathbf{r};\omega)\cdot\mathbf{E}(\mathbf{r};\omega),$
(16b)

where we assume that the electric susceptibility is local in space. This means
that there are two origins for the quantum torque: field fluctuations and
dipole fluctuations. We evaluate the two contributions to the torque by use of
the FDT:

$\displaystyle\langle
SE_{i}(\mathbf{r};\omega)E_{j}(\mathbf{r}^{\prime};\nu)\rangle$
$\displaystyle=$ $\displaystyle
2\pi\delta(\omega+\nu)(\Im\bm{\Gamma})_{ij}(\mathbf{r,r^{\prime}};\omega)\coth\frac{\beta\omega}{2},\quad\beta=\frac{1}{T},$
(17a) $\displaystyle\langle
SP_{i}(\mathbf{r};\omega)P_{j}(\mathbf{r}^{\prime};\nu)\rangle$
$\displaystyle=$ $\displaystyle
2\pi\delta(\omega+\nu)\delta(\mathbf{r-r^{\prime}})(\Im\bm{\chi})_{ij}(\mathbf{r};\omega)\coth\frac{\beta^{\prime}\omega}{2},\quad\beta^{\prime}=\frac{1}{T^{\prime}},$
(17b)

where $S$ indicates that the symmetrized expectation values are used.
Therefore, the last term in Eq. (LABEL:torque2) vanishes, because the sum of
the two frequencies is zero, leaving us with111This is made up of the
“internal” and “external” torques, as given in Ref. ce , Eqs. (4.47) and
(4.46).

$\bm{\tau}=\int(d\mathbf{r})\frac{d\omega}{2\pi}\frac{d\nu}{2\pi}e^{-i(\omega+\nu)t}\left[\mathbf{P(r;\omega)\times
E(r;\nu)+P(r;\omega)\cdot(r\times\bm{\nabla})\cdot E(r;\nu)}\right].$ (18)

Here, the notation in the last term signifies that the free vector index is on
the angular momentum operator; the $\mathbf{P,E}$ are dotted together.

Using Eqs. (16b) and (17a), we find for the EE contribution to the torque

$\bm{\tau}_{i}^{\rm
EE}=\int(d\mathbf{r})\frac{d\omega}{2\pi}\coth\frac{\beta\omega}{2}\epsilon_{ijk}\left[\chi_{jl}(\mathbf{r};\omega)(\Im\bm{\Gamma})_{lk}(\mathbf{r,r};\omega)+\chi_{lm}(\mathbf{r};\omega)x_{j}\nabla^{\prime}_{k}(\Im\bm{\Gamma})_{ml}(\mathbf{r,r^{\prime}};\omega)\bigg{|}_{\mathbf{r-r^{\prime}\equiv
R}\to\bm{0}}\right].$ (19)

Here, $\bm{\Gamma}$ is taken to be the usual vacuum retarded Green’s dyadic,
$\bm{\Gamma}^{0}$, the divergenceless part of which can be written as222The
omitted term is proportional to $\delta(\mathbf{r-r^{\prime}})$, which does
not contribute in the point-split limit.

$\bm{\Gamma}^{0\prime}(\mathbf{r,r^{\prime}};\omega)=(\bm{\nabla\nabla-1}\nabla^{2})\frac{e^{i\omega
R}}{4\pi R}=\left[\bm{\hat{R}\hat{R}}(3-3i\omega
R-\omega^{2}R^{2})-\bm{1}(1-i\omega R-\omega^{2}R^{2})\right]\frac{e^{i\omega
R}}{4\pi R^{3}},$ (20)

where $R=|\mathbf{r-r^{\prime}}|$ and $\hat{\mathbf{R}}=\mathbf{R}/R$. It is
evident that the second term in Eq. (19) vanishes here because
$\operatorname{Im}e^{i\omega R}/R$ is a function of $R$. When (20) is
rotationally averaged in the coincidence limit ($\mathbf{R}\to\bm{0}$), we
obtain

$\bm{\Gamma}^{0\prime}(\mathbf{r,r^{\prime}};\omega)\to\bm{1}\left(\frac{\omega^{2}}{6\pi
R}+i\frac{\omega^{3}}{6\pi}+O(R)\right).$ (21)

Therefore, we are left with only a single term for the torque:

$\tau_{i}^{\rm
EE}=\int\frac{d\omega}{2\pi}\coth\frac{\beta\omega}{2}\epsilon_{ijk}\operatorname{Re}\alpha_{jk}(\omega)\frac{\omega^{3}}{6\pi},$
(22)

where the mean polarizability333 The nonlinear effects occurring in the
Lorenz-Lorentz law relate the polarizability to the permittivity; the
polarizability is implied thereby to be linear—see Ref. entropy , and Sec. X.
of the body444Note that there is no requirement that the body be spherical, so
any rotation would be observable. is given by

$\alpha_{jk}(\omega)=\int(d\mathbf{r})\chi_{jk}(\mathbf{r};\omega),$ (23)

the real part of which is picked out by the necessity of the integrand being
even in $\omega$. Thus, nonreciprocity is necessary for a vacuum torque in
first order.

For the PP fluctuation part, Eq. (18) is written as

$\bm{\tau}^{\rm
PP}=\int(d\mathbf{r})(d\mathbf{r^{\prime}})\frac{d\omega}{2\pi}\frac{d\nu}{2\pi}e^{-i(\omega+\nu)t}\left[\mathbf{P(r};\omega)\times\bm{\Gamma}(\mathbf{r,r^{\prime}};\nu)\cdot\mathbf{P(r}^{\prime};\nu)+\mathbf{P(r};\omega)\cdot(\mathbf{r}\times\bm{\nabla})\cdot\bm{\Gamma}(\mathbf{r,r^{\prime}};\nu)\cdot\mathbf{P(r^{\prime}};\nu)\right],$
(24)

to which the FDT (17b) is to be applied. Again, for vacuum, the coincidence
limit of the second term is zero, and because of the diagonal form of the
limit of the Green’s dyadic, only the antisymmetric part of the susceptibility
survives upon use of Eq. (17b). That is, the $i$th component of the quantity
in square brackets in Eq. (24) becomes

$\delta(\mathbf{r-r^{\prime}})2\pi\delta(\omega+\nu)\epsilon_{ijk}(\Im\bm{\chi})_{jk}(\mathbf{r};\omega)\left(-i\frac{\omega^{3}}{6\pi}\right)\coth\frac{\beta^{\prime}\omega}{2}.$
(25)

Here we have recognized that the antisymmetric part of $\Im\bm{\chi}$ is even
in $\omega$, according to Eq. (4b), and therefore only the odd part of the
vacuum Green’s dyadic survives. The resulting torque is thus of the same form
as for the EE contribution, except for the sign, and the replacement
$\beta\to\beta^{\prime}$.

The combination of the two contributions thus yields the torque on a
nonreciprocal body in vacuum, when the temperature of the body, $T^{\prime}$,
differs from that of the blackbody radiation, $T$, due to PP and EE
fluctuations:

$\tau_{i}=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\frac{\omega^{3}}{6\pi}\epsilon_{ijk}\operatorname{Re}\alpha_{jk}(\omega)\left[\coth\frac{\beta\omega}{2}-\coth\frac{\beta^{\prime}\omega}{2}\right].$
(26)

This result exactly agrees with that of Guo and Fan fan , for zero rotational
velocity, and with that of Strekha et al. strekha . However, there is no
quantum vacuum force in this static situation, at least in first order, which
we will demonstrate in Sec. VII. (Actually, this is already evident from the
vanishing of the external torque contribution.)

Let us use the model in Eq. (13) to give an estimate of the size of the
torque. Inserting this into Eq. (26) and letting $\omega=\eta x$ gives

$\tau_{z}=\frac{\eta\omega_{c}\omega_{p}^{2}V}{3\pi^{2}}\int_{-\infty}^{\infty}dx\frac{x^{3}}{(x^{2}+1)^{2}}\left(\coth\frac{\beta^{\prime}\eta
x}{2}-\coth\frac{\beta\eta
x}{2}\right)=\frac{4\eta\omega_{c}\omega_{p}^{2}V}{3\pi^{2}}[I_{2}(\beta^{\prime}\eta)-I_{2}(\beta\eta)],$
(27)

where $V$ is the volume of the body, and the integrals are defined by Eq. (93)
of Appendix A. As expected, this is positive if $T^{\prime}>T$. These
integrals are readily evaluated in Eq. (94):

$\tau_{z}=\frac{\eta\omega_{c}\omega_{p}^{2}V}{3\pi^{2}}\left[\frac{\pi}{\eta}(T-T^{\prime})+2\ln\frac{T^{\prime}}{T}+2\psi\left(\frac{\eta}{2\pi
T}\right)-2\psi\left(\frac{\eta}{2\pi T^{\prime}}\right)+\frac{\eta}{2\pi
T}\psi^{\prime}\left(\frac{\eta}{2\pi T}\right)-\frac{\eta}{2\pi
T^{\prime}}\psi^{\prime}\left(\frac{\eta}{2\pi T^{\prime}}\right)\right].$
(28)

The torque, and the approximations, are shown in Fig. 1.

Figure 1:  The quantum vacuum torque, apart from the prefactor, in Eqs. (27)
or (28) is shown as a function of $T^{\prime}/\eta$ by the solid line for
$T/\eta=0.714$, appropriate for an environment at room temperature, 300 K, and
a gold body. The dashed blue line is the high-temperature approximation (96b),
while the low-temperature approximation (98b) is shown by the dashed red line.
In both limiting cases, the enviromental temperature is treated exactly.

For a gold ($\omega_{p}=9$ eV, $\eta=0.035$ eV) nanosphere of radius 100 nm,
with $\omega_{c}=10^{-4}$ eV, the prefactor in Eq. (28) is $8\times 10^{-25}$
Nm.

##  V Torque on a rotating body

Of course, a torque on a body will cause it to rotate. So, what is the torque
on a rotating body? Naturally, there should be a vacuum torque on a rotating
body made of ordinary (reciprocal) material, just as there is quantum vacuum
friction on a linearly moving body guo ; guo2 . The nonreciprocal aspect of
this torque was first treated in Ref. fan .

We consider a body rotating about the $z$ axis passing through its center of
mass with angular velocity $\Omega$. The formula (18) should still apply, with
the external torque (the second term there) still not contributing if the
background is vacuum. However, the polarization and electric fields should now
refer to the body (rotating) frame, denoted by a prime subsequently. For low
velocities, these are related to those in the blackbody (unprimed) frame by a
rotation:

$\displaystyle E_{x}^{\prime}(\mathbf{r^{\prime}},t)$ $\displaystyle=$
$\displaystyle E_{x}(\mathbf{r},t)\cos\Omega t+E_{y}(\mathbf{r},t)\sin\Omega
t,$ (29a) $\displaystyle E_{y}^{\prime}(\mathbf{r^{\prime}},t)$
$\displaystyle=$ $\displaystyle E_{y}(\mathbf{r},t)\cos\Omega
t-E_{x}(\mathbf{r},t)\sin\Omega t,$ (29b)

which means for the frequency transforms,

$\displaystyle E^{\prime}_{x}(\mathbf{r}^{\prime};\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[E_{x}(\mathbf{r},\omega_{+})+E_{x}(\mathbf{r};\omega_{-})\right]+\frac{1}{2i}\left[E_{y}(\mathbf{r};\omega_{+})-E_{y}(\mathbf{r};\omega_{-})\right],$
(30a) $\displaystyle E^{\prime}_{y}(\mathbf{r}^{\prime};\omega)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left[E_{y}(\mathbf{r},\omega_{+})+E_{y}(\mathbf{r};\omega_{-})\right]-\frac{1}{2i}\left[E_{x}(\mathbf{r};\omega_{+})-E_{x}(\mathbf{r};\omega_{-})\right],$
(30b)

where $\omega_{\pm}=\omega\pm\Omega$. $\mathbf{P}$ transforms in the same way.

The strategy followed to calculate the quantum rotational friction is the same
as that used for quantum rectilinear friction guo2 . There are two
contributions: field fluctuations and dipole fluctuations. For the former we
use Eq. (16b) to replace the polarization by the electric field, but now
understood in the body frame. Then we have to transform both electric fields
to the blackbody frame. For vacuum friction we can use Eq. (21) for the
Green’s dyadic that appears when the FDT is employed for the fields. We also
only keep half the terms: only those proportional to
$\delta(\omega_{\pm}+\nu_{\mp})$ do not average to zero in time. The result of
a straightforward calculation is

$\tau_{z}^{\rm
EE}=\int\frac{d\omega}{2\pi}\frac{\omega^{3}}{6\pi}\left[\operatorname{Im}(\alpha_{xx}+\alpha_{yy})(\omega_{-})+\operatorname{Re}(\alpha_{xy}-\alpha_{yx})(\omega_{-})\right]\coth\frac{\beta\omega}{2},$
(31)

where we have also noticed that under $\omega\to-\omega$,
$\omega_{+}\to-\omega_{-}$.

The procedure for the PP fluctuations is similar, except now we replace
$\mathbf{E}$ by $\mathbf{P}$ according to Eq. (16a). This holds in the
blackbody frame, so $\mathbf{P}$ must be transformed back to the body
(rotating) frame. Simplifications occur as before for the vacuum case, and we
find after a bit of algebra

$\tau_{z}^{\rm
PP}=-\int\frac{d\omega}{2\pi}\frac{\omega^{3}}{6\pi}\left[\operatorname{Im}(\alpha_{xx}+\alpha_{yy})(\omega_{-})+\operatorname{Re}(\alpha_{xy}-\alpha_{yx})(\omega_{-})\right]\coth\frac{\beta^{\prime}\omega_{-}}{2}.$
(32)

Thus, when the two contributions are added, we find for the torque on a
(slowly) rotating body:

$\tau_{z}=\frac{1}{12\pi^{2}}\int_{-\infty}^{\infty}d\omega\,\omega_{+}^{3}\left[\operatorname{Im}(\alpha_{xx}+\alpha_{yy})(\omega)+\operatorname{Re}(\alpha_{xy}-\alpha_{yx})(\omega)\right]\left(\coth\frac{\beta\omega_{+}}{2}-\coth\frac{\beta^{\prime}\omega}{2}\right).$
(33)

This is precisely the torque found in Ref. fan  (recall
$4\pi\alpha^{G}=\alpha^{\rm HL}$). This result for an isotropic (reciprocal)
particle was given in Ref. pan , which further considered the effect of a
magnetic field.555Earlier related works on forces and torques on bodiess with
various kinds of asymmetries include Refs. Reid ; Ott ; maghrebi ; Khandekar .
Note that if $\Omega=0$ the first term involving the diagonal polarizabilities
vanishes because the integrand is odd, and the second term reproduces Eq.
(26).

It is illuminating to expand this expression to leading order in the
rotational velocity $\Omega$ (this is the adiabatic approximation):

$\displaystyle\tau_{z}$ $\displaystyle=$
$\displaystyle\frac{1}{12\pi^{2}}\int_{-\infty}^{\infty}d\omega\,\omega^{3}\bigg{\\{}\operatorname{Re}\left(\alpha_{xy}-\alpha_{yx}\right)(\omega)\left[\coth\frac{\beta\omega}{2}-\coth\frac{\beta^{\prime}\omega}{2}\right]$
(34)
$\displaystyle\quad\mbox{}-\Omega\frac{3}{\omega}\operatorname{Im}\left(\alpha_{xx}+\alpha_{yy}\right)(\omega)\left[\coth\frac{\beta^{\prime}\omega}{2}-\coth\frac{\beta\omega}{2}\right]-\Omega\frac{\beta}{2}\operatorname{Im}\left(\alpha_{xx}+\alpha_{yy}\right)(\omega)\operatorname{csch}^{2}\frac{\beta\omega}{2}\bigg{\\}}.$

The first term here is the (nonreciprocal) quantum vacuum torque (26), the
second is the nonequilibrium contribution to the ordinary (reciprocal) quantum
vacuum frictional torque, and the third term is the analog of the Einstein-
Hopf quantum vacuum friction. The sum of the two frictional terms is a drag if
$T^{\prime}>T$. If $T^{\prime}<T$, the angular velocity changes sign, so
initially the friction remains a drag, but for sufficiently low temperatures,
$T^{\prime}$, the second term in Eq. (34) will dominate, and exponential
growth of the angular velocity will ensue, insofar as the low-velocity
approximation remains valid.

Note, that if the last two terms constitute a drag, the nonreciprocal torque
found here will lead to the body rotating with a constant terminal angular
velocity. Writing Eq. (34) in the abbreviated form,666That is,
$\tau_{0}=\tau_{z}(\Omega=0)$ and
$\tau_{1}^{\prime}=-\frac{d\tau_{z}}{d\Omega}(\Omega=0)$.

$\tau_{z}=I\dot{\Omega}=\tau_{0}-\Omega\tau^{\prime}_{1},$ (35)

where $I$ is the moment of inertia of the body, we immediately obtain

$\Omega(t)=\frac{\tau_{0}}{\tau^{\prime}_{1}}\left(1-e^{-\tau^{\prime}_{1}t/I}\right),$
(36)

if the body is not rotating at time $t=0$. The terminal velocity is
$\Omega_{T}=\Omega(t\to\infty)=\tau_{0}/\tau^{\prime}_{1}$, which might be
expected to be small. (This, of course, assumes that the particle and
environmental temperatures do not change. We will address the tendency toward
thermal equilibrium in Sec. IX.)

To proceed, let us again use the model (10), with $\omega_{0}=0$ and
$\omega_{c}\ll\eta$. Then we have

$\operatorname{Im}(\alpha_{xx}+\alpha_{yy})(\omega)=\frac{2\omega_{p}^{2}\eta}{\omega(\omega^{2}+\eta^{2})}V.$
(37)

The result is ($x=\omega/\eta$)

$\displaystyle\tau_{1}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{2\omega_{p}^{2}\eta
V}{3\pi^{2}}\left(3+\beta\frac{\partial}{\partial\beta}\right)\int_{0}^{\infty}dx\frac{x}{(x^{2}+1)}\left[\frac{1}{e^{\beta^{\prime}\eta
x}-1}-\frac{1}{e^{\beta\eta x}-1}\right]$ (38) $\displaystyle=$
$\displaystyle\frac{\omega_{p}^{2}\eta
V}{3\pi^{2}}\left\\{\frac{\pi}{\eta}(2T-3T^{\prime})+3\ln\frac{T}{T^{\prime}}-1+3\psi\left(\frac{\eta}{2\pi
T}\right)-3\psi\left(\frac{\eta}{2\pi T^{\prime}}\right)+\frac{\eta}{2\pi
T}\psi^{\prime}\left(\frac{\eta}{2\pi T}\right)\right\\}.$

Let us write, from Eqs. (28) and (38),

$\tau_{0}=\frac{\eta\omega_{c}\omega_{p}^{2}V}{\pi^{2}}f(T,T^{\prime}),\quad\tau^{\prime}_{1}=\frac{\eta\omega_{p}^{2}V}{\pi^{2}}g(T,T^{\prime}).$
(39)

Then, the terminal angular velocity is

$\Omega_{T}=\frac{\tau_{0}}{\tau^{\prime}_{1}}={\omega_{c}}\frac{f(T,T^{\prime})}{g(T,T^{\prime})}\sim\omega_{c}\sim
10^{-4}\mbox{eV}\sim 10^{11}\,\mbox{s}^{-1},$ (40)

perhaps surprisingly high, but very small compared to atomic frequencies. (The
terminal circumferential speed in this case is $\Omega_{T}R\sim 10^{4}$ m/s,
for a gold nanosphere of radius $R=100$ nm.) The relaxation time required to
reach this velocity is

$t_{0}=\frac{I}{\tau_{1}^{\prime}}\sim\frac{MR^{2}}{\omega_{p}^{2}\eta V}\sim
10^{6}\,\mbox{s},$ (41)

for the same parameters. The temperature dependence of the terminal angular
velocity, $\Omega_{T}$, is shown in Fig. 2.

Figure 2:  The terminal angular velocity of the nonreciprocal body when it is
hotter or colder, than its environment, which is taken to be at room
temperature. The limit of $\Omega_{T}$ for high particle temperature is
universally $\omega_{c}/3$, independent of the background temperature. When
the temperature is lower than that of the background, the frictional torque
initially acts as a drag, but for sufficiently low temperature, the
“frictional” terms change sign, and the angular velocity increases
exponentially without bound.

Note that if the temperature of the body is lower than that of the
environment, the terminal angular velocity is negative. Whether the body is
hotter or not too much colder than the environment, it will reach a terminal
velocity if the temperature difference is maintained. If the temperature of
the body is much colder than that of the environment (for the example shown,
about half room temperature), the frictional torques reverse sign, and no
bound to the angular velocity can be reached. A negative $\tau^{\prime}_{1}$
means exponential growth. Of course, before the angular velocity gets too
large, the nonrelativistic approximation used here breaks down, even if the
temperature difference can be maintained by some external or internal agent.

However, it is not necessary to wait a long time to reach the terminal
velocity, because the initial angular acceleration,

$\dot{\Omega}(0)=\frac{\Omega_{T}}{t_{0}}\sim 10^{5}\,\mbox{s}^{-2},$ (42)

should be easily discernible.

##  VI Torque in presence of a dielectric plate

What happens if the background is less trivial, say consisting of an isotropic
dielectric plate filling the halfspace $z<0$, while the body lies a distance
$a$ above it? Then, of course, there will be a torque on a nonspherical body
as well as a force, due to ordinary Casimir forces, even when the body is made
of ordinary reciprocal material. What would be unusual is if there were a
torque around the $z$ axis, since the environment possesses rotational
symmetry about that direction. For simplicity, we will assume that the entire
background, vacuum plus dielectric plate, is in equilibrium at temperature
$T$, while the nanoparticle has temperature $T^{\prime}$. In that case, the
$z$-component of the torque coming from field fluctuations is given by Eq.
(19). Now both terms can contribute, but only for a nonreciprocal body, which
might be irregularly shaped. For such an body, where

$\hat{\chi}_{ij}(\mathbf{r};\omega)=\operatorname{Re}(\chi_{ij}-\chi_{ji})(\mathbf{r};\omega)$
(43)

is nonzero, the normal torque component can be computed from the explicit
construction of the Green’s dyadic (for example, as given in Ref. guo ). Using
the Fourier representation (7), the integration over $k_{x}$ and $k_{y}$ will
vanish except for the $g_{xx}$ and $g_{yy}$ terms, for the first term in Eq.
(19), while $g_{xz,zx}$ and $g_{yz,zy}$ contribute for the second term,
yielding the scattering part of the torque:

$\displaystyle\tau_{z}^{s}$ $\displaystyle=$
$\displaystyle\int(d\mathbf{r})\frac{d\omega}{2\pi}\left[\coth\frac{\beta\omega}{2}-\coth\frac{\beta^{\prime}\omega}{2}\right]\frac{1}{4}\int\frac{(d\mathbf{k}_{\perp})}{(2\pi)^{2}}\bigg{\\{}\hat{\chi}_{xy}(\mathbf{r};\omega)\operatorname{Im}\left[\left(\kappa
r^{H}+\frac{\omega^{2}}{\kappa}r^{E}\right)e^{-2\kappa z}\right]$ (44)
$\displaystyle\quad\mbox{}+[\hat{\chi}_{yz}(\mathbf{r};\omega)x-\hat{\chi}_{xz}(\mathbf{r};\omega)y]k^{2}\operatorname{Im}\left[r^{H}e^{-2\kappa
z}\right]\bigg{\\}},$

where $\kappa=\sqrt{k^{2}-\omega^{2}}$, and the transverse magnetic and
electric reflection coefficients are

$r^{H}=\frac{\kappa-\kappa^{\prime}/\varepsilon(\omega)}{\kappa+\kappa^{\prime}/\varepsilon(\omega)},\quad
r^{E}=\frac{\kappa-\kappa^{\prime}}{\kappa+\kappa^{\prime}},\quad\kappa^{\prime}=\sqrt{k^{2}-\omega^{2}\varepsilon(\omega)}.$
(45)

Now the torque depends on the distribution of the anisotropic material across
the body, so is not describable by simply an effective nonreciprocal
polarizability.

Note further that the second two terms in Eq. (44), proportional to $x$ and
$y$, respectively, depend on the position of the body as well as the
distribution of material within the body. If we write
$\mathbf{r=R+r^{\prime}}$, where $\mathbf{R}$ locates the center of mass of
the body, we can read off the force on the center of mass of the body from
$\tau_{z}=XF_{y}-YF_{x}$, so that

$F_{x}=\int(d\mathbf{r})\int_{0}^{\infty}\frac{d\omega}{2\pi}\left[\frac{1}{e^{\beta\omega}-1}-\frac{1}{e^{\beta^{\prime}\omega}-1}\right]\hat{\chi}_{xz}(\mathbf{r};\omega)\int\frac{(d\mathbf{k}_{\perp})}{(2\pi)^{2}}k^{2}\operatorname{Im}\left[r^{H}e^{-2\kappa
z}\right].$ (46)

We will derive this result directly in Sec. VIII.

###  VI.1 Torque for a nanoparticle above a perfectly conducting plate

A very simple example is provided by a perfectly conducting surface lying in
the $z=0$ plane. This means $r^{H,E}=\pm 1$. Consider only the case with
$\hat{\chi}_{xy}\neq 0$, that is, for our model, the magnetic field lying in
the $z$ direction. The imaginary part comes only from the region where
$\omega^{2}>k^{2}$, where

$\operatorname{Im}\kappa=-\operatorname{sgn}(\omega)\sqrt{\omega^{2}-k^{2}},$
(47)

and then the integral in Eq. (44) over transverse wavenumbers is (provided the
body is of negligible extent, a nanoparticle, so $z=a$)

$\displaystyle\frac{1}{(2\pi)^{2}}\int_{0}^{|\omega|}dk\,k\int_{0}^{2\pi}d\theta\left[-\sqrt{\omega^{2}-k^{2}}-\frac{\omega^{2}}{\sqrt{\omega^{2}-k^{2}}}\right]\cos\left[2a\sqrt{\omega^{2}-k^{2}}\right]$
$\displaystyle=$
$\displaystyle-\frac{\omega^{3}}{2\pi}\int_{0}^{1}dy(y^{2}+1)\cos 2\omega ay$
(48) $\displaystyle=$ $\displaystyle-\frac{1}{\pi(2a)^{2}}[u\cos
u+(u^{2}-1)\sin u].$

Here we have defined $u=2\omega a$. Then, we write the scattering part of the
torque in the direction perpendicular to the plate as

$\tau^{s}_{z}=\frac{2}{\pi^{2}}V\omega_{c}\omega_{p}^{2}\eta\int_{0}^{\infty}du\frac{u\cos
u+(u^{2}-1)\sin u}{(u^{2}+(2\eta a)^{2})^{2}}\left[\frac{1}{e^{\beta
u/(2a)}-1}-\frac{1}{e^{\beta^{\prime}u/(2a)}-1}\right].$ (49)

Close to the plate, $2\eta a\ll 1$, $\beta/(2a)\gg 1$, $\beta^{\prime}/(2a)\gg
1$, where the integral is dominated by small values of $u$, for which the
numerator in the integrand approaches

$u\cos u+(u^{2}-1)\sin u\sim\frac{2}{3}u^{3},$ (50)

we obtain precisely the negative of the torque coming from the vacuum
contribution (27). That the total normal torque vanishes as the perfectly
conducting plate is approached is expected, because the tangential electric
field must vanish there.

The total torque is then

$\tau^{\rm vac}_{z}+\tau_{z}^{\rm
s}=\frac{4}{3\pi^{2}}V\omega_{c}\omega_{p}^{2}\eta\int_{0}^{\infty}du\frac{u^{3}}{(u^{2}+(2a\eta)^{2})^{2}}\left[1-\frac{3}{2}\frac{u\cos
u+(u^{2}-1)\sin
u}{u^{3}}\right]\left[\frac{1}{e^{\beta^{\prime}u/(2a)}-1}-\frac{1}{e^{\beta
u/(2a)}-1}\right].$ (51)

This is plotted in Fig. 3.

Figure 3:  Torque (apart from the prefactor in Eq. (51)) as a function of
separation of the nanoparticle $a$ from the perfectly conducting plate in
$\mu$m, for a damping parameter of $\eta=0.035$ eV. The temperatures are taken
to be $T=300$ K and $T^{\prime}=600$ K. The torque vanishes, as expected,
close to the plate, and approaches the vacuum value far from the plate. Most
interesting is the appearance of a very weak maximum at about 14.4 $\mu$m,
just before the decrease to the vacuum torque value, as displayed in the
inset.

##  VII Quantum vacuum force

Let us start by writing the force on a dielectric body on which an electric
field is impressed, writing the fields in terms of their frequency transforms:

$\displaystyle\mathbf{F}$ $\displaystyle=$
$\displaystyle\int(d\mathbf{r})\int\frac{d\omega}{2\pi}\frac{d\nu}{2\pi}e^{-i(\omega+\nu)t}\left\\{-\left[\bm{\nabla}\cdot\mathbf{P}(\mathbf{r};\omega)\right]\mathbf{E}(\mathbf{r};\nu)-i\omega\mathbf{P}(\mathbf{r};\omega)\times\left[\frac{1}{i\nu}\bm{\nabla}\times\mathbf{E}(\mathbf{r};\nu)\right]\right\\}$
(52) $\displaystyle=$
$\displaystyle\int(d\mathbf{r})\int\frac{d\omega}{2\pi}\frac{d\nu}{2\pi}e^{-i(\omega+\nu)t}\left\\{-\left[\bm{\nabla}\cdot\mathbf{P}(\mathbf{r};\omega)\right]\mathbf{E}(\mathbf{r};\nu)\left(1+\frac{\omega}{\nu}\right)-\frac{\omega}{\nu}\mathbf{P}(\mathbf{r};\omega)\cdot(\bm{\nabla})\cdot\mathbf{E}(\mathbf{r};\nu)\right\\},$

where in the second line we integrated spatially by parts. Unlike for the
torque, the total divergence does not contribute. Now we either expand
$\mathbf{P}$ in terms of $\mathbf{E}$, using Eq. (16b), or $\mathbf{E}$ in
terms of $\mathbf{P}$, using Eq. (16a), and then use the fluctuation-
dissipation theorem on the two parts. This yields rather immediately for the
force on the body777This general formula can be derived immediately from the
torque on the center of mass inferred from Eq. (19).

$F_{i}^{\rm EE}+F_{i}^{\rm
PP}=\int\frac{d\omega}{2\pi}(d\mathbf{r})\left[\chi_{jl}(\mathbf{r};\omega)\nabla^{\prime}_{i}(\Im\bm{\Gamma})_{lj}(\mathbf{r,r^{\prime}};\omega)\big{|}_{\mathbf{r^{\prime}=r}}\coth\frac{\beta\omega}{2}+(\Im\bm{\chi})_{jl}(\mathbf{r};\omega)\nabla_{i}\Gamma_{lj}(\mathbf{r,r^{\prime}};-\omega)\big{|}_{\mathbf{r^{\prime}=r}}\coth\frac{\beta^{\prime}\omega}{2}\right].$
(53)

However, for vacuum, Eq. (21) describes the vacuum Green’s dyadic, So, again
in the coincidence limit, it is then clear that the gradient of the Green’s
dyadic vanishes, and thus there is no vacuum force. The conclusion appears to
be opposite to that of Refs. muller ; Reid , but evidently the self-propulsion
found there arises as a second-order effect, with which we will deal in a
later paper. There, nonreciprocity is not required. The only necessary
conditions are that the system be out of equilibrium, and that the body be
extended and inhomogeneous.

##  VIII Transverse force on a nonreciprocal nanoparticle induced by a
dielectric surface

In contrast to the result found in the previous section, a nonreciprocal body
does experience, in first order, a force transverse to another ordinary body,
even when both bodies are at rest, provided they are not in thermal
equilibrium with each other. This was observed in Ref. muller  and more
recently in Refs. gelb ; Khandekar .

This is still described by the formula (53), but requires the scattering part
of the Green’s function. We will consider the second, ordinary body to be a
planar dielectric, with permittivity $\varepsilon(\omega)$, lying in the
halfspace $z<0$, while the nonreciprocal body lies at a distance $z=a$ above
the plane. It is convenient then to introduce a two-dimensional Fourier
transform in the transverse coordinates, $x$ and $y$. Then the force in the
$x$ direction, say, is

$F_{x}=-\int(d\mathbf{r})\frac{d\omega}{2\pi}\frac{(d\mathbf{k}_{\perp})}{(2\pi)^{2}}ik_{x}\left[\chi_{jk}(\mathbf{r};\omega)(\Im\mathbf{g}^{s})_{kj}(z,z;\omega,\mathbf{k}_{\perp})\coth\frac{\beta\omega}{2}-(\Im\bm{\chi})_{jk}(\mathbf{r};\omega)g^{s}_{kj}(z,z;-\omega,\mathbf{k}_{\perp})\coth\frac{\beta^{\prime}\omega}{2}\right].$
(54)

Here, the $s$ superscripts on the reduced Green’s functions represent the
scattering parts, since it is evident that the bulk (vacuum) part does not
contribute, as already demonstrated in the previous section. Now the
integration over $k_{x}$ and $k_{y}$ will vanish except for $g_{xz}$ and
$g_{zx}$.888For the bulk (vacuum) contribution, the force would involve the
symmetric limit $\lim_{z\to z^{\prime}}\operatorname{sgn}(z-z^{\prime})=0$.
Hence, unlike for the torque, only the TM Green’s function contributes. Using
the properties of $\Im\bm{\chi}$ given in Sec. II, we immediately obtain

$F_{x}=2\int_{0}^{\infty}\frac{d\omega}{2\pi}\int\frac{(d\mathbf{k}_{\perp})}{(2\pi)^{2}}(d\mathbf{r})\hat{\chi}_{xz}(\mathbf{r};\omega)k_{x}^{2}\operatorname{Im}\left(r^{H}e^{-2\kappa
z}\right)\left[\frac{1}{e^{\beta\omega}-1}-\frac{1}{e^{\beta^{\prime}\omega}-1}\right].$
(55)

This force is precisely that inferred from the torque in Eq. (46).

As with Casimir friction, this force will vanish unless dissipation occurs
somewhere. This could be due to dissipation in the dielectric slab, or to
radiation. We will consider these in the following subsections.

###  VIII.1 Dissipation in a metallic slab

We will describe the metallic substrate by a Drude model,

$\varepsilon(\omega)=1-\frac{\omega^{2}_{p}}{\omega^{2}+i\omega\nu},$ (56)

where $\omega_{p}$ is the plasma frequency, and $\nu$ the damping parameter.
For simplicity, we will consider the regime

$\nu\ll\omega\ll\omega_{p},\quad\omega\ll k,$ (57)

so that

$\operatorname{Im}\varepsilon(\omega)\approx\frac{\omega_{p}^{2}\nu}{\omega^{3}},$
(58)

and then, for low frequencies,

$\operatorname{Im}r^{H}=\operatorname{Im}\frac{\kappa-\kappa^{\prime}/\varepsilon}{\kappa+\kappa^{\prime}/\varepsilon}\approx\operatorname{Im}\frac{\varepsilon-1}{\varepsilon+1}\approx\frac{2\omega\nu}{\omega_{p}^{2}}.$
(59)

Thus, in this approximation, where we crudely replace $\kappa$ and
$\kappa^{\prime}$ by $k$, the force on a nanoparticle of negligible extent is

$\displaystyle F_{x}$ $\displaystyle=$
$\displaystyle-2\frac{\nu}{\omega_{p}^{2}}\int\frac{d\omega}{2\pi}\frac{(d\mathbf{k}_{\perp})}{(2\pi)^{2}}\hat{\alpha}_{xz}(\omega)\omega
k_{x}^{2}e^{-2\kappa
a}\left(\frac{1}{e^{\beta\omega}-1}-\frac{1}{e^{\beta^{\prime}\omega}-1}\right)$
(60) $\displaystyle\approx$
$\displaystyle-\frac{12}{(2\pi)^{2}}\frac{1}{(2a)^{4}}\frac{\nu}{\omega_{p}^{2}}\int_{0}^{\infty}d\omega\,\hat{\alpha}_{xz}(\omega)\omega\left(\frac{1}{e^{\beta\omega}-1}-\frac{1}{e^{\beta^{\prime}\omega}-1}\right).$

Now for the nonreciprocal polarizability, let us use the model (13), where we
now assume that the magnetic field (confined to the particle) lies in the $y$
direction. This leads directly to the following formula for the force,
($x=\omega/\eta$)

$F_{x}=\frac{3V}{4\pi^{2}a^{4}}\frac{\nu}{\eta}\omega_{c}f(\beta\eta,\beta^{\prime}\eta),\quad
f(\beta\eta,\beta^{\prime}\eta)=-\int_{0}^{\infty}dx\frac{x}{(x^{2}+1)^{2}}\left[\frac{1}{e^{\beta\eta
x}-1}-\frac{1}{e^{\beta^{\prime}\eta x}-1}\right],$ (61)

where the integral follows from Appendix A,

$\int_{0}^{\infty}dx\frac{x}{(x^{2}+1)^{2}}\frac{1}{e^{\beta\eta
x}-1}=\frac{1}{4}\left[1+\frac{\pi}{\beta\eta}-\frac{\beta\eta}{2\pi}\psi^{\prime}\left(\frac{\beta\eta}{2\pi}\right)\right].$
(62)

The dimensionless force $f$ is plotted in Fig. 4, and compared to the high-
temperature approximation. The prefactor, for $\nu=\eta$, $a=1\,\mu$m, and the
radius of the nanosphere being $100$ nm, is $5\times 10^{-21}$ N.

Figure 4:  Force between a nonreciprocal nanoparticle and a metal plate with
finite conductivity out of thermal equilibrium. The temperature of both the
plate and the background electromagnetic field is fixed at room temperature,
which corresponds to $T/\eta=0.714$ for gold. The polarizability of the
nanoparticle is described by the model (13). The temperature of the
nanoparticle, $T^{\prime}$, is given in units of the damping parameter for
gold, $\eta$. For comparison, the straight line shows the force when both
temperatures are large.

###  VIII.2 Transverse force in presence of a perfectly conducting plate

If the slab is a perfect conductor, with $r^{H}=1$, the formula for the
transverse force simplifies considerably. The imaginary part of the Green’s
function then requires that $\omega^{2}>k^{2}$, and so the integral over the
wavenumber is

$\int\frac{(d\mathbf{k}_{\perp})}{(2\pi)^{2}}k_{x}^{2}\operatorname{Im}e^{-2\kappa
a}=-\frac{1}{4\pi}\frac{1}{(2a)^{4}}[6u\cos u+2(u^{2}-3)\sin u],$ (63)

where $u=2\omega a$. Then, the transverse force is

$F_{x}=\frac{\omega_{c}\eta\omega_{p}^{2}V}{2\pi^{2}a}f_{0}(\epsilon,b,b^{\prime}),\quad
f_{0}(\epsilon,b,b^{\prime})=\int_{0}^{\infty}du\frac{6u\cos u+2(u^{2}-3)\sin
u}{(u^{2}+\epsilon^{2})^{2}}\left[\frac{1}{e^{ub}-1}-\frac{1}{e^{ub^{\prime}}-1}\right],$
(64)

where $\epsilon=2\eta a$, $b=1/(2aT)$, and $b^{\prime}=1/(2aT^{\prime})$. The
integral $f$ is plotted in Fig. 5 as a function of the nanoparticle
temperature $T^{\prime}$ for the environment at room temperature, for a
separation of $a=100$ nm, with a damping parameter appropriate for gold,
$\eta=0.035$ eV. Note, that the high-temperature limit for the force is given
by $f_{0}\sim\frac{\pi}{8b^{\prime}}$, for $b^{\prime}\ll 1/\epsilon,b$, but
this limit requires very high temperatures which are not accessible in
practice. It is noteworthy that Figs. 4 and 5 are qualitatively (but not
quantitatively) similar, given that the physical mechanisms invoked are rather
different.

Figure 5:  The transverse force given in Eq. (64) as a function of the
temperature of the nonreciprocal nanoparticle relative to room temperature,
300 K, for the environment and perfectly conducting plate at room temperature.
Here we take the separation $a$ of the nanoparticle and the plate to be 100
nm, and the damping to be that appropriate for gold, 0.035 eV. In this case,
the prefactor in the force in Eq. (64), for a gold nanosphere of 10 nm radius
is $1.2\times 10^{-20}$ N, so this would be challenging to observe.

It is easily seen that the lateral force rapidly vanishes as $a\to\infty$,
consistent with the absence of a quantum vacuum force.

We expect, as we saw for the torque, that this force will be resisted by the
quantum vacuum friction in the presence of the plate, which for low velocities
will lead to a terminal velocity, according to

$m\frac{dv}{dt}=F_{0}-vF_{1}^{\prime}\Rightarrow
v(t)=\frac{F_{0}}{F_{1}^{\prime}}\left(1-e^{-F_{1}^{\prime}t/m}\right).$ (65)

We require, then, the nonequilibrium frictional force in the presence of a
conducting plate, which we derive in Appendix B. Using the same model for the
permittivity of the nanoparticle, the linear term in the friction is

$F_{1}^{\prime}=\frac{\omega_{p}^{2}\eta
V}{\pi^{2}(2a)^{2}}f_{1}(\epsilon,b,b^{\prime}),$ (66)

where

$\displaystyle f_{1}(\epsilon,b,b^{\prime})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}du\frac{u^{3}}{u^{2}+\epsilon^{2}}\bigg{\\{}\left[1-\frac{2\cos
u+(u^{2}-2)\sin
u}{u^{3}}\right]\left[\frac{1}{e^{b^{\prime}\\!u}-1}-\frac{1}{e^{bu}-1}\right]$
(67)
$\displaystyle\quad\mbox{}+\frac{1}{12}\frac{bu}{\sinh^{2}(bu/2)}\left[1-3\frac{-u(u^{2}-12)\cos
u+(5u^{2}-12)\sin u}{u^{5}}\right]\bigg{\\}}.$

Indeed, $f_{1}$ is always positive, corresponding to a frictional drag, and
the corresponding terminal velocity is

$v_{T}=\frac{F_{0}}{F_{1}^{\prime}}=2\omega_{c}a\frac{f_{0}}{f_{1}}.$ (68)

The scale factor here is small compared to the speed of light: for a particle
100 nm above the plate, $2\omega_{c}a=10^{-4}$. The ratio of $f_{0}/f_{1}$ is
shown in Fig. 6.

Figure 6:  Terminal velocity of a nonreciprocal nanoparticle near a perfectly
conducting plate in units of $2\omega_{c}a$. Here it is assumed that the plate
and the background are at temperature $T=300$ K, that the particle is made of
gold, and it is $a=100$ nm above the surface of the plate. The temperature of
the nanoparticle is $T^{\prime}=1/(2ab^{\prime})$. Thus, the highest particle
temperature displayed on the graph is 1200 K.

The apparent saturation of the terminal velocity near 0.2 is illusory; for
still larger temperatures, the terminal velocity tends to zero, since the
frictional force rapidly increases with temperature. However, for these
nominal values, the damping time is long:

$t_{0}=\frac{m}{F_{1}^{\prime}}\sim 2\times
10^{3}\,\mbox{s}\frac{1}{f_{1}}\sim 10^{6}\,\mbox{s},$ (69)

if the particle is at twice room temperature.

##  IX Relaxation to Thermal Equilibrium

All of the above considerations assumed that the temperatures of the body and
of the background are constant. Of course, this will not be so unless some
mechanism keeps the system out of thermal equilibrium. Here, we will calculate
the time it would take for such a body at rest to come to thermal equilibrium
with its environment. We cannot regard the body to be a black body, but we can
calculate the rate at which it loses heat from the power (for an isotropic
body)999For the purpose of the rough calculation presented here, we will
ignore the small nonreciprocal effects. guo2

$P(T,T^{\prime})=\frac{1}{\pi^{2}}\int_{0}^{\infty}d\omega\,\omega^{4}\operatorname{Im}\alpha(\omega)\left[\frac{1}{e^{\beta\omega}-1}-\frac{1}{e^{\beta^{\prime}\omega}-1}\right]=\frac{dQ}{dt}.$
(70)

This is related to the rate of change of temperature of the body by its heat
capacity:

$\frac{dQ}{dt}=C_{V}(T^{\prime})\frac{dT^{\prime}}{dt}.$ (71)

Thus, the time it takes for the body to cool from temperature $T^{\prime}_{0}$
to temperature $T^{\prime}_{1}$, where $T^{\prime}_{0}>T^{\prime}_{1}>T$, is

$t=\int_{T^{\prime}_{0}}^{T^{\prime}_{1}}dT^{\prime}\frac{C_{V}(T^{\prime})}{P(T,T^{\prime})}.$
(72)

To proceed, we need a model for the heat capacity of the body. Such is
provided by the Debye model,101010We consider bulk effects only, and ignore
surface effects, for the purpose of a rough estimate. which is satisfactory
for simple crystals (see Ref. LL ):

$C_{V}(T)=9N\\!\left(\frac{T}{\Theta}\right)^{\\!3}\\!\\!\int_{0}^{\Theta/T}\\!\\!dx\frac{x^{4}e^{x}}{(e^{x}-1)^{2}}.$
(73)

where $N$ is the number of atoms constituting the body, and $\Theta$ is the
Debye temperature. This interpolates between the low- and high-temperature
limits:

$C_{V}(T)\sim 3N\left\\{\begin{array}[]{cc}1,&T\gg\Theta,\\\
\frac{4\pi^{4}}{5}\\!\\!\left(\frac{T}{\Theta}\right)^{\\!3},&T\ll\Theta.\end{array}\right.$
(74)

Since the Debye temperature for gold is about $\Theta=170$ K, the high-
temperature approximation would seem appropriate for an estimate at room
temperature and above.

We finally need a model for the imaginary part of the polarizability of the
body. The Lorenz-Lorentz model would give

$\operatorname{Im}\alpha(\omega)=\frac{V\omega_{p}^{2}\omega\eta}{(\omega_{1}^{2}-\omega^{2})^{2}+\omega^{2}\eta^{2}}\approx\frac{V\omega_{p}^{2}\omega\eta}{\omega_{1}^{4}},$
(75)

where, for a metal (Drude model), $\omega_{1}=\omega_{p}/\sqrt{3}$. The
approximation here is appropriate if, as expected, $\omega_{1}\gg\omega,\eta$.
Inserting this approximation into the formula (70) we obtain

$P(T,T^{\prime})\approx\frac{8\pi^{4}}{7}\frac{V\eta}{\omega_{p}^{2}}(T^{6}-T^{\prime
6}).$ (76)

Now we compute the cooling time from Eq. (72):

$t=t_{0}\int_{T^{\prime}_{0}/T}^{T^{\prime}_{1}/T}\\!du\frac{1}{1-u^{6}},\quad
T^{\prime}_{0}>T^{\prime}_{1}>T,\quad\mbox{where}\quad
t_{0}=\frac{21}{8\pi^{4}}n\frac{\omega_{p}^{2}}{\eta}\frac{1}{T^{5}}$ (77)

Here, $n$ is the number density of atoms in the body. The relaxation scale,
$t_{0}$, is independent of the volume of the particle, and is about $10^{4}$
seconds for gold, for an environmental temperature of 300 K. The cooling time
diverges as $T^{\prime}_{1}\to T$, but cooling to a temperature slightly above
the environmental temperature takes a finite time. The integral here is
elementary, but the resulting expression is not very illuminating. We content
ourselves by showing some representative values in left panel of Fig. 7.

Figure 7: Time required for a body to cool from temperature $T^{\prime}_{0}$
to temperature $T^{\prime}_{1}$ for different environmental temperatures $T$.
Here $T^{\prime}_{0}>T^{\prime}_{1}>T$. The left panel is for the
environmental temperature $T=300$ K, and the right is for $T=1$ K. On the
left, the upper set of curves is for $T^{\prime}_{1}/T=1.1$, and the lower
curves are for $T^{\prime}_{1}/T=1.5$. On the right, the three curves from top
to bottom are for $T^{\prime}_{1}/T=1.05$, 1.1, and 1.5. The times are scaled
by the prefactor, $t_{0}$, in Eq. (77), which for a gold body evaluates to
about $10^{4}$ s, for the environment at room temperature, and by
$\tilde{t}_{0}$ (Eq. (79)), which is about $10^{11}$ s, for $T=1$ K.

It will be seen that if $T^{\prime}_{0}$ is appreciably larger than $T$ the
cooling time rapidly saturates to an asymptotic value. If we then take
$T^{\prime}_{0}$ to be large, Fig. 8 shows how long it will take to reach a
multiple of the environmental temperature.

Figure 8:  The cooling time, in terms of the scale $t_{0}$, for the
nanoparticle to cool from a high temperature to $T_{1}$. It takes an
increasingly long time to get very close to the environmental temperature.

Thus, we see that the terminal angular velocity seen in Eq. (40) and the
terminal linear velocity obtained in Eq. (68) will not be achievable unless
some mechanism maintains the thermal imbalance, because the time scales for
achieving those velocities, Eqs. (41) and (69), are much longer than the
cooling time found here.

If the environmental temperature is very low, $T\ll\Theta$, the cooling time
is very much longer. The analysis proceeds as above, using the low temperature
limit in Eq. (74), with the result for the cooling time being

$t=\tilde{t}_{0}\int_{(T^{\prime}_{0}/T)^{2}}^{(T^{\prime}_{1}/T)^{2}}\\!\\!dy\frac{y}{1-y^{3}},\quad\tilde{t}_{0}=\frac{21}{20}n\frac{\omega_{p}^{2}}{\eta}\left(\frac{T}{\Theta}\right)^{\\!3}\\!\frac{1}{T^{5}}.$
(78)

The integral, which has a relatively simple analytic form, is shown in the
right panel of Fig. 7. The ratio of the time scales in the two cases is

$\frac{\tilde{t}_{0}}{t_{0}}=\frac{2\pi^{4}}{5}\left(\frac{T_{\rm
low}}{\Theta}\right)^{\\!\\!3}\left(\frac{T_{\rm high}}{T_{\rm
low}}\right)^{\\!\\!5}\sim 10^{7},$ (79)

for $T_{\rm low}=1$ K, $T_{\rm high}=300$ K, and $\Theta=170$ K for gold, so
terminal velocities might be achievable.

##  X Lorenz-Lorentz correction

Hithertofore, we have ignored the Lorenz-Lorentz correction familiar in
passing from the permittivity of a body to its polarizability. This was
because the forces and torques were derived directly from the macroscopic
susceptibilities appropriate for a dissipative metal body. In the case of
small bodies, we could always pass from the susceptibility to the mean
polarizability by integrating over the volume of the body. However, as is
evident from the discussion in the preceding section [see Eq. (75)], the
effect of the medium on the local electric field can result in a large
correction in the case of metal bodies.

The difficulty is that the simple Lorenz-Lorentz model is ordinarily derived
using spherical symmetry. The relation between the electric polarizability and
the permittivity is, in the isotropic case, in HL units,

$\alpha=\frac{\varepsilon-1}{\varepsilon+2}4\pi a^{3}.$ (80)

This is not valid for a nonsymmetric permittivity; for example, see Ref. Pal .
However, for $\omega_{c}$ small, the nonsymmetric nature is small, so, for the
purpose of an estimate, we use, as in Refs. fan ; strekha , the matrix
generalization of the above:

$\bm{\alpha}=(\bm{\varepsilon-1})(\bm{\varepsilon+2})^{-1}4\pi a^{3}.$ (81)

It is quite straightforward to compute the components of this matrix: the term
we need for the torque in Eq. (26) is, for $\omega_{p}\gg\omega\sim T$,

$\operatorname{Re}\alpha_{xy}\approx
54V\frac{\omega^{2}\omega_{c}\eta}{\omega_{p}^{4}}.$ (82)

When this is inserted into Eq. (26), we obtain

$\tau_{z}=\frac{32}{7}\pi^{4}V\frac{\omega_{c}\eta}{\omega_{p}^{4}}T^{6}\left[1-\left(\frac{T^{\prime}}{T}\right)^{\\!\\!6}\right].$
(83)

Putting in the numbers for a 100 nm gold nanosphere, the coefficient of
$(1-T^{\prime 6}/T^{6})$ is about $5\times 10^{-36}$ Nm, some 11 orders of
magnitude smaller than that found at the end of Sec. IV.

We can also repeat the calculation of the terminal angular velocity in Sec. V
in this Lorenz-Lorentz model. Then,

$\tau^{\prime}_{1}=-\frac{1}{6\pi^{2}}\left(3+\beta\frac{\partial}{\partial\beta}\right)\int_{0}^{\infty}d\omega\,\omega^{2}\operatorname{Im}(\alpha_{xx}+\alpha_{yy})(\omega)\left(\frac{1}{e^{\beta^{\prime}\\!\omega}-1}-\frac{1}{e^{\beta\omega}-1}\right),$
(84)

where in our model

$\operatorname{Im}(\alpha_{xx}+\alpha_{yy})\approx
18V\frac{\omega\eta}{\omega_{p}^{2}},$ (85)

which implies

$\tau_{1}^{\prime}=\frac{2\pi^{2}}{5}\frac{V\eta}{\omega_{p}^{2}}T^{4}\left[1+3\left(\frac{T^{\prime}}{T}\right)^{\\!\\!4}\right].$
(86)

Note that $\tau_{1}^{\prime}$ is always positive, indicating that it always
opposes the rotation. The corresponding terminal angular velocity is

$\Omega_{T}=\frac{\tau_{0}}{\tau_{1}^{\prime}}=\frac{80}{7}\pi^{2}\frac{\omega_{c}}{\omega_{p}^{2}}T^{2}\frac{1-\frac{T^{\prime
6}}{T^{6}}}{1+3\frac{T^{\prime 4}}{T^{4}}},$ (87)

where the prefactor, independent of $T^{\prime}$, implies a substantial
angular velocity for gold at room temperature: $\sim 10^{8}$ s-1. Although the
time required to reach such a velocity is very long,
$t_{0}=I/\tau_{1}^{\prime}\sim 10^{13}$ s, the initial angular acceleration is
not so small,

$\dot{\Omega}(0)=\frac{\Omega_{T}}{t_{0}}\sim 10^{-5}\,\mbox{s}^{-2}.$ (88)

While this angular acceleration is 10 orders of magnitude smaller than that
found without the Lorenz-Lorentz correction in Eq. (42), the body will acquire
a measurable angular velocity after a relatively small period of observation.

However, is this correction valid or even necessary for a metal nanoparticle?
The discussion in Secs. IV–VIII is based on describing the susceptiblity of a
metal by the phenomenological Drude model, which should include,
approximately, all internal effects. There is a large literature on the
subject of ordinary polarizabilities of metal nanoparticles—see Refs. Puska ;
Kheirandish , for example—where it is seen that both classical and approximate
quantum mechanical treatments are inadequate. We are unaware of comparable
work in the nonreciprocal case. So, to some extent, the issue of applying the
Lorenz-Lorentz correction remains open. In this paper we are interested in the
interaction between the electromagnetic field and the body, the
electromagnetic properties of which are specified by a given susceptibility,
so that the crude models for the latter should only be taken as illustrative.

##  XI Conclusions

In this paper we have concentrated on analysis to first order in the
susceptibility, to better understand the effects of nonreciprocal materials on
torque and on forces for bodies out of thermal equilibrium with their
environment. Time-reversal symmetry is broken by these materials, so
spontaneous forces and torques are possible. Of course, time reversal symmetry
is not broken by electrodynamics, whether classical or quantum; rather the
nonreciprocity is a consequence of an external agent, such as a magnetic
field, that is encoded in the dielectric response of the materials.

Interestingly, potentially observable phenomena are, nevertheless, predicted.
A nonreciprocal body out of thermal equilibrium will spontaneously start to
rotate, and reach a substantial terminal angular velocity. Such a body will
not feel a net force to first order in the susceptibility. However, if the
body is placed near a translationally invariant surface, even a perfect
conductor, then a force parallel to the surface would arise. The presence of
such a surface would tend to suppress the vacuum torque. A potentially
observable terminal linear velocity arises here as well, although the time
scales are such that it would be difficult to keep the system out of thermal
equilibrium. A possible drastic reduction in the strength of these
nonreciprocal effects, due to the Lorenz-Lorentz correction for dielectric
susceptibilities, is discussed in the penultimate section, although it seems
the angular and linear accelerations might still be amenable to observation.

Elsewhere, we will examine higher-order effects, to see how phenomena such as
vacuum self-propulsion can arise, even for a reciprocal body Reid ; muller .

###### Acknowledgements.

The work of KAM and XG was supported in part by a grant from the US National
Science Foundation, No. 2008417. The work of GK and KAM was supported in part
by a grant from the US National Science Foundation, No. PHY-1748958. In
particular, GK and KAM thank KITP for its hospitality, and Benjamin Strekha
for conversations there, which stimulated the research carried out here. We
thank Steve Fulling, Li Yang, Prachi Parashar, Shadi Rezaei, and Venkat
Abhignan for collaborative assistance. This paper reflects solely the authors’
personal opinions and does not represent the opinions of the authors’
employers, present and past, in any way. For the purpose of open access, the
authors have applied a CC BY public copyright license to any Author Accepted
Manuscript version arising from this submission.

##  Appendix A Evaluation and expansion of integrals

Here, we express the integrals encountered in the text in terms of the digamma
and trigamma functions, $\psi(z)$ and $\psi^{\prime}(z)$, and provide
corresponding low- and high-temperature expansions.

Differentiation of Binet’s second integral representation of the log gamma
function immediately yields the integral representation

$\psi(z)=\ln
z-\frac{1}{2z}-2\int_{0}^{\infty}dx\,\frac{x}{x^{2}+1}\frac{1}{e^{2\pi
zx}-1}.$ (89)

Thus,

$I_{1}(\beta\eta)\equiv\int_{0}^{\infty}dx\,\frac{x}{x^{2}+1}\frac{1}{e^{\beta\eta
x}-1}=\frac{1}{2}\left[-\frac{\pi}{\beta\eta}+\ln\left(\frac{\beta\eta}{2\pi}\right)-\psi\left(\frac{\beta\eta}{2\pi}\right)\right].$
(90)

Since

$\beta\frac{\partial}{\partial\beta}\frac{1}{e^{\beta\eta
x}-1}=\eta\frac{\partial}{\partial\eta}\frac{1}{e^{\beta\eta
x}-1}=x\frac{\partial}{\partial x}\frac{1}{e^{\beta\eta x}-1},$ (91)

it follows that

$\beta\frac{\partial}{\partial\beta}I_{1}(\beta\eta)=\eta\frac{\partial}{\partial\eta}I_{1}(\beta\eta)=\int_{0}^{\infty}dx\,\frac{x^{2}}{x^{2}+1}\frac{\partial}{\partial
x}\frac{1}{e^{\beta\eta
x}-1}=2\int_{0}^{\infty}dx\,\left[\frac{x^{3}}{(x^{2}+1)^{2}}-\frac{x}{x^{2}+1}\right]\frac{1}{e^{\beta\eta
x}-1}.$ (92)

Thus,

$I_{2}(\beta\eta)\equiv\int_{0}^{\infty}dx\,\frac{x^{3}}{(x^{2}+1)^{2}}\frac{1}{e^{\beta\eta
x}-1}=\left[1+\frac{\beta}{2}\frac{\partial}{\partial\beta}\right]I_{1}(\beta\eta)=\left[1+\frac{\eta}{2}\frac{\partial}{\partial\eta}\right]I_{1}(\beta\eta).$
(93)

Hence, from Eq. (90),

$I_{2}(\beta\eta)=\frac{1}{4}\left[-\frac{\pi}{\beta\eta}+2\ln\left(\frac{\beta\eta}{2\pi}\right)+1-2\psi\left(\frac{\beta\eta}{2\pi}\right)-\frac{\beta\eta}{2\pi}\psi^{\prime}\left(\frac{\beta\eta}{2\pi}\right)\right].$
(94)

Using the series representation

$\psi(z)=-\frac{1}{z}-\gamma_{E}+z\sum_{k=1}^{\infty}\frac{1}{k(z+k)},$ (95)

where $\gamma_{E}$ is the Euler-Mascheroni constant, we readily obtain from
Eqs. (90) and (94) the small $\beta\eta$ (or high-temperature) expansions

$I_{1}(\beta\eta)\sim\frac{1}{2}\left[\frac{\pi}{\beta\eta}+\ln\left(\frac{\beta\eta}{2\pi}\right)+\gamma_{E}\right]\quad\quad(\beta\eta\to
0)$ (96a) and
$I_{2}(\beta\eta)\sim\frac{1}{4}\left[\frac{\pi}{\beta\eta}+2\ln\left(\frac{\beta\eta}{2\pi}\right)+1+2\gamma_{E}\right]\quad\quad(\beta\eta\to
0).$ (96b)

Likewise, using the asymptotic representation

$\psi(z)\sim\ln
z-\frac{1}{2z}-\sum_{k=1}^{\infty}\frac{B_{2k}}{2kz^{2k}}\quad\quad(z\to\infty),$
(97)

there follow from Eqs. (90) and (94) the large $\beta\eta$ (or low-
temperature) expansions

$I_{1}(\beta\eta)\sim\frac{\pi^{2}}{6\beta^{2}\eta^{2}}-\frac{\pi^{4}}{15\beta^{4}\eta^{4}}\quad\quad(\beta\eta\to\infty)$
(98a) and
$I_{2}(\beta\eta)\sim\frac{\pi^{4}}{15\beta^{4}\eta^{4}}\quad\quad(\beta\eta\to\infty).$
(98b)

##  Appendix B Out-of-equilibrium frictional force near perfectly conducting
plate

Following the discussion in Ref. guo2 , it is easy to derive the expression
for the frictional force used in Sec. VIII.2. The general expression for the
force is

$F=\int\frac{d\omega}{2\pi}\frac{d\mathbf{k}_{\perp}}{(2\pi)^{2}}(k_{x}+\omega
v)\operatorname{tr}\Im\bm{\alpha}(\omega)\Im\mathbf{g^{\prime}}(\omega,\mathbf{k}_{\perp})\left[\coth\frac{\beta\gamma(\omega+k_{x}v)}{2}-\coth\frac{\beta^{\prime}\omega}{2}\right].$
(99)

Here $\mathbf{g^{\prime}}$ is the reduced Green’s function in the rest frame
of the particle. For a perfectly conductimg plate, as with the vacuum,
$\mathbf{g^{\prime}}=\mathbf{g}$, the Green’s function in the frame of the
plate and vacuum. From this we can calculate both the frictional force and the
propulsive force. The latter is present at $v=0$, and arises from the
antisymmetric parts of the polarizability and the Green’s dyadic, and is given
in Sec. VIII. The diagonal parts of both these tensors correspond to friction.
According to the model (10) with $\omega_{c}$ neglected, the diagonal terms of
the imaginary part of the polarizability are all equal to

$\operatorname{Im}\alpha_{d}=\frac{V\omega_{p}^{2}\omega\eta}{(\omega_{0}^{2}-\omega^{2})^{2}+\omega^{2}\eta^{2}}.$
(100)

That leaves us with the imaginary part of the trace of $\mathbf{g}$:

$\operatorname{Im}\operatorname{tr}\mathbf{g}=\operatorname{sgn}(\omega)\left[\frac{\omega^{2}}{\sqrt{\omega^{2}-k^{2}}}-\sqrt{\omega^{2}-k^{2}}\cos\left(2a\sqrt{\omega^{2}-k^{2}}\right)\right]\theta(\omega^{2}-k^{2}).$
(101)

Now when we expand Eq. (99) to first order in $v$ we obtain two terms:

$F=F^{(1)}+F^{(2)}.$ (102)

Here the first term comes from expanding the hyperbolic cotangent:

$F^{(1)}=-\frac{\beta
v}{8\pi^{2}}\int_{0}^{\infty}d\omega\operatorname{Im}\alpha_{d}(\omega)\frac{\omega^{5}}{\sinh^{2}(\beta\omega/2)}\left\\{\frac{2}{3}-\frac{2}{x^{5}}\left[-x(x^{2}-12)\cos
x+(5x^{2}-12)\sin x\right]\right\\},$ (103)

where we’ve carried out the elementary integration over $\mathbf{k}_{\perp}$.
Note that the first term here corresponds to the usual Einstein-Hopf effect.
The $F^{(2)}$ contribution to the force corresponds to the $\omega v$
prefactor in Eq. (99), and is a non-equilibrium friction contribution:

$F^{(2)}=\frac{v}{2\pi^{2}}\int_{0}^{\infty}d\omega\operatorname{Im}\alpha_{d}(\omega)\omega^{4}\left(\coth\frac{\beta\omega}{2}-\coth\frac{\beta^{\prime}\omega}{2}\right)\left[1-\frac{1}{x^{3}}\left(2x\cos
x+(x^{2}-2)\sin x\right)\right],$ (104)

again after carrying out the wavenumber integration. The sum of these two
terms is Eq. (67) if we set $\omega_{0}=0$, appropriate for a metal
nanoparticle.

## References

  * (1) E. V. Teodorovich, Contribution of macroscopic van der Waals interactions to frictional force. Proc. R. Soc. Lond. A 362, 71–77 (1978).
  * (2) L. S. Levitov, Van der Waals friction. Europhys. Lett. 8, 499–504 (1989).
  * (3) J. S. Høye, I. Brevik and K. A. Milton, The reality of Casimir friction [arXiv:1508.00626], Symmetry 8 (5), 29 (2016), doi:10.3390/sym8050029.
  * (4) J. B. Pendry, Shearing the vacuum—quantum friction, J. Phys. Condens. Matter 9, 10301 (1997).
  * (5) A. I. Volokitin and B. N. J. Persson, Theory of friction: The contribution from a fluctuating electromagnetic field, J. Phys. Condens. Matter 11, 345 (1999).
  * (6) G. V. Dedkov and A. A. Kyasov, The relativistic theory of fluctuation electromagnetic interactions of moving neutral particles with a flat surface, Phys. Solid State 45, 1815 (2003).
  * (7) P. C. W. Davies, Quantum vacuum friction, J. Opt. B 7, S40 (2005).
  * (8) G. Barton, On van der Waals friction. II: Between atom and half-space, New J. Phys. 12, 113045 (2010).
  * (9) F. Intravaia, V. E. Mkrtchian, S. Y. Buhmann, S. Scheel, D. A. R. Dalvit, and C. Henkel, Friction forces on atoms after acceleration, J. Phys. Condens. Matter 27, 214020 (2015).
  * (10) M. Oelschläger, Fluctuation-induced phenomena in nanophotonic systems, Ph.D. thesis, Institut für Physik, Humboldt Universität zu Berlin, 2020.
  * (11) A. Einstein and L. Hopf, Statistische Untersuchung der Bewegung eines Resonators in einem Strahlungsfeld, Ann. Phys. (Leipzig) 338, 1105 (1910).
  * (12) X. Guo, K. A. Milton, G. Kennedy, W. P. McNulty, N. Pourtolami, and Y. Li, Energetics of quantum vacuum friction: Field fluctuations, [arXiv:2108.01539] Phys. Rev. D 104, 116006 (2021).
  * (13) X. Guo, K. A. Milton, G. Kennedy, W. P. McNulty, N. Pourtolami, and Y. Li, Energetics of quantum vacuum friction. II: Dipole fluctuations and field fluctuations, [arXiv:2204.10886], Phys. Rev. D 106, 016008 (2022).
  * (14) D. Reiche, F. Intravaia, J.-T. Hsiang, K. Busch, and B.-L. Hu, Nonequilibrium thermodynamics of quantum friction, Phys. Rev. A 102, 050203(R) (2020).
  * (15) G. V. Dedkov and A. A. Kyasov, Tangential force and heating rate of a neutral relativistic particle mediated by equilibrium background radiation, Nucl. Instrum. Methods Phys. Res., Sect. B 268, 599 (2010).
  * (16) G. Pieplow and C. Henkel, Fully covariant radiation force on a polarizable particle, New J. Phys. 15, 023027 (2013).
  * (17) A. I. Volokitin and B. N. J. Persson, Electromagnetic Fluctuations at the Nanoscale (Springer, Berlin, 2017).
  * (18) M. Krüger, T. Emig, and M. Kardar, Nonequilibrium Electromagnetic Fluctuations: Heat Transfer and Interactions, Phys. Rev. Lett. 106, 210404 (2011).
  * (19) A. Ott, P. Ben-Abdallah, and S.-A. Biehs, Circular heat and momentum flux radiated by magneto-optical nanoparticles, Phys. Rev. B 97, 205414 (2018).
  * (20) M. F. Maghrebi, A. V. Gorshkov, and J. D. Sou, Fluctuation-Induced Torque on a Topological Insulator out of Thermal Equilibrium, Phys. Rev. Lett. 123, 055901 (2019).
  * (21) C. Khandekar and Z. Jacob, Thermal spin photonics in the near-field of nonreciprocal media, New J. Phys. 21, 103030 (2019)
  * (22) C. Khandekar, S. Buddhiraju, P. R. Wilkinson, J. K. Gimzewski, A. W. Rodriguez, C. Chase, S. Fan, Nonequilibrium lateral force and torque by thermally excited nonreciprocal surface electromagnetic waves, Phys. Rev. B 104, 245443 (2021).
  * (23) C. Guo and S. Fan, Theoretical constraints on reciprocal and non-reciprocal many-body radiative heat transfer, [arXiv:2007.13274] Phys. Rev. B 102, 085401 (2020).
  * (24) Y. Guo and S. Fan, A single gyrotropic particle as a heat engine, [arXiv:2007.11234] ACS Photonics 8, 1623–1629 (2021).
  * (25) B. Strekha, S. Molesky, P. Chao, M. Krüger, and A. Rodriguez, Trace expressions and associated limits for non-equilibrium Casimir torque, [arXiv:2207.13646] Phys. Rev. A 106, 042222 (2022).
  * (26) J. Schwinger, L. L. DeRaad, Jr., K. A. Milton, and W.-y. Tsai, Classical Electrodynamics (Perseus/Taylor and Francis, New York, 1998).
  * (27) Y. Li, K. A. Milton, P. Parashar, G. Kennedy, N. Pourtolami, and X. Guo, “Casimir self-entropy of nanoparticles with classical polarizabilities: Electromagnetic field fluctuations,” Phys. Rev. D 106, 036002 (2022).
  * (28) D. Pan, H. Xu, and F. J. García de Abajo, Magnetically activated rotational vacuum friction, Phys. Rev. A 99, 062509 (2019).
  * (29) M. T. H. Reid, O. D. Miller, A. G. Polimeridis, A. W. Rodriguez, E. M. Tomlinson, and S. G. Johnson, Photon torpedoes and Rytov pinwheels: Integral-equation modeling of non-equilibrium fluctuation-induced forces and torques on nanoparticles, arXiv:1708.01985 (2017)
  * (30) B. Müller and M. Krüger, Anisotropic particles near surfaces: Propulsion force and friction, Phys. Rev. A 93, 032511 (2016).
  * (31) D. Gelbwaser-Klimovsky, N. Graham, M. Kardar, and M. Krüger, Near Field Propulsion Forces from Nonreciprocal Media, Phys. Rev. Lett. 126, 170401 (2021).
  * (32) L. D. Landau and E. M. Lifshitz, Statistical Physics, Part 1 (Pergamon, Oxford, 1980), 3rd edition, Sec. 66.
  * (33) K. A. Milton, P. Kalauni, P. Parashar, and Y. Li, Casimir self-entropy of a spherical electromagnetic $\delta$-function shell, [arXiv: 1707.09840], Phys. Rev. D 96, 085007 (2017).
  * (34) P. Palffy-Muhoray and D. A. Balzarini, The Clausius-Mossotti relation for anisotropic molecular fluids, Can. J. Phys. 59, 375 (1981).
  * (35) M. J. Puska, R. M. Nieminen, and M. Manninen, Electronic polarizability of small metal spheres, Phys. Rev. B 31, 3486 (1985).
  * (36) A. Kheirandish, N. S. Javan, and H. Mohammadzedah, Modified Drude model for small gold nanoparticles surface plasmon resonancce based on the role of classical confinement, Sci. Rep. (2000) 10:6517 [https://doi.org/10.1038/s41598-020-63066-9].


# Buckling of graphene/MoS2 van der Waals heterostructures: the misfit strain
effect

Run-Sen Zhang  Shanghai Key Laboratory of Mechanics in Energy Engineering,
Shanghai Institute of Applied Mathematics and Mechanics, School of Mechanics
and Engineering Science, Shanghai University, Shanghai 200072, People’s
Republic of China  Jin-Wu Jiang  Corresponding author<EMAIL_ADDRESS>Shanghai Key Laboratory of Mechanics in Energy Engineering, Shanghai Institute
of Applied Mathematics and Mechanics, School of Mechanics and Engineering
Science, Shanghai University, Shanghai 200072, People’s Republic of China

###### Abstract

Van der Waals heterostructures are constructed by stacking different atomic
layers and can inherit many novel electronic and optical properties from the
constituting atomic layers. Mechanical stability is of key importance for the
high performance of nano devices based on the van der Waals heterostructure.
In particular, buckling instability is a critical mechanical issue for the
heterostructure due to its two-dimensional nature. Using graphene/MoS2
heterostructure as an example, the present work demonstrates the relationship
between the buckling instability and the inevitable misfit strain in the
heterostructure by molecular dynamics simulations. The misfit strain has
rather different effects on the buckling phenomenon depending on the magnitude
of the misfit strain. (1). For negative misfit strain, the buckling stability
of the heterostructure is reduced by the misfit strain. It is because the
graphene layer, which initiates the buckling process in the heterostructure of
negative misfit strain, is pre-compressed by the misfit strain that
accelerates the buckling of the graphene layer. (2). For small positive misfit
strain, the buckling stability for the graphene/MoS2 heterostructure is
elevated. The underlying mechanism is that the graphene layer initiates the
buckling process of the heterostructure and is pre-stretched by the small
positive misfit strain, which decelerates the buckling of the graphene layer.
(3). For large positive misfit strain, the graphene layer is pre-stretched
while the MoS2 layer is considerably pre-compressed, so the buckling of the
heterostructure is initiated by the MoS2 layer. As a consequence, the buckling
stability of the graphene/MoS2 heterostructure is reduced by the increase of
large positive misfit strain. These findings shall be valuable for
understanding mechanical properties of van der Waals heterostructures.

van der Waals heterostructure, misfit strain, buckling

##  I Introduction

Two-dimensional (2D) van der Waals heterostructures have attracted great
scientific research interest in recent years due to their outstanding
functional properties inherited from different constitute layers.Novoselov et
al. (2016); Liu et al. (2016) Compared with bulk materials, 2D materials have
ultra low bending rigidity and can be easily buckled or wrinkled under
compressive load in the in-plane direction. Taking advantage of the buckling
or wrinkling properties, 2D materials are promising candidates for stretchable
and flexible electronics, functional surface bionics, soft sensors, actuators,
and etcKo et al. (2020); Ying (2016).

The buckling of single atomic layers like graphene or MoS2 have been
extensively studied. It has been illustrated that buckling properties of
graphene are sensitive to various effects, including the initial ripple,Xiang
and Shen (2016) the defect,Tserpes and Vatistas (2015) the size effect and
biaxial compression,Pradhan and Murmu (2010); Shi et al. (2011) and the shear
loading.Wen et al. (2011); Zhang et al. (2012); Zhou et al. (2015) Jiang
investigated the effect of strain rate and temperature on the critical
buckling strain of the single layer MoS2 sheet under uniaxial
compression.Jiang (2014) The fan-blade shaped wrinkling of single layer MoS2
under circular torsion was studied by Bao et al..Bao et al. (2017) The
intrinsic thermal vibration induced ripples are closely related to the
stiffness and bending modulus for the MoS2.Singh et al. (2015) The defect can
reduce the buckling stability of MoS2 under uniaxial compressive loadings.Li
et al. (2019)

As a characteristic feature of the van der Waals heterostructure, misfit
strain between the constituting atomic layers is inevitable, due to different
lattice constants of the atomic layers.Jiang and Park (2014); Lin et al.
(2014); Li and Zhao (2014) The misfit strain in the heterostructure can be
controlled through various methods, including the substrate,Ni et al. (2008)
the structural design,Chaste et al. (2018) and the photoelectrochemical
etching process.Shivaraman et al. (2013) The misfit strain can cause direct
effects on physical and mechanical properties for the van der Waals
heterostructure. For instance, it was found that the misfit strain induces
additional energy dissipation for the mechanical resonant oscillation of the
graphene/MoS2 heterostructure.He and Jiang (2019)

Although the buckling behavior of individual atomic layers has been
extensively studied, the buckling phenomenon for the van der Waals
heterostructures is not well studied. In particular, the effect of the
inevitable misfit strain on the buckling of van der Waals heterostructure is
still unclear. The mechanical stability, especially the buckling stability, is
a key factor for stable performance of nano devices based on the van der Waals
heterostructure. It is thus an urgent task to investigate the mechanism for
the buckling of heterostructure with different misfit strains.

In this paper, we perform molecular dynamics (MD) simulations to investigate
the buckling behavior of graphene/MoS2 van der Waals heterostructures
subjected to compressive loadings. We find that the misfit strain within the
graphene/MoS2 heterostructure has significant effect on the buckling
phenomenon. There are three typical ranges for the misfit strain. For negative
misfit strain, the graphene layer in the heterostructure is pre-compressed, so
the graphene layer initiates the buckling process of the heterostructure.
Hence, the critical buckling strain for the heterostructure is considerably
reduced by the negative misfit strain. For small positive misfit strain, the
buckling of the heterostructure is still initiated by the graphene layer which
is slightly pre-stretched by the misfit strain, so the critical buckling
strain for the heterostructure is enhanced by small positive misfit strain.
For large positive misfit strain, the MoS2 layer is pre-compressed seriously
and it initiates the buckling process of the heterostructure, so the critical
buckling strain will be reduced by increasing misfit strain in this range. We
also show that the buckling direction for the heterostructure is also affected
by the misfit strain.

##  II Structure and simulation details

Figure 1: Bilayer structures studied in the present work. (a) The
graphene/graphene structure is compressed and buckled. (b) The MoS2/MoS2
structure is compressed and buckled. (c) The graphene/MoS2 heterostructure.
The graphene and MoS2 have a difference of $\Delta l$ in their length,
resulting in a misfit strain in the relaxed configuration. The heterostructure
is compressed and buckled.

In the present work, three structures shown in Fig.1 are comparatively
investigated. The graphene/graphene bilayer structure has the AB stacking
order. The MoS2/MoS2 bilayer structure is in its lowest-energy stacking order.
These structures are compressed by uniaxial compressive strains along the
length direction, and will be buckled when the strain is above a critical
value. The buckling property is sensitive to the length of the structure,
while the width is not important for the buckling phenomenon. We thus choose a
constant width for graphene and MoS2 for all structures simulated in the
present work. The width of graphene is 12.3 Å and the width of MoS2 is 12.5 Å.
These widths are properly chosen, so that the misfit strain in the width
direction is minimum. The length is along the armchair direction, and the
width is along the zigzag direction.

All MD simulations are performed by the Large-scale Atomic/Molecular Massively
Parallel Simulator (LAMMPS) code,Plimpton (1995) and the OVITO package is
employed for visualization.Stukowski (2009) The standard Newton equations of
motion are integrated in time using the velocity Verlet algorithm with a time
step of 1 fs. All simulations are performed at the low temperature of 4.2 K,
so that the effect of thermal vibration is less important. Periodic boundary
conditions are applied in the two in-plane directions, while the free boundary
condition is applied in the out-of-plane direction. The structures are
uniaxially compressed along horizontal length direction with a strain rate of
$10^{7}$ s-1, which is a typical value used in many previous works.Jiang
(2014); Tserpes and Vatistas (2015); Xiang and Shen (2016)

The carbon-carbon interactions are described by the second generation Brenner
potential,Brenner et al. (2002) while the MoS2 interatomic interactions are
described by the Stillinger-Weber potential.Jiang et al. (2013a) The inter-
layer interactions between graphene/graphene, MoS2/MoS2 and graphene/MoS2 are
described by the Lennard-Jones potential. The distance and energy parameters
in the Lennard-Jones potential are listed in table 1.

Table 1: Lennard-Jones parameters used in the present work for graphene/MoS2,Jiang and Park (2014) graphene/graphene,Jiang and Park (2015) and MoS2/MoS2.Jiang and Park (2015) structure |  $\epsilon$ (meV) |  $\sigma$ (Å) | cut off (Å)
---|---|---|---
graphene/MoS2 | 3.95 | 3.625 | 10.0
graphene/graphene | 2.96 | 3.382 | 10.0
MoS2/MoS2 | 23.6 | 3.18 | 10.0

Our simulations are performed as follows. The system is first thermalized for
1 ns within the NPT (i.e. the particles number N, the pressure P and the
temperature T of the system are constant) ensemble by the Nose-Hoover
thermostat.Nosé (1984); Hoover (1985) Then the structures are compressed
uniaxially in the horizontal length direction, while the lateral directions
are allowed to be fully relaxed during the compression process. A value for
the thickness is required in the computation of the stress. The thickness is
not well defined for atomic thick materials like graphene and MoS2, so we have
used the same thickness of 10.0 Å for all of these three types of structures
studied in the present work.

##  III Results and discussion

Figure 2: The stress-strain relations for graphene/graphene, MoS2/MoS2, and
graphene/MoS2 heterostructures under uniaxial compression.

Figure 3: The length dependence for the critical buckling strain of
graphene/graphene, MoS2/MoS2, and the graphene/MoS2 heterostructure with
different misfit strain $\epsilon_{m}$.

Figure 2 shows the stress-strain relation for these three structures with
similar length around 47 Å. The buckling strains are 0.22%, 5.5%, and 1.1% for
the graphene/graphene bilayer, MoS2/MoS2 bilayer, and the graphene/MoS2
heterostructure, respectively. The buckling strain for the MoS2/MoS2 structure
is considerably larger than the graphene/graphene structure. This large
difference can be understood by the Euler buckling theory. According to Euler
buckling theory, the critical buckling strain of a thin plate isTimoshenko and
Woinowsky-Krieger (1987),

$\displaystyle\epsilon_{c}=-\frac{4\pi^{2}D}{EL^{2}},$ (1)

where $D$ is the bending modulus, $E$ is the in-plane stiffness, and $L$ is
the length. By linear fitting of the stress-strain relation within [0, 0.2%]
in Fig. 2, we get the stiffness of 552 GPa and 210 GPa for the
graphene/graphene and MoS2/MoS2, respectively. The bending modulus for a
single MoS2 layer is larger than a single graphene layer by a factor of
seven.Jiang et al. (2013b) As a result, we find that the critical buckling
strain for the MoS2/MoS2 is much larger than the graphene/graphene structure.

Figure 4: The misfit strain in the graphene/MoS2 heterostructures.

The length dependence for the buckling strain of these three types of
structures are shown in Fig. 3. The buckling strain for the graphene/graphene
and MoS2/MoS2 structures can be described by the Euler buckling theory. Some
deviations from the Euler buckling theory are probably due to the high strain
rate as discussed previously.Jiang (2014) For the graphene/MoS2
heterostructure, an overall information from Fig. 3 is that the critical
buckling strain for the graphene/MoS2 structure is between the values of
MoS2/MoS2 and graphene/graphene, and are closely related to the misfit strain.
To disclose the mechanism for the misfit strain effect on the buckling
phenomenon, we first derive the analytic formula for the misfit strain in the
graphene/MoS2 heterostructure. The total strain energy of a heterostructure
related to the misfit strain isAlred et al. (2015)

$\displaystyle
U=\frac{1}{2}E_{1}\epsilon_{1}^{2}+\frac{1}{2}E_{2}\epsilon_{2}^{2},$ (2)

where $\epsilon_{1}=\frac{l-l_{1}}{l_{1}}$ and
$\epsilon_{2}=\frac{l-l_{2}}{l_{2}}$ are the strain for these two atomic
layers. Here $l$ is the final length of the heterostructure structure, while
$l_{1}$ and $l_{2}$ are the original lengths for graphene and MoS2 layer,
respectively. The Young’s moduli are $E_{1}$ and $E_{2}$. The final length $l$
can be obtained by minimizing the total strain energy with respective to $l$,

$\displaystyle
l=\frac{E_{1}l_{2}+E_{2}l_{1}}{E_{1}\frac{l_{2}}{l_{1}}+E_{2}\frac{l_{1}}{l_{2}}}.$
(3)

The misfit strain is defined to be the strain in the graphene layer as follows

$\displaystyle\epsilon_{m}=\epsilon_{1}=-\frac{l_{1}-l}{l_{1}}.$ (4)

Fig. 4 shows that the misfit strain predicted by the analytic expression in
Eq. (4) agrees quite well with numerical simulation results. For comparison,
the strain in the MoS2 layer, $\epsilon_{2}=-\frac{l_{2}-l}{l_{2}}$, is also
plotted in the figure.

Figure 5: The relation between the reduced critical buckling strain and the
misfit strain of the graphene/MoS2 heterostructure. See text for the
definition of the reduced critical buckling strain.

Figure 6: Buckling modes for the graphene/MoS2 heterostructure with misfit
strain of large magnitudes (a) -1.4% and (b) 2.0%.

To demonstrate the misfit strain effect on the buckling phenomenon, we perform
a set of simulations for the graphene/MoS2 heterostructure with length varied
within the range of [40, 160] Å. The misfit strain for these structures is
within the range of [-2%, 2%]. Length and misfit strain both have effects on
the critical buckling strain. To reveal the misfit strain effect separately,
we get rid of the length effect by using the length $L$ to scale the critical
strain, i.e., $\epsilon^{L}_{c}=\epsilon_{c}L^{2}$. The quantity
$\epsilon^{L}_{c}$ should be length-independent according to the Euler
buckling theory, so this quantity solely reflects the misfit strain effect. We
further introduce a reduced critical strain
$\tilde{\epsilon}_{c}=\epsilon^{L}_{c}/\epsilon^{L}_{c0}$, in which
$\epsilon^{L}_{c0}$ is the scaled critical strain for the heterostructure
structure with minimum (almost zero) misfit strain. With this definition, the
critical buckling strain is reduced (upgraded) by the misfit strain if
$\tilde{\epsilon}_{c}<1$ ($\tilde{\epsilon}_{c}>1$). There is no effect from
the misfit strain for $\tilde{\epsilon}_{c}=1$.

Figure 5 shows the relation between the reduced critical buckling strain and
the misfit strain. The reduced critical buckling strain varies from 0.6 to 1.2
for different misfit strains; i.e., the misfit strain can affect the buckling
strain by a factor of two. There are several interesting features in the
figure. For negative misfit strain, $\epsilon_{m}\in[-2\%,0]$, graphene is
pre-compressed while MoS2 is pre-stretched according to the definition of the
misfit strain in Eq. (4). The critical buckling strain decreases with the
increase of the magnitude of the misfit strain. It should be noted that the
buckling strain for graphene is much smaller than MoS2. The buckling of the
graphene/MoS2 heterostructure thus initiates from the graphene layer. For
negative misfit strain, the graphene layer in the graphene/MoS2
heterostructure is pre-compressed by the misfit strain, so the graphene layer
can be buckled even more easily, resulting in the reduction of the critical
buckling strain. If the magnitude of the misfit strain is too large, then the
graphene layer is seriously pre-compressed while the MoS2 layer is
considerably pre-stretched. As a result, the pre-compressed graphene layer is
readily to be buckled, but the pre-stretched MoS2 is difficult to be buckled.
Eventually, the decoupling between these two atomic layers occurs as shown in
Fig. 6 (a).

If the misfit strain is in the small positive range,
$\epsilon_{m}\in[0,0.5\%]$, the graphene layer in the graphene/MoS2
heterostructure is pre-stretched and the MoS2 layer is pre-compressed by the
misfit strain. The buckling of the heterostructure still initiates from the
buckling of the graphene layer. However, the pre-stretched graphene layer
becomes more difficult to be buckled, i.e., larger compressive strain is
needed to cause the buckling of the graphene layer. As a result, the critical
buckling strain increases with the increase of the magnitude of the misfit
strain in the small positive region [0, 0.5%].

For misfit strain in a large positive range, $\epsilon_{m}\in[0.5\%,2\%]$, the
graphene layer in the heterostructure is pre-stretched while the MoS2 layer is
pre-compressed by the misfit strain. The pre-compression in the MoS2 layer is
so large that the buckling of the graphene/MoS2 heterostructure is initiated
by the MoS2 layer instead of the graphene layer. For larger misfit strain in
this positive range, the MoS2 layer is pre-compressed more seriously, so the
buckling of the MoS2 layer becomes easier. As a result, the critical buckling
strain for the graphene/MoS2 heterostructure decreases with increasing misfit
strain in [0.5%, 2%]. For positive misfit strain of large magnitude, the
decoupling can also happen as shown in Fig. 6 (b). It is because the seriously
pre-stretched graphene layer is difficult to be buckled, while the pre-
compressed MoS2 layer is easily to be buckled. The competition between these
two atomic layers leads to the decoupling phenomenon.

The above effects can be discussed qualitatively in a more general scheme. For
$\epsilon_{m}\in[-2\%,0]$ and $\epsilon_{m}\in[0,0.5\%]$, the buckling of the
graphene/MoS2 heterostructure initiates from the graphene layer, so the
buckling strain is mainly governed by the graphene layer, i.e.,
$\tilde{\epsilon}_{c}\propto(\epsilon_{c}^{\rm gra}+\epsilon_{m})$ with
$\epsilon_{c}^{\rm gra}$ as the buckling strain for graphene. For
$\epsilon_{m}\in[0.5\%,2\%]$, the buckling of the graphene/MoS2
heterostructure initiates from the MoS2 layer, and the buckling strain is
mainly governed by the MoS2 layer, i.e.,
$\tilde{\epsilon}_{c}\propto(\epsilon_{c}^{\rm MoS_{2}}-\epsilon_{m})$ with
$\epsilon_{c}^{\rm MoS_{2}}$ as the buckling strain for MoS2. Indeed, Fig. 5
shows that the numerical data can be well fitted to two linear functions in
different ranges.

Figure 7: Probability for the buckling toward the graphene side (purple
histograms) or the MoS2 side (green histograms) in the graphene/MoS2
heterostructure with different misfit strain. Top inset shows the buckling
shape towards the upper graphene side. Bottom inset shows the buckling shape
towards the lower MoS2 side.

The misfit strain also affects the buckling direction for the graphene/MoS2
heterostructure as shown in Fig. 7. More specifically, the heterostructure can
buckle toward either the upper side or the lower side. This buckling direction
depends on the misfit strain for the graphene/MoS2 heterostructure. For each
structure, we perform twenty MD simulations with different random
distributions for the initial atomic velocities, while all other simulation
parameters are the same. We count the number for the structure to buckle
upwards or downwards. For graphene/graphene and MoS2/MoS2 structures, the
buckling shape has equal probability to buckle upwards or downwards. However,
for graphene/MoS2 heterostructure of large negative misfit strain, the
buckling shape prefers to buckle toward the graphene side, because the pre-
compressed graphene layer is easily to be buckled while the pre-stretched MoS2
layer is difficult to be buckled. In contrast, for large positive misfit
strain, the buckling shape prefers to buckle toward the MoS2 side. The
probability for the buckling direction varies gradually with increasing misfit
strain. Generally, the heterostructure prefers to buckle towards the atomic
layer that is pre-compressed by the misfit strain.

##  IV Conclusion

To summarize, we have performed MD simulations to study the buckling
phenomenon of graphene/MoS2 van der Waals heterostructures under uniaxial
compression. We find that the misfit strain can affect both the critical
buckling strain and the buckling direction for the heterostructure. The misfit
strain has different effects in three strain ranges. For negative misfit
strain with graphene layer pre-compressed, the critical buckling strain can be
considerably reduced by increasing the magnitude of the misfit strain. In
small positive misfit strain region, the critical buckling strain is upgraded
by increasing misfit strain. The critical buckling strain is reduced by the
increase of the magnitude of the misfit strain for large positive misfit
strains. The effects of the misfit strain on the buckling phenomenon are
discussed based on the competition between the buckling of the graphene layer
and the MoS2 layer.

Acknowledgment The work is supported by the National Natural Science
Foundation of China (Grant Nos. 11822206 and 12072182) and Innovation Program
of the Shanghai Municipal Education Commission (Grant No.
2017-01-07-00-09-E00019).

## References

  * Novoselov et al. (2016) K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. C. Neto, Science 353, aac9439 (2016).
  * Liu et al. (2016) Y. Liu, N. O. Weiss, X. Duan, H. C. Cheng, and X. Duan, Nature Reviews Materials 1, 16042 (2016).
  * Ko et al. (2020) J. S. Ko, H. S. Dong, W. J. Lee, W. J. Chan, and S. H. Choi, Journal of Alloys and Compounds 864, 158118 (2020).
  * Ying (2016) L. Ying, Soft Matter 12, 1 (2016).
  * Xiang and Shen (2016) Y. Xiang and H. S. Shen, International Journal of Structural Stability and Dynamics 16, 1550071 (2016).
  * Tserpes and Vatistas (2015) K. I. Tserpes and I. Vatistas, Mechanics Research Communications 64, 50 (2015).
  * Pradhan and Murmu (2010) S. C. Pradhan and T. Murmu, Computational Materials Science 47, 268 (2010).
  * Shi et al. (2011) J. X. Shi, Q. Q. Ni, X. W. Lei, and T. Natsuki, Computational Materials Science 50, 3085 (2011).
  * Wen et al. (2011) H. D. Wen, G. Kai, and W. Quan, Carbon 49, 3107 (2011).
  * Zhang et al. (2012) Z. Zhang, W. H. Duan, and C. M. Wang, Nanoscale 4, 5077 (2012).
  * Zhou et al. (2015) Y. Zhou, Y. Chen, B. Liu, S. Wang, and H. Ming, Carbon 84, 263 (2015).
  * Jiang (2014) J. W. Jiang, Nanotechnology 25, 355402 (2014).
  * Bao et al. (2017) H. Bao, Y. H. Huang, Z. Yang, Y. Miao, Y. Sun, Y. Bai, K. Xu, and F. Ma, Physical Chemistry Chemical Physics 19, 25360 (2017).
  * Singh et al. (2015) S. K. Singh, M. Neek-Amal, S. Costamagna, and F. M. Peeters, Physical Review B 91, 014101 (2015).
  * Li et al. (2019) Y. Li, P. Chen, C. Zhang, J. Peng, F. Gao, and H. Liu, Computational Materials Science 162, 116 (2019).
  * Jiang and Park (2014) J. W. Jiang and H. S. Park, Applied Physics Letters 105, 033108 (2014).
  * Lin et al. (2014) Y. C. Lin, C. Y. Chang, R. K. Ghosh, J. Li, H. Zhu, R. Addou, B. Diaconescu, T. Ohta, X. Peng, N. Lu, et al., Nano Letters 14, 6936–6941 (2014).
  * Li and Zhao (2014) L. Y. Li and M. W. Zhao, Journal of Physical Chemistry C 118, 19129 (2014).
  * Ni et al. (2008) Z. H. Ni, W. Chen, X. F. Fan, J. L. Kuo, and Z. X. Shen, Physical Review B 77, 115416 (2008).
  * Chaste et al. (2018) J. Chaste, A. Missaoui, A. Saadani, D. G. Sanchez, and A. Ouerghi, ACS Applied Nano Materials 1, 6752 (2018).
  * Shivaraman et al. (2013) S. Shivaraman, J. Jobst, D. Waldmann, H. B. Weber, and M. G. Spencer, Physical Review B 87, 195425 (2013).
  * He and Jiang (2019) J. D. He and J. W. Jiang, Nanotechnology 30, 265701 (2019).
  * Plimpton (1995) S. Plimpton, Journal of Computational Physics 117, 1 (1995).
  * Stukowski (2009) A. Stukowski, Modelling and Simulation in Materials Science and Engineering 18, 015012 (2009).
  * Brenner et al. (2002) D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott, Journal of Physics:Condensed Matter 14, 783 (2002).
  * Jiang et al. (2013a) J.-W. Jiang, H. Park, and T. Rabczuk, Journal of Applied Physics 114, 064307 (2013a).
  * Jiang and Park (2015) J. W. Jiang and H. S. Park, Journal of Applied Physics 117, 124304 (2015).
  * Nosé (1984) S. Nosé, The Journal of Chemical Physics 81, 511 (1984).
  * Hoover (1985) W. G. Hoover, Physical review A 31, 1695 (1985).
  * Timoshenko and Woinowsky-Krieger (1987) S. Timoshenko and S. Woinowsky-Krieger, _Theory of Plates and Shells, 2nd ed_ (McGraw-Hill, New York, 1987).
  * Jiang et al. (2013b) J.-W. Jiang, Z. Qi, H. S. Park, and T. Rabczuk, Nanotechnology 24, 435705 (2013b).
  * Alred et al. (2015) J. M. Alred, Z. H. Zhang, Z. L. Hu, and B. I. Yakobson, Nano Research 8, 2015 (2015).