# Bi-level Trajectory Optimization on Uneven Terrains with Differentiable
Wheel-Terrain Interaction Model.

Amith Manoharan1 Aditya Sharma2 Himani Belsare2 Kaustab Pal2 K. Madhava
Krishna2 Arun Kumar Singh1 1 are with the Institute of Technology, University
of Tartu, Tartu, Estonia. 2 are with RRC, IIIT Hyderabad, India.

###### Abstract

Navigation of wheeled vehicles on uneven terrain necessitates going beyond the
2D approaches for trajectory planning. Specifically, it is essential to
incorporate the full $6dof$ variation of vehicle pose and its associated
stability cost in the planning process. To this end, most recent works aim to
learn a neural network model to predict the vehicle evolution. However, such
approaches are data-intensive and fraught with generalization issues.

In this paper, we present a purely model-based approach that just requires the
digital elevation information of the terrain. Specifically, we express the
wheel-terrain interaction and $6dof$ pose prediction as a non-linear least
squares (NLS) problem. As a result, trajectory planning can be viewed as a bi-
level optimization. The inner optimization layer predicts the pose on the
terrain along a given trajectory, while the outer layer deforms the trajectory
itself to reduce the stability and kinematic costs of the pose.

We improve the state-of-the-art in the following respects. First, we show that
our NLS based pose prediction closely matches the output from a high-fidelity
physics engine. This result coupled with the fact that we can query gradients
of the NLS solver, makes our pose predictor, a differentiable wheel-terrain
interaction model. We further leverage this differentiability to efficiently
solve the proposed bi-level trajectory optimization problem. Finally, we
perform extensive experiments, and comparison with a baseline to showcase the
effectiveness of our approach in obtaining smooth, stable trajectories.

##  I Introduction

Applications like forestry, construction, and search and rescue require
wheeled vehicles to navigate over uneven terrains. Thus, for automation in
these use cases, we need a trajectory planner that can account for the vehicle
tip-over stability over uneven terrains. This in turn requires explicitly
predicting the vehicle’s $6dof$ pose on the terrain and the kinematic and
stability costs associated with those poses.

There are broadly two ways to predict the pose of a wheeled vehicle on an
uneven terrain. The first approach is to use high-fidelity physics simulators
that can model complex wheel-terrain interaction [1], often in the form of
complex differential equations [2]. However, these simulators are often too
slow to be used during online planning. Alternately, wheel-terrain interaction
and pose prediction can be learned from recorded vehicle motions on different
terrains [3]. However, purely data-driven approaches are data-hungry and often
struggle to generalize to novel scenarios.

Our main motivation in this paper is to make the physics/model-based approach
more computationally tractable by finding the appropriate simplified
abstraction for the wheel-terrain interaction. Moreover, we aim to make the
interaction model differentiable and parallelizable to ensure online planning.
We achieve these features through some core algorithmic innovations, which are
summarized below, along with their benefits.

Figure 1: Husky navigating an uneven terrain. The black points represent the
entire terrain, whereas the coloured points represent the terrain’s local
patch. All points within a given radius from the Husky are used to model the
terrain, which is used to compute a stable trajectory to the goal. As seen
above, the trajectory (shown in pink) successfully avoids the ditch and goes
around it to reach the goal.

Algorithmic Contribution: We jointly view the wheeled vehicle and the
underneath terrain through the lens of closed-loop kinematics, typically used
for modeling parallel manipulators. This in turn leads to a set of coupled
non-linear equations that dictate how the vehicle’s $6dof$ pose depends on the
terrain surface geometry. The pose prediction can then be obtained by
minimizing the residual of the non-linear equations, which, in our case, is
done through a least-squares optimization problem. Importantly, our set-up
allows for differentiating through the least squares solver using implicit
function theorem to ascertain how the pose of the vehicle will vary if its
position on the terrain is perturbed. Thus, in essence, our non-linear least
squares (NLS) approach provides a differentiable wheel-terrain interaction
model, albeit at the kinematic level.

The differentiability feature also opens up efficient ways of performing
trajectory planning on uneven terrain. Specifically, for the first time, we
show that the planning problem can be viewed as a bi-level optimization. Its
inner layer (NLS) predicts pose along a given trajectory while the outer layer
perturbs the trajectory itself to lower the kinematic and stability costs. The
implicit gradients from the inner layer can be used within any momentum-based
gradient descent approach to efficiently solve the bi-level problem.

State-of-the-Art Performance: We show that the prediction from our NLS-based
approach closely matches the output from a high-fidelity physics engine. As a
result, the former provides a computationally cheap alternative for learning-
based approaches that require tele-operating vehicles on potentially dangerous
terrains. We show that by modeling the stability costs through force-angle
measure [4], our trajectory optimizer produces trajectories that successfully
avoid risky areas such as ditches and valleys. Finally, we show that our
gradient-based approach is competitive with more computationally demanding
sampling-based optimization, such as the Cross Entropy Method [5] that does
not require access to implicit gradients.

##  II Related Works

#### Learning Based Wheel-Terrain Interaction

An unsupervised method for the classification of a ground rover’s slip events
based on proprioceptive signals is proposed in [6]. The method is able to
automatically discover and track various degrees of slip on the given terrain.
A meta-learning-based approach is proposed in [7] to predict future dynamics
by leveraging rover-terrain interaction data. In [8], an approach for learning
and predicting the motion resistance encountered by a vehicle while traversing
is studied using the information obtained from a stereovision device. A
learning-based approach for predicting the pitch and roll angle from a depth
map was proposed in [9]. Authors in [10] show that purely data-driven
approaches struggle to accurately predict the vehicle pose on the terrain.
Thus, they embed physics-based priors into the learning pipeline. Our current
approach is closely related [10], as our NLS-based pose predictor can also be
potentially used as a physics-based prior.

#### Motion Planning on Uneven Terrains

Model-based planning on uneven terrains requires a pose predictor module. This
can be either a learned model as used in [3], or that obtained from first
principles of physics and differential equations [11]. A hybrid trajectory
optimization method for generating stable and smooth trajectories for
articulated tracked vehicles is proposed in [12]. A vehicle-terrain contact
model is developed to divide the vehicle’s motion into hybrid modes of driving
and traversing, which is converted into a nonlinear programming problem to be
solved in a moving-horizon manner. A terrain pose mapping to describe the
impact of terrain on the vehicle is proposed in [13], based on which a
trajectory optimization framework for car-like vehicles on uneven terrain is
built. Our proposed method differs from these cited works in how we leverage
implicit gradients from a first principle-based wheel-terrain interaction
model for trajectory optimization.

#### End-to-End Learning Based Approaches

A learning-based end-to-end navigation for planetary rovers is given in [14].
A deep reinforcement learning-based navigation method is proposed to
autonomously guide the rover towards goals through paths with low wheel slip
ratios. An end-to-end network architecture is proposed, in which the visual
perception and the wheel-terrain interaction are combined to learn the terrain
properties implicitly. In [15], a sim-to-real pipeline for a vehicle to learn
how to navigate on rough terrains is presented. A deep reinforcement learning
architecture is used to learn a navigation policy from data obtained from
simulated environments. An approach that employs a fully-trained deep
reinforcement learning network that uses elevation maps of the environment,
vehicle pose, and goal as inputs to compute an attention mask of the
environment is proposed in [16]. The attention mask is then used to identify
reduced stability regions in the terrain.

#### Connection to Author’s Prior Work:

The proposed NLS-based pose predictor builds upon and improves our prior work
[17]. Herein, small angle approximation and local planar assumption on the
underlying terrain were used to obtain closed-form solutions for pose
prediction. However, such an approach is not appropriate for highly uneven
terrains.

##  III Preliminaries

(a)

(b)

Figure 2: (a) A four-wheeled vehicle with the geometry vectors describing the
holonomic constraints. (b) Definition of the vectors associated with tip-over
stability.

#### Symbols and Notation

Scalars are denoted by normal lowercase letters, whereas bold font variants
are used for vectors. Matrices are represented by uppercase bold font letters.
The superscript $T$ is used to define the transpose of a matrix or a vector.

###  III-A Primer on Chain Rule of Differentiation

This section covers some basic notations from multi-variable calculus that
will be crucial for developing the implicit differentiation rules later in the
paper. Most of the developments follow from the work [18].

Consider a function $\psi(\lambda,\phi,\mu)$ where $\phi$ and $\mu$ are
themselves functions of $\lambda$. In such a scenario, the chain rule of
differentiation takes the following form.

$\frac{\mathrm{d}}{\mathrm{d}\lambda}\psi=\frac{\partial\psi}{\partial\lambda}\frac{\mathrm{d}\lambda}{\mathrm{d}\lambda}+\frac{\partial\psi}{\partial\phi}\frac{\mathrm{d}\phi}{\mathrm{d}\lambda}+\frac{\partial\psi}{\partial\mu}\frac{\mathrm{d}\mu}{\mathrm{d}\lambda}.$
(1)

For a scalar function over multiple variables,
$\psi(\lambda_{1},\lambda_{2},\dots\lambda_{n}):\mathbb{R}^{n}\rightarrow\mathbb{R}$,
the derivative vector can be written as

$\mathrm{D}\psi=\begin{bmatrix}\frac{\partial\psi}{\partial\lambda_{1}},\frac{\partial\psi}{\partial\lambda_{2}},\dots,\frac{\partial\psi}{\partial\lambda_{n}}\end{bmatrix}\in\mathbb{R}^{1\times
n}.$ (2)

Analogously, for vector-valued functions
$\boldsymbol{\psi}:\mathbb{R}\rightarrow\mathbb{R}^{m}$, defined over a scalar
argument, we can write

$\mathrm{D}\psi=\begin{bmatrix}\frac{\mathrm{d}\psi_{1}}{\mathrm{d}\lambda},\dots,\frac{\mathrm{d}\psi_{m}}{\mathrm{d}\lambda}\end{bmatrix}\in\mathbb{R}^{m\times
1}.$ (3)

The general definition of $\mathrm{D}\boldsymbol{\psi}$ of
$\boldsymbol{\psi}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ can be written as
an $m\times n$ matrix with entries

$\left(\mathrm{D}\boldsymbol{\psi}(\lambda)\right)_{ij}=\frac{\partial\boldsymbol{\psi}_{i}}{\partial\lambda_{j}}(\lambda).$
(4)

Now, the chain rule for $h(\lambda)=\zeta(\psi(\lambda))$ can be written as

$\mathrm{D}h(\lambda)=\mathrm{D}\zeta(\psi(\lambda))\mathrm{D}\psi(\lambda).$
(5)

For the partial derivatives, the variable over which the derivative is
computed is written as a subscript. For example,
$\mathrm{D}_{\lambda}\psi(\lambda,\phi)$ denotes partial derivative with
respect to $\lambda$. Similarly, $\mathrm{D}_{\lambda\phi}^{2}\psi$ means
$\mathrm{D}_{\lambda}(\mathrm{D}_{\phi}\psi)^{T}$.

Finally, when subscripts are not given for multi-variate functions, the
operator $\mathrm{D}$ means the total derivative with respect to the
independent variables. So, with $\lambda$ as the independent variable, the
vector form of Equation (1) is represented by

$\mathrm{D}\psi=\mathrm{D}_{\lambda}\psi+\mathrm{D}_{\phi}\psi\mathrm{D}\phi+\mathrm{D}_{\mu}\psi\mathrm{D}\mu.$
(6)

###  III-B Trajectory Parametrization

We assume that the vehicle follows a bi-cycle model. Moreover, we can leverage
its differential flatness property [19] to plan directly in the positional
space ($x_{k},y_{k}$). The control inputs, for example, linear velocity and
curvature, can then be obtained as a function of positional variables and
their higher-order derivatives. We parametrize the position-level trajectory
of the vehicle in terms of polynomials in the following form:

$\displaystyle\begin{bmatrix}x_{1},\dots,x_{n}\end{bmatrix}=\mathbf{W}\mathbf{c}_{x},\begin{bmatrix}y_{1},\dots,y_{n}\end{bmatrix}=\mathbf{W}\mathbf{c}_{y},$
(7)

where, $\mathbf{W}$ is a matrix formed with time-dependent polynomial basis
functions and ($\mathbf{c}_{x},\mathbf{c}_{y}$) are the coefficients of the
polynomial. We can also express the derivatives in terms of
$\dot{\mathbf{W}},\ddot{\mathbf{W}}$.

###  III-C Tip-over stability criteria

A vehicle’s stability at a given point on the terrain can be calculated by the
force angle measure or tip-over stability criteria proposed in [20, 4]. Fig.
2(b) shows the polygon formed with the ground contact points of a four-wheel
vehicle. The edges of the polygon given by the vectors $\mathbf{e}_{i}$
represent the tip-over axes, i.e., the axes about which the vehicle could
possibly tip-over. These axes can be computed in the following way:

$\displaystyle\mathbf{e}_{i}$
$\displaystyle=\mathbf{p}_{oc,i+1}-\mathbf{p}_{oc,i},\forall i=1,2,3,$ (8)
$\displaystyle\mathbf{e}_{4}$
$\displaystyle=\mathbf{p}_{oc,1}-\mathbf{p}_{oc,4}.$

We define the unit vectors
$\mathbf{\hat{e}}_{i}=\frac{\mathbf{e}_{i}}{\lVert\mathbf{e}_{i}\rVert}$,
$\boldsymbol{\hat{\pi}}_{i}=\frac{\boldsymbol{\pi}_{i}}{\lVert\boldsymbol{\pi}_{i}\rVert}$,
and
$\boldsymbol{\hat{\nu}}_{i}=\frac{\boldsymbol{\nu}_{i}}{\lVert\boldsymbol{\nu}_{i}\rVert}$,
where $\boldsymbol{\pi}_{i}$ are the tip-over axis normals which intersect the
vehicle’s center of mass (COM) and $\boldsymbol{\nu}_{i}$ are the components
of the force acting on the vehicle which contribute to the tip-over. These
vectors can be calculated as:

$\displaystyle\boldsymbol{\pi}_{i}$
$\displaystyle=(\mathrm{I}-\mathbf{\hat{e}}_{i}\mathbf{\hat{e}}_{i}^{T})(\mathbf{p}_{oc,i+1}-\mathbf{p}_{og}),$
(9) $\displaystyle\boldsymbol{\nu}_{i}$
$\displaystyle=(\mathrm{I}-\mathbf{\hat{e}}_{i}\mathbf{\hat{e}}_{i}^{T})\boldsymbol{\nu},$
(10)

where $\mathrm{I}$ is the identity matrix and $\boldsymbol{\nu}$ is the total
force acting on the vehicle. The angles between the forces and the tip-over
axis normals can be computed as:

$\theta_{i}=\sigma_{i}\arccos{(\boldsymbol{\hat{\nu}}_{i}\cdot\boldsymbol{\hat{\pi}}_{i})},$
(11)

where,

$\sigma_{i}=\begin{cases}+1,&(\boldsymbol{\hat{\pi}}_{i}\times\boldsymbol{\hat{\nu}}_{i})\cdot\mathbf{\hat{e}}_{i}<0\\\
-1,&\text{otherwise},\end{cases}$ (12)

from which the tip-over stability criteria can be defined as:

$\theta_{i}>0\implies\min\theta_{i}>0,\forall i=1,2,3,4.$ (13)

We define the costs associated with the tip-over stability criteria as:

$\displaystyle\theta_{c,i}$ $\displaystyle=\max(0,-\theta_{i}+\epsilon),$ (14)
$\displaystyle\theta_{\Delta c,i}$
$\displaystyle=(\theta_{i}-\theta_{i+1})^{2},$ (15)

where $\epsilon$ is a small number. The total stability cost for the vehicle
is defined as:

$c_{s}=\sum_{i}(\theta_{c,i}+w_{\theta}\,\theta_{\Delta c,i}),$ (16)

where $w_{\theta}$ is a scaling factor.

##  IV Main Algorithmic Results

Our trajectory planning pipeline has the following components. First, we fit a
functional form to the terrain digital elevation data, which is then used to
formulate a differentiable wheel-terrain interaction model. Subsequently, we
formulate a trajectory planning problem that leverages the differentiability
of the wheel-terrain interaction.

###  IV-A Functional Form for the Terrain Model

Imagine that we are given the digital elevation data, which describes the
height $z^{j}$ of the terrain at any given point $(x^{j},y^{j})$. For example,
these can be obtained online from the point clouds obtained with on-board
LiDAR. In this section, our aim is to fit a functional form to this data. That
is, we want to obtain analytical relationships of the form:

$z^{j}=f(x^{j},y^{j}).$ (17)

Our approach is built on approximating $f$ in terms of Fourier basis
functions. That is,

$f(x^{j},y^{j})=\sum_{n=1}^{N}a_{n}\cos{(\omega_{1,n}x^{j}+\omega_{2,n}y^{j})}+b_{n}\sin{(\omega_{3,n}x^{j}+\omega_{4,n}y^{j})},$
(18)

where $\omega_{1,n},\omega_{2,n},\omega_{3,n},\omega_{4,n}$ are the
frequencies of the Fourier basis functions, $N$ is the number of frequencies,
and $a_{n},b_{n}$ are the weights associated with each functions. These are
obtained by solving the following regression problem

$\displaystyle\sum_{j=1}^{M}\|f(x^{j},y^{j})-z^{j}\|_{2}^{2},$ (19)

where $M$ is the total number of points in the digital elevation data or
point-cloud. Fig.3 shows the fit obtained over circular terrain patches of
radius $7m$.

Practical Considerations: The number of frequencies needed to get a good fit
depends on the area of the terrain. Our implementation, which leverages warm-
started gradient descent on GPUs, takes on an average of $0.5s$ to fit the
Fourier approximation to a circular patch of radius $7m$. During practical
implementation, we can adopt a receding horizon approach, where the vehicle
can fit the terrain model over a certain distance and then plan a stable
trajectory over it. Such an approach would ensure that the terrain fit is
called sparingly.

(a)

(b)

(c)

(d)

Figure 3: Red points represent the ground truth terrain patch, and green
points represent the terrain patch with predicted z coordinates.

###  IV-B Differentiable Wheel-Terrain Interaction

This section presents the first core result of the paper. Let the yaw plane
states of the vehicle at time $k$ can be represented by the vector
$\mathbf{x}_{k}=\begin{bmatrix}x_{k}&y_{k}&\alpha_{k}\end{bmatrix}^{T}$, where
$(x_{k},y_{k})$ is the position, and $\alpha_{k}$ is the heading angle of the
vehicle. For wheeled vehicles with no-active suspension mechanism, the
evolution of a vehicle’s yaw plane configuration, i.e., $\mathbf{x}_{k}$, can
be directly controlled. The rest of the pose variables, i.e., $z_{k}$
coordinate, roll $\beta_{k}$, and pitch $\gamma_{k}$, will be a function of
the yaw plane configuration and the terrain geometry. Mathematically, this
dependency can be represented in the following manner.

$\displaystyle z_{k}=s_{1}(x_{k},y_{k},\alpha_{k}),$
$\displaystyle\beta_{k}=s_{2}(x_{k},y_{k},\alpha_{k}),$ (20)
$\displaystyle\gamma_{k}=s_{3}(x_{k},y_{k},\alpha_{k}).$

In this section, we are interested in obtaining the functions $s_{i}$ that
decide the wheel-terrain interaction. To this end, we view the wheeled vehicle
on uneven terrain as a parallel manipulator [21] and derive the so-called loop
closure equations. To this end, refer to Fig. 2(a) and imagine that we have a
global reference frame $\\{O\\}$ and a reference frame $\\{G\\}$ that moves
along with the vehicle and is attached to its center. The reference frame
$\\{L\\}$ is also similar, but its orientation is the same as the vehicle.
With these notations in place, the vector describing the wheel-ground contact
point $\mathbf{p}_{oc,i}$ can be directly obtained from the terrain equation.
Alternately, the contact points can also be reached by first following through
$\mathbf{p}_{og}$ that identifies the center of the vehicle and then moving
along $\mathbf{p}_{gc,i}$. This insight can be encoded into the following
equation, which resembles the loop closure equations obtained for parallel
manipulators [21] [17].

$\mathbf{p}_{og}+\mathbf{p}_{gc,i}=\mathbf{p}_{oc,i},$ (21)

$\displaystyle\mathbf{p}_{gc,i}$
$\displaystyle=\mathbf{R}\begin{bmatrix}\delta_{i}h&r_{i}w&-(l_{i})\end{bmatrix},\forall
i=1,2,3,4,$ (22) $\displaystyle\delta_{i}$
$\displaystyle=\begin{cases}1,&i=1,4,\\\ -1,&i=2,3,\end{cases}$ $\displaystyle
r_{i}$ $\displaystyle=\frac{2.5-i}{|2.5-i|},$ $\displaystyle\mathbf{p}_{og}$
$\displaystyle=\begin{bmatrix}x_{k}&y_{k}&z_{k}\end{bmatrix},$ (23)
$\displaystyle\mathbf{p}_{oc,i}$
$\displaystyle=\begin{bmatrix}x_{{c,i}_{k}}&y_{{c,i}_{k}}&z_{{c,i}_{k}}\end{bmatrix}.$
(24)

The term $\mathbf{p}_{gc,i}$ is expressed as the vector describing the wheel
contact point in the vehicle local frame $\\{L\\}$ multiplied with the
rotation matrix between the $\\{O\\}$ and $\\{L\\}$ (refer Fig. 2(a)). The
variables $r_{i}$ and $\delta_{i}$ have been incorporated to ensure the proper
sign of $w$ and $h$ corresponding to each vertex of the chassis. $l_{i}$ are
the equivalent leg lengths, including the wheels’ radius. The constants $h$
and $w$ are the half-width and half-breadth of the chassis.

We also note that $\mathbf{p}_{oc,i}$ satisfy the terrain equation (17) and
thus, it is possible to write:

$z_{{c,i}_{k}}=f(x_{{c,i}_{k}},y_{{c,i}_{k}}).$ (25)

We can expand (21) using the Euler-angle parametrization for the rotation
matrix $\mathbf{R}$. When the resulting expressions are stacked alongside
(25), we obtain a set of coupled non-linear equations in the following form:

$g_{j}(\mathbf{x}_{k},\mathbf{u}_{k})=0,\quad j=1\dots 16,$ (26)

where
$\mathbf{u}_{k}(\mathbf{x}_{k})=\begin{bmatrix}z_{k}&\beta_{k}&\gamma_{k}&x_{{c,i}_{k}}&y_{{c,i}_{k}}&z_{{c,i}_{k}}\end{bmatrix}^{T}$.
The pose prediction can now be formulated as a non-linear least squares (NLS)
problem:

$\mathbf{u}_{k}^{*}(\mathbf{x}_{k})=\operatorname*{arg\,min}_{\mathbf{u}_{k}}\sum_{j}g_{j}(\mathbf{x}_{k},\mathbf{u}_{k})^{2}.$
(27)

###### Remark 1

The residual associated with the $\mathbf{u}_{k}^{*}(\mathbf{x}_{k})$ can be
used to ascertain whether it is possible for the wheeled vehicle to ensure
contact on all four wheels on a particular patch of terrain.

Implicit Differentiation For an efficient solution of the trajectory
optimization problem introduced later, it is important to compute the Jacobian
of $\mathbf{u}_{k}^{*}(\mathbf{x}_{k})$ with respect to its input parameter
$\mathbf{x}_{k}$. However, the relationship between $\mathbf{u}_{k}^{*}$ and
$\mathbf{x}_{k}$ does not have an analytical form and thus requires tools from
implicit differentiation. The derivation is based on Dini’s implicit function
theorem [22] applied to the first-order optimality condition. We define the
following proposition that has been adapted from [18] for our pose prediction
NLS.

###### Proposition 1

Consider the NLS problem (27). We can define
$\mathbf{H}=\mathrm{D}_{\mathbf{u}_{k}\mathbf{u}_{k}}^{2}\mathbf{g}(\mathbf{x}_{k},\mathbf{u}_{k}(\mathbf{x}_{k}))\in\mathbb{R}^{m\times
m}$, where $\mathbf{g(.)}$ is obtained by stacking $g_{j}(.)$ and
$\mathbf{B}=\mathrm{D}_{\mathbf{x_{k}u_{k}}}^{2}\mathbf{g}(\mathbf{x}_{k},\mathbf{u}_{k}(\mathbf{x}_{k}))\in\mathbb{R}^{m\times
n}$. Then, the Jacobian of the optimal pose is

$\mathrm{D}\mathbf{u}^{*}_{k}(\mathbf{x}_{k})=-\mathbf{H}^{-1}\mathbf{B},$

with the assumption that $\mathbf{H}$ is non-singular.

###### Proof:

The first-order optimality condition of our NLS is given by
$\mathrm{D}_{\mathbf{u}_{k}}\mathbf{g}(\mathbf{x}_{k},\mathbf{u}_{k})=\mathbf{0}_{1\times
m}$. Differentiating both sides of this optimality condition with respect to
$\mathbf{x}_{k}$ provides us the following equation

$\displaystyle\mathbf{0}_{m\times n}$
$\displaystyle=\mathrm{D}(\mathrm{D}_{\mathbf{u}_{k}}\mathbf{g}(\mathbf{x}_{k},\mathbf{u}_{k}))^{T}$
$\displaystyle=\mathrm{D}_{\mathbf{x}_{k}\mathbf{u}_{k}}^{2}\mathbf{g}(\mathbf{x}_{k},\mathbf{u}_{k})+\mathrm{D}_{\mathbf{u}_{k}\mathbf{u}_{k}}^{2}\mathbf{g}(\mathbf{x}_{k},\mathbf{u}_{k})\mathrm{D}\mathbf{u}_{k}(\mathbf{x}_{k}),$

which can be rearranged to

$\mathrm{D}\mathbf{u}_{k}(\mathbf{x}_{k})=-(\mathrm{D}_{\mathbf{u}_{k}\mathbf{u}_{k}}^{2}\mathbf{g}(\mathbf{x}_{k},\mathbf{u}_{k}))^{-1}\mathrm{D}_{\mathbf{x}_{k}\mathbf{u}_{k}}^{2}\mathbf{g}(\mathbf{x}_{k},\mathbf{u}_{k}),$
(28)

when
$\mathrm{D}_{\mathbf{u}_{k}\mathbf{u}_{k}}^{2}\mathbf{g}(\mathbf{x}_{k},\mathbf{u}_{k})$
is non-singular. ∎

(a)

(b)

(c)

(d)

(e)

Figure 4: (a) Trajectory obtained from manually driving Husky on a synthetic
terrain in Gazebo. (b) $z$ Ground truth and predicted values. (c) $\beta$ and
$\gamma$ Ground truth vs predicted. (d) Error statistics for $z$. (e) Error
statistics for $\beta$ and $\gamma$.

###  IV-C Bi-Level Optimization Based Trajectory Planning

We assume a differentially flat non-holonomic vehicle. Thus, its heading,
velocity, and curvature can be uniquely defined by the position-level
trajectory. For example, the heading can be obtained as
$\arctan(\dot{y}_{k},\dot{x}_{k})$. With this insight, we propose the
following bi-level optimization problem

$\displaystyle\min\sum_{k}$
$\displaystyle\left(c_{r}(\ddot{x}_{k},\ddot{y}_{k},\dot{x}_{k},\dot{y}_{k})+c_{s}(\mathbf{u}^{*}_{k}(\textbf{x}_{k}))\right),$
(29a)
$\displaystyle(x_{0},y_{0},\dot{x}_{0},\dot{y}_{0},\ddot{x}_{0},\ddot{y}_{0})=\mathbf{b}_{0},$
(29b)
$\displaystyle(x_{n},y_{n},\dot{x}_{n},\dot{y}_{n},\ddot{x}_{n},\ddot{y}_{n})=\mathbf{b}_{n},$
(29c) $\displaystyle\underline{x}\leq x_{k}\leq\overline{x},$ (29d)
$\displaystyle\underline{y}\leq y_{k}\leq\overline{y},$ (29e)
$\displaystyle\mathbf{u}_{k}^{*}(\mathbf{x}_{k})=\operatorname*{arg\,min}_{\mathbf{u}_{k}}\sum_{j}g_{j}(\mathbf{x}_{k},\mathbf{u}_{k})^{2}.$
(29f) $\displaystyle
c_{a}(\ddot{x}_{k},\ddot{y}_{k})=\ddot{x}_{k}^{2}+\ddot{y}_{k}^{2},$ (29g)
$\displaystyle
c_{c}(\ddot{x}_{k},\ddot{y}_{k},\dot{x}_{k},\dot{y}_{k})=\ddot{y}\dot{x}-\ddot{x}\dot{y}/(\dot{x}^{2}+\dot{y}^{2}+\epsilon)^{\frac{3}{2}},$
(29h) $\displaystyle c_{r}(.)=c_{a}(.)+c_{c}(.).$ (29i)

The first term ($c_{r}(.)$) in the cost function (29a) is the kinematics cost
that ensures smoothness in the planned trajectory by penalizing high
accelerations and sharp turns in the trajectory. The last term ($c_{s}(.)$) is
the tip-over stability cost defined in (16). Equality constraints (29b) and
(29c) ensure that the planned trajectory satisfies the initial and final
boundary conditions. Consequently, $\mathbf{b}_{0},\mathbf{b}_{n}$ are the
stacked vectors of initial and final positions, velocities, and accelerations.
The constants $\underline{x},\underline{y},\overline{x},\overline{y}$ are the
minimum and maximum bounds for $x_{k},y_{k}$. The position level is critical
for restricting the trajectory to lie within the patch that is currently
observable by the vehicle’s sensors.

Using the trajectory parametrization given in Sec III-B, the Bi-level
optimization problem can be written in a compact form as given below.

$\displaystyle\min_{\boldsymbol{\xi}}c_{r}(\boldsymbol{\xi})+c_{s}(\mathbf{u}^{*}(\boldsymbol{\xi})),$
(30a)
$\displaystyle\mathbf{u}(\boldsymbol{\xi})=\min_{\mathbf{u}}\left\lVert\mathbf{g}(\boldsymbol{\xi},\mathbf{u})\right\rVert_{2}^{2},$
(30b)
$\displaystyle\textbf{A}_{eq}\boldsymbol{\xi}=\textbf{b}_{eq},\quad\textbf{A}\boldsymbol{\xi}\leq\textbf{b}.$
(30c)

where $\mathbf{u}^{*}$ is formed by stacking $\mathbf{u}_{k}^{*}$ at different
instants of time. The matrices $\mathbf{A}_{eq},\mathbf{A}$ and the vectors
$\mathbf{b}_{eq},\mathbf{b}$ comes from (29b)-(29e).

###  IV-D Projected gradient descent

The optimal value for the parameter $\boldsymbol{\xi}_{k}$ is obtained through
projected gradient descent. The gradient descent update rule is given as

${}^{t+1}\overline{\boldsymbol{\xi}}=\,^{t}\boldsymbol{\xi}-\eta\left(\nabla_{\boldsymbol{\xi}}c_{r}({{}^{t}}\boldsymbol{\xi})+\mathrm{D}c_{s}\left(\mathbf{u}^{*}\left({{}^{t}}\boldsymbol{\xi}\right)\right)\right),$
(31)

where $\eta$ is the learning rate. The parameter $\boldsymbol{\xi}$ is
projected onto the constraint set using the following Quadratic Programming
problem at each iteration $t$ of the gradient descent.

${}^{t+1}\boldsymbol{\xi}=$
$\displaystyle\arg\min_{\boldsymbol{\xi}}\frac{1}{2}\left\lVert{}^{t+1}\overline{\boldsymbol{\xi}}-\boldsymbol{\xi}\right\rVert_{2}^{2},$
(32) $\displaystyle\textbf{A}_{eq}\,^{t+1}\boldsymbol{\xi}=\textbf{b}_{eq},$
(33) $\displaystyle\textbf{A}\,^{t+1}\boldsymbol{\xi}\leq\textbf{b}.$ (34)

The term
$\mathrm{D}c_{s}\left(\mathbf{u}^{*}\left({{}^{t}}\boldsymbol{\xi}\right)\right)$
will be computed through implicit differentiation discussed in the previous
section. However, since the implicit gradient is computed for a scalar cost
$c_{s}$, we never actually have to explicitly form the Jacobian of
$\textbf{u}^{*}$ with respect to $\boldsymbol{\xi}$. Rather, we just need the
Jacobian-vector product that can be efficiently computed through existing
autodifferentiaton libraries like JAX [23].

(a) Example 1

(b) Example 2

(c) Example 3

(d) Example 4

Figure 5: Trajectories with and without stability cost. Black colour
represents trajectory with stability cost and pink colour without the
stability cost.

##  V Results

In this section, we validate our proposed wheel-terrain interaction and pose
predictor along with the bi-level optimizer through extensive simulations. We
mainly focus on answering the following research questions with the help of
qualitative and quantitative (statistical) results.

  * •

Q1: How well does our NLS based pose predictor fare vis-a-vis high-fidelity
physics simulator Gazebo [24]?

  * •

Q2: How does the stability cost affect the trajectory taken by the vehicle?

  * •

Q3: How does our gradient descent approach compare to the state-of-the-art
Cross Entropy Method (CEM) [5] sampling based planner?

(a) Example 1

(b) Example 2

(c) Example 3

(d) Example 4

Figure 6: Roll and pitch angle variations with and without the stability cost.

###  V-A Implementation Details

The trajectory planner comprising of the bi-level optimization (30a)-(30c) and
projection (33)-(34) is implemented in Python using the JAX [23] library for
GPU-accelerated computing. The matrix W in (7) is constructed from a $10^{th}$
order polynomial. Our simulation pipeline consists of open-loop analyses
carried out with stand-alone Python code, and the closed-loop simulations were
performed on the Gazebo simulator [24]. The terrains were created using
Blender [25] and were converted to the required formats for the Python and
Gazebo simulations. All simulations were conducted on a computer with Intel i9
CPU and Nvidea RTX 3080 GPU. For modeling the terrain at each instant of time,
we take a patch of radius $7m$ centered at the vehicle position. The number of
Fourier frequencies $N$ in the terrain model was taken as $100$ to balance the
trade-off between computation time and accuracy.

###  V-B Validation using the Gazebo simulator

In this section, we analyze the accuracy of our pose predictor by comparing
the ground truth data obtained from Gazebo with the values obtained by solving
the NLS problem (27) with the help of (17). To this end, we drive the Husky
robot in Gazebo on synthetic terrain models and collect the odometry data, as
shown in Fig. 4(a). We then compare the ground truth values $z_{gt}$,
$\beta_{gt}$, and $\gamma_{gt}$ from Gazebo with the predicted $z_{pred}$,
$\beta_{pred}$, and $\gamma_{pred}$ from our pose optimizer. Several runs of
the robot on different terrains were performed, and some qualitative results
are shown in Figures 4(b) and 4(c). We can see that the predicted pose values
closely match the ground truth values. The statistical results presented in
Figures 4(d) and 4(e) show the minimum, maximum, and median error obtained
across several runs. It can be seen that our NLS pose optimizer is highly
accurate in predicting the pose of the robot on different terrains, and the
error values are minimal.

###  V-C Validation of the Safe Planning on Uneven Terrain

In this subsection, we evaluate the performance of our trajectory optimizer by
considering the tip-over stability metric as an indicator of safe
trajectories. Fig 5 shows four example trajectories generated by the planner
with and without the stability cost $c_{s}(.)$ defined in (16). From the
figures, we can see that the trajectory computed by the planner with stability
cost closely follows the terrain gradients, while the one without the
stability cost cuts across hills and valleys in the terrain. The former
behavior results in a safer and more stable trajectory toward the goal point,
whereas the latter is prone to tipping along the way. The variations in the
roll angle $\beta$ and the pitch angle $\gamma$ for the trajectories are
plotted in Fig. 6. It can be seen that the magnitude of the angles is higher
for the trajectory without the stability cost. Also, there are sharp
variations in the angles, which means that the vehicle is traversing through a
very uneven terrain. Sudden changes in the vehicle’s orientation may lead to
the loss of stability. With the stability cost in effect, the magnitude and
variations in the angles become lower, and the vehicle becomes more stable.

Next, we present a statistical result in Fig. 7(a). Here, we also highlight
the significance of the weighting term $w_{\theta}$ in the stability cost 16.
Several trajectories were simulated on different terrains, and the average
worst-case tip-over angles were recorded. It can be seen that for the
trajectories without stability cost, the value is lower. In the case of tip-
over angle, higher positive values are preferable, as shown in (13). The
values with stability cost are far higher than those without stability cost,
which confirms that the trajectories taken by the vehicle with stability cost
have a higher safety margin. It can also be seen that as the weight
$w_{\theta}$ increased from 0.05 to 0.2, the values of the angles increased,
which in turn increased the tip-over stability margin.

###  V-D Comparison with CEM sampling-based planner

This section investigates the importance of differentiability of the wheel-
terrain interaction/pose predictor model. If the implicit gradients are not
available, then one can use a black-box sampling-based optimizer such as CEM.
However, it will require an order of magnitude larger calls of the NLS pose
predictor. In fact, it may not even be possible to solve our bi-level
optimizer through a sampling-based approach on resource-constrained hardware.
Nevertheless, in this section, our aim is to show that our gradient-based
approach is competitive with a more computationally demanding approach. The
results are summarized in Fig 7(b), which shows the statistical results of the
average worst-case tip-over angle obtained with our gradient-based approach
and CEM. Two different batch sizes (100, 20) were used for CEM. The results
show that our gradient descent-based planner performs as well as a CEM planner
with batch size 100 and is slightly better than a CEM with batch size 20. The
average computation time for our approach was $0.42s$, whereas the CEM with
batch sizes 100 and 20 took $2.12s$ and $1.03s$, respectively. It is also to
be noted that the CEM planners are memory-hungry; thus, implementation on
onboard computers might be difficult, if at all possible.

(a) a: with $c_{s}$, $w_{\theta}=0.2$, b: with $c_{s}$, $w_{\theta}=0.05$, c:
without $c_{s}$

(b) gradient descent compared with CEM with batch sizes 100, 20.

Figure 7: (a) Comparison of the effect of the stability cost, $c_{s}$. (b)
Comparison between gradient descent and CEM.

##  VI conclusions

We presented, for the first-time, a differentiable wheel-terrain
interaction/pose predictor model that is derived purely from the first
principles and yet generalizes to arbitrary terrains. The fidelity of our
model was shown to be close to that of a state-of-the-art physics engine. This
opens up exciting new possibilities. For example, our pose predictor can be
used as a differentiable and parallelizable world model within trajectory
optimization or even reinforcement learning pipelines. We presented one such
approach based on bi-level optimization, wherein we leveraged implicit
gradients from the NLS solver to efficiently compute the optimal trajectory.
The proposed optimizer can accommodate arbitrary stability metrics that depend
on the vehicle’s pose on the uneven terrain. As an example, we used the force-
angle measure and demonstrated stable trajectory planning on uneven terrains.
Our future endeavors are focused towards extending the NLS-based approach to
handle wheel-terrain dynamics.

## References

  * [1] V. Wiberg, E. Wallin, T. Nordfjell, and M. Servin, “Control of rough terrain vehicles using deep reinforcement learning,” _IEEE robotics and automation letters_ , vol. 7, no. 1, pp. 390–397, 2021.
  * [2] A. Gattupalli, V. P. Eathakota, A. K. Singh, and K. Madhava Krishna, “A simulation framework for evolution on uneven terrains for synchronous drive robot,” _Advanced Robotics_ , vol. 27, 2013.
  * [3] R. Agishev, T. Petříček, and K. Zimmermann, “Trajectory optimization using learned robot-terrain interaction model in exploration of large subterranean environments,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 3365–3371, 2022.
  * [4] E. Papadopoulos and D. A. Rey, “The force-angle measure of tipover stability margin for mobile manipulators,” _Vehicle System Dynamics_ , vol. 33, no. 1, pp. 29–48, 2000.
  * [5] R. Rubinstein, “The cross-entropy method for combinatorial and continuous optimization,” _Methodology and computing in applied probability_ , vol. 1, pp. 127–190, 1999.
  * [6] M.-R. Bouguelia, R. Gonzalez, K. Iagnemma, and S. Byttner, “Unsupervised classification of slip events for planetary exploration rovers,” _Journal of Terramechanics_ , vol. 73, pp. 95–106, 2017.
  * [7] S. Banerjee, J. Harrison, P. M. Furlong, and M. Pavone, “Adaptive meta-learning for identification of rover-terrain dynamics,” _arXiv preprint arXiv:2009.10191_ , 2020.
  * [8] A. Ugenti, F. Vulpi, A. Milella, and G. Reina, “Learning and prediction of vehicle-terrain interaction from 3d vision,” in _Multimodal Sensing and Artificial Intelligence: Technologies and Applications II_ , vol. 11785. SPIE, 2021, pp. 167–173.
  * [9] A. Datar, C. Pan, and X. Xiao, “Learning to model and plan for wheeled mobility on vertically challenging terrain,” _arXiv preprint arXiv:2306.11611_ , 2023.
  * [10] V. Šalanskỳ, K. Zimmermann, T. Petříček, and T. Svoboda, “Pose consistency kkt-loss for weakly supervised learning of robot-terrain interaction model,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 3, pp. 5477–5484, 2021.
  * [11] V. Eathakota, G. Aditya, and M. Krishna, “Quasi-static motion planning on uneven terrain for a wheeled mobile robot,” in _2011 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2011, pp. 4314–4320.
  * [12] Z. Xu, Y. Chen, Z. Jian, J. Tan, X. Wang, and B. L. Liang, “Hybrid trajectory optimization for autonomous terrain traversal of articulated tracked robots,” _IEEE Robotics and Automation Letters_ , vol. 9, no. 1, pp. 755–762, 2023.
  * [13] L. Xu, K. Chai, Z. Han, H. Liu, C. Xu, Y. Cao, and F. Gao, “An efficient trajectory planner for car-like robots on uneven terrain,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2023, pp. 2853–2860.
  * [14] W. Feng, L. Ding, R. Zhou, C. Xu, H. Yang, H. Gao, G. Liu, and Z. Deng, “Learning-based end-to-end navigation for planetary rovers considering non-geometric hazards,” _IEEE Robotics and Automation Letters_ , 2023.
  * [15] H. Hu, K. Zhang, A. H. Tan, M. Ruan, C. Agia, and G. Nejat, “A sim-to-real pipeline for deep reinforcement learning for autonomous robot navigation in cluttered rough terrain,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 4, pp. 6569–6576, 2021.
  * [16] K. Weerakoon, A. J. Sathyamoorthy, U. Patel, and D. Manocha, “Terp: Reliable planning in uneven outdoor environments using deep reinforcement learning,” in _International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 9447–9453.
  * [17] A. K. Singh and K. M. Krishna, “Feasible acceleration count: A novel dynamic stability metric and its use in incremental motion planning on uneven terrain,” _Robotics and Autonomous Systems_ , vol. 79, pp. 156–171, 2016\.
  * [18] S. Gould, R. Hartley, and D. Campbell, “Deep declarative networks,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 44, no. 8, pp. 3988–4004, 2021.
  * [19] Z. Han, Y. Wu, T. Li, L. Zhang, L. Pei, L. Xu, C. Li, C. Ma, C. Xu, S. Shen _et al._ , “An efficient spatial-temporal trajectory planner for autonomous vehicles in unstructured environments,” _IEEE Transactions on Intelligent Transportation Systems_ , 2023.
  * [20] D. Rey and E. Papadoupoulos, “Online automatic tipover prevention for mobile manipulators,” in _Proc. of the IEEE International Conference on Intelligent Robot and Systems_ , vol. 3, 1997, pp. 1273–1278.
  * [21] N. Chakraborty and A. Ghosal, “Kinematics of wheeled mobile robots on uneven terrain,” _Mechanism and machine theory_ , vol. 39, no. 12, pp. 1273–1287, 2004.
  * [22] A. L. Dontchev and R. T. Rockafellar, _Implicit functions and solution mappings_. Springer, 2009, vol. 543.
  * [23] “Jax,” https://github.com/google/jax.
  * [24] N. Koenig and A. Howard, “Design and use paradigms for gazebo, an open-source multi-robot simulator,” in _IEEE/RSJ international conference on intelligent robots and systems_ , 2004, pp. 2149–2154.
  * [25] “Blender,” https://www.blender.org/.


# Does an accelerated mirror suffer hindrance from vacuum?

E. Sadurní<EMAIL_ADDRESS>Instituto de Física, Benemérita Universidad
Autónoma de Puebla, Apartado Postal J-48, 72570 Puebla, México  M. A. Estévez
Instituto de Física, Benemérita Universidad Autónoma de Puebla, Apartado
Postal J-48, 72570 Puebla, México  J. L. Díaz-Cruz  Facultad de Ciencias
Físico Matemáticas, Benemérita Universidad Autónoma de Puebla, 72570 Puebla,
México

###### Abstract

We study the accelerated motion of mirrors and the photonic Unruh effect. The
solutions of Maxwell’s equations with appropriate boundary conditions in
Rindler coordinates are found. The canonical quantization of the field is
carried out properly. Important consequences arise from the presence of
mirrors or polarizers. It is shown that a single reflective surface produces
frequency quantization due to Dirichlet conditions applied to the $-1/(\alpha
x)^{2}$ anomalous potential for non-inertial observers. The number of photons
per accelerated mode is estimated via Hankel functions.

Accelerated mirrors in a flat background possess a notion of quantum vacuum
that differs significantly from that of inertial observers. Such is the lesson
that we have learnt from the Unruh effect, originally studied for scalar
particles unruh1976 ; crispino2008 . By studying the mechanism of particle
production –Bogoliubov transformation– we may conclude that the outcome should
be similar for other types of quantized fields, either bosonic or fermionic,
including Maxwell’s theory gill2010 ; narozhny2001 ; tulipant2021 ;
schleich2019 . Therefore, an accelerated mirror or polarizer interacting with
the inertial vacuum state must feel an effective radiation field, and it
should modify its motion as a reaction to radiation pressure. In other words,
an accelerated object collides with real photons, damping its motion. Similar
statements would apply when the mirror moves along the geodesics of a
background curved space, e.g. a blackhole.

The Unruh effect has spurred many experimental jiaz2019 ; uliana2019 ;
lynch2020  and theoretical carballo2019 ; cisco2020 ; ivette2011 ;
vanzella2017  investigations on field quantization in curved spaces, as well
as some polemical discussions on plausible detection schemes ford2006 . For
gravitational effects on interference patterns, see condado2018 . Although
these discussions delve deep into the very nature of particle detectors in
motion, it should be clear that quantum-mechanical averages also admit the
interpretation of a fixed –previously prepared– state and a transformed
observable that belongs to the non-inertial observer. In this case, the
observable is the particle number, and it is transformed together with the
corresponding detector to the non-inertial frame. Therefore, if an accelerated
observer starts abruptly its motion from rest and without visible radiation,
an effective state preparation at inertial vacuum will take place. Moreover,
if this happens at the beginning of our gedankenexperiment, the accelerated
trajectory will force the detectors to operate on the previous vacuum state
and not on the new one. The overall perception is that the transformed
observable $U^{\dagger}NU$ will produce non-vanishing averages $\langle
0|U^{\dagger}NU|0\rangle\neq 0$. In the present work we are interested in
studying the Unruh effect for photons interacting with mirrors, or polarizers,
and solving the corresponding mathematical intricacies of such scenario, which
to our knowledge has not been done before. In particular, we identify the
corresponding boundary conditions and present a consistent solution (exact and
approximated), which turns out to be non-trivial. We also discuss the physical
interpretation of our findings. In our exploration, we shall find the
existence of continuous and discrete frequency spectra of states localized by
$-1/x^{2}$ potentials, also known as quantum anomalous interactions with fall-
to-the-origin behavior tkachuk , shape invariance under scale transformations
serg  and Efimov states efimov . These potentials also appear naturally in
black hole horizons govin2000 . We shall verify this result in the
corresponding wave equation after the use of Rindler coordinates and the
Rindler metric, but with a slight variation in the case of photons as opposed
to scalars. In the former case, the radiation gauge can be applied to the
vector potential with the aim of eliminating polarization dependent
interactions. Boundary conditions describing the position of one or two
mirrors facing each other will also play an important role, as their presence
will modify the nature of the frequency spectrum, i.e. from continuous to
point-like. In this sense, accelerated Fabry-Pérot interferometers can be
regarded as photonic traps or bottles with a natural infrarred cutoff,
providing thus the ideal conditions for a detector. In addition, photons do
not have a natural length scale, such as $\lambda_{\rm Compton}=\hbar/mc$. As
a result, non-inertial observers need not surpass large critical accelerations
$\alpha_{c}$ in order to produce massless particles. The critical acceleration
scale is the overwhelming quantity $mc^{3}/\hbar$ (for the Higgs boson, the
acquired mass $m$ yields $\alpha_{c}\sim 10^{34}$m$/$s2 and a more
conservative example such as the pion yields $\alpha_{c}\sim 10^{31}$m$/$s2).

Figure 1:  (a) Frontal view of an accelerated mirror with polarization axis
$\hat{z}$ and perfect reflectivity along $\hat{x}$. The green signal
corresponds to the chirped wave given by Hankel functions; the blue signal
passes by. (b) and (c) Hankel functions for several values of $\omega/\alpha$,
green curves. A comparison between densities $|H^{(1)}_{iy}(ix)|^{2}$ and the
approximation (15) is included. (d) Frequency quantization (8) and (9) for one
and two mirrors.

We start our treatment with the following notation (units $\hbar=c=1$). The
Rindler coordinates of a right wedge involving $t$ and $x$ only, are given by

$\displaystyle T=x\sinh\alpha t,\quad X=x\cosh\alpha t,$ (1)

and the corresponding metric is

$\displaystyle g_{00}=\alpha^{2}x^{2},\quad g_{ij}=-\delta_{ij},\quad
g_{0i}=0.$ (2)

The world line of an accelerated object in $(T,X)$ spacetime is described by
$x=$ constant and $0<t<\infty$. In this respect, we shall always remember that
in this system, the low acceleration limit is reached by letting
$\alpha\rightarrow 0$ as well as $x\rightarrow 1/\alpha$, which corresponds to
the asymptotically flat region $\alpha x=1$. Some authors prefer the
redefinition $x\mapsto x+1/\alpha$ in order to recover the flat limit
seamlessly, but this is merely a translation in $x$. The Klein-Gordon equation
has at least two equivalent forms (here, the scalar curvature is $R=0$)

$\displaystyle\left[(|g|^{-1/2})\partial_{\mu}(|g|^{1/2}g^{\mu\nu}\partial_{\nu})+m^{2}\right]\phi=0,$
$\displaystyle\left[\square+M^{2}(x)\right]\psi=0,$ (3)

where $\square\equiv\partial_{\mu}g^{\mu\nu}\partial_{\nu}$ (not to be
confused with $\nabla_{\mu}\nabla^{\mu}$) and $\psi=|g|^{1/4}\phi$. The second
form can be useful, as it allows to understand the acceleration effect as a
redefinition of the mass $M^{2}=m^{2}-1/(2x)^{2}$, while the wave $\psi$
incorportates automatically the square root of the Jacobian. The wave equation
coming from Maxwell’s theory suffers a slightly different modification due to
the new curvilinear coordinates, as it deals with a vector field instead of a
scalar. We have an extra term coupled to polarization in the following
equation

$\displaystyle\frac{1}{\sqrt{g}}\partial_{\mu}\left(\sqrt{g}\partial^{\mu}A_{\nu}\right)-\Gamma^{\sigma}_{\mu\nu}\partial^{\mu}A_{\sigma}=0,$
(4)

where we have made the choice of a covariant Lorentz gauge, i.e.
$\nabla_{\mu}A^{\mu}=0=(1/\sqrt{g})\partial_{\mu}(\sqrt{g}A^{\mu})$. Now we
solve (4) by separation of variables; the Lorentz gauge helps to eliminate the
$\Gamma$\- dependent polarization terms. In Rindler spacetime
$\Gamma^{1}_{00}=\alpha x$, $\Gamma^{0}_{10}=\Gamma^{0}_{01}=1/(\alpha x)$ and
all other Christoffel symbols vanish. Therefore

$\displaystyle-\Gamma^{\sigma}_{\mu\nu}\partial^{\mu}A_{\sigma}=\left(-\alpha
x\partial^{0}A_{1}-\frac{\partial^{1}A_{0}}{\alpha
x},-\frac{\partial^{0}A_{0}}{\alpha x},0,0\right).$ (5)

Meanwhile, separable solutions of the form $A_{\mu}(x_{\nu})=e^{-i\omega
t}F_{\mu}(x)e^{i{\bf k}_{\perp}\cdot{\bf x}}$ can be proposed. In fig. 1(a),
we introduce the notation ${\bf k}_{\rm net}={\bf k}_{\perp}+{\bf
k}_{\parallel}$ as a total wave vector, where ${\bf k}_{\perp}\cdot{\bf
x}=k_{y}y+k_{z}z$ is perpendicular to the accelerated motion and ${\bf
k}_{\parallel}=\hat{x}k_{\parallel}$ is defined along the mirror’s trajectory,
with magnitude $k_{\parallel}\equiv\sqrt{\omega^{2}-k_{\perp}^{2}}$. For
linear polarizers along ${\bf k}_{\perp}$, a field polarization vector
$\hat{\varepsilon}\parallel\hat{y}$ suffers reflection, which is the situation
of interest.

Upon substitution of $A_{\nu}$ into
$\partial_{\mu}(\sqrt{g}g^{\mu\nu}A_{\nu})=0$ we find the usual transverse
condition $A_{0}=A_{1}=0$ (no longitudinal photons for rectilinear accelerated
motion) and the field is reduced to

$\displaystyle A_{2}=A_{2}(0)e^{-i\omega t}e^{i{\bf k}_{\perp}\cdot{\bf
x}}F(x),$ $\displaystyle A_{3}=A_{3}(0)e^{-i\omega t}e^{i{\bf
k}_{\perp}\cdot{\bf x}}F(x),$ $\displaystyle{\bf k}_{\perp}\cdot{\bf
A}=k_{y}A_{2}(0)+k_{z}A_{3}(0)=0.$ (6)

As expected, the polarization of ${\bf A}$ (and ${\bf E}$) in the $(y,z)$
plane is also orthogonal to oblique wave vectors for this set of solutions, as
shown in the green chirped pulses in fig. 1(a). It is important to note that
there is another set of solutions for which $A_{1}\neq 0$, leading to the
gauge restriction $\partial_{1}(xA_{1})-ix{\bf k}\cdot{\bf A}=0$ and a non-
vanishing $\Gamma$ term in (4). These solutions correspond to polarizations in
the $x$ axis. If a beam splitter with perfect reflectivity for $y$-polarized
waves is accelerated, only (6) will have the chance to interact with it. The
other type of photons will pass by undetected, and they shall be ignored in
what follows for simplicity.

Finally, we can reduce Maxwell’s wave equation for (6) to the massless K.-G.
theory in curved space for each component of $A$, i.e.
$\left[\square-1/(4x^{2})\right]|g|^{1/4}A_{\mu}=0$. With this, the
undetermined function $F(x)$, defined through $\psi=|g|^{1/4}F$, is shown to
satisfy

$\displaystyle\left[-\frac{d^{2}}{dx^{2}}-\frac{(\omega/\alpha)^{2}+1/4}{x^{2}}\right]\psi=-k_{\perp}^{2}\psi.$
(7)

Indeed, the quantum anomalous potential along the direction of motion appears
in (7). It is identical to a radial problem $-1/r^{2}$ for S waves; many
important mathematical results have been obtained for this kind of eigenvalue
problem. They are indispensable, as the completeness of the accelerated modes
is required in the computation of relativistic inner products, which in turn
play the role of coefficients in the Bogoliubov transformation. One of such
properties is the scale invariance of the eigenvalue equation under
$x\mapsto\lambda x,k_{\perp}\mapsto k_{\perp}/\lambda$, for any $\lambda\neq
0$, leading to a continuous spectrum in the absence of additional boundary
conditions.

A more delicate analysis begins here: It is well known that $F$ can be put in
terms of Bessel functions with the help of the definitions $\nu\equiv
i\omega/\alpha$, $\zeta\equiv ik_{\perp}x$ for the order $\nu$ and the
argument $\zeta$ (both are imaginary). The Hankel functions are better suited
for this study, as we can identify immediately that
$H_{\nu}^{(1)}(ik_{\perp}x)$ satisfies the condition $|F|\rightarrow 0$ as
$x\rightarrow\infty$, which is the only possibility for appropriate quantum-
mechanical interpretations. This boundary condition at infinity does not
quantize the frequency, as we can see from the scaling properties above. The
effect of this restriction reduces to the elimination of
$H^{(2)}_{\nu}(ik_{\perp}x)$, for it grows exponentially as
$x\rightarrow\infty$.

Now we arrive at the crucial point; if a mirror accelerates in vacuo, the
world line $x=x_{\rm m}=$ constant imposes a boundary condition on the
photon’s wave function precisely at this location. A mirror may interact with
the produced radiation field in many forms, but if we think of this detector
as a scatterer, there will be a resonance of the photon field when the phase
shift of $\pi$ is attained. The reflective properties of a polarizer are not
too different, if a specific polarization is chosen. The wave function at rest
is $\sin\left[q(x-x_{\rm m})\right]$ for $x\geq x_{m}$ and $0$ otherwise. For
the accelerated observer, the Dirichlet boundary condition still holds, but
the wave function is now $\psi=\mbox{$\cal N\,$}H^{(1)}_{\nu}(ik_{\perp}x)$,
so we must have the following exact  quantization condition:

$\displaystyle H^{(1)}_{i\omega/\alpha}(ik_{\perp}x_{\rm
m})=0,\,\,k_{\perp}x_{\rm m}=r_{n}(\nu),\,\,n=0,1,2,...$ (8)

where $r_{n}(\nu)$ are the (real) roots of the Hankel function of imaginary
variable and index. A density plot in fig. 1(b) and (c) reveals the shape of
these functions.

In this kind of anomalous potential, it is even more striking that a second
boundary condition does not lead to an overspecified system. A photonic
’bottle’ can be built with two mirrors facing each other at positions
$x_{1},x_{2}$. Only certain modes can be established inside of this device,
which corresponds to a Fabry-Pérot interferometer in motion (initially empty).
The reason for frequency quantization is the allowed presence of $H^{(2)}$, as
the full wavefunction vanishes identically outside of the bottle and becomes
square integrable in the domain $x_{1}<x<x_{2}$. The eigenmodes and exact
quantization condition are, respectively,

$\displaystyle\psi$ $\displaystyle=$ $\displaystyle\mbox{$\cal
N\,$}\sqrt{x}\left[H_{i\omega/\alpha}^{(1)}(ik_{\perp}x_{1})H_{i\omega/\alpha}^{(2)}(ik_{\perp}x)\right.$
$\displaystyle-$
$\displaystyle\left.H_{i\omega/\alpha}^{(2)}(ik_{\perp}x_{1})H_{i\omega/\alpha}^{(1)}(ik_{\perp}x)\right],$
$\displaystyle 0$ $\displaystyle=$ $\displaystyle
H_{i\omega/\alpha}^{(1)}(ik_{\perp}x_{1})H_{i\omega/\alpha}^{(2)}(ik_{\perp}x_{2})$
(9) $\displaystyle-$ $\displaystyle
H_{i\omega/\alpha}^{(2)}(ik_{\perp}x_{1})H_{i\omega/\alpha}^{(1)}(ik_{\perp}x_{2}).$

Now we proceed to quantize the e.m. field that interacts with the mirror.
Proper care must be taken in the discussion of polarization and Lorentz
invariance (Gupta-Bleuler quantization). Firstly, from the superposition of
solutions parallel and orthogonal to the mirror, i.e. ${\bf A}_{\rm
total}={\bf A}_{\parallel}+{\bf A}_{\perp}$, we make sure that the field
bouncing off the mirror (or inside the interferometer) is ${\bf
A}_{\parallel}\equiv{\bf A}$. Secondly, the electric field is recovered here
with ${\bf E}=-\partial_{t}{\bf A}$ after $A_{0}=0$ has been chosen, so the
field ’coordinate’ and the field ’momentum’ share the same polarization.
Thirdly, the travelling mirror breaks both spatial isotropy and global Lorentz
covariance (this is not Minkowski). Therefore, we employ (6) alone in the
canonical quantization scheme, i.e. only one polarization
$\hat{\varepsilon}\perp{\bf k}_{\perp}$ interacts with the mirror and the
other set of solutions for which $A_{1}\neq 0$ becomes irrelevant. We write
the commutators at equal times as

$\displaystyle\left[\hat{\varepsilon}\cdot{\bf A}({\bf
x},t),g_{00}^{-1/2}\partial_{t}\hat{\varepsilon}\cdot{\bf A}({\bf
x}^{\prime},t)\right]$ $\displaystyle=$ $\displaystyle i\delta^{3}({\bf
x}-{\bf x}^{\prime})$ $\displaystyle g^{00}\left[\hat{\varepsilon}\cdot{\bf
A}({\bf x},t),\hat{\varepsilon}\cdot{\bf A}({\bf x}^{\prime},t)\right]$
$\displaystyle=$ $\displaystyle 0.$ (10)

The mode decomposition of these operators rests on completeness relations, so
first we define the ’free’ modes as plane waves with appropriate metric
normalization: $g^{1/4}e^{i{\bf k}\cdot{\bf x}}/(2\pi)^{3/2}$. Then we
incorporte boundary conditions with the combination
$\phi_{q}=g^{1/4}e^{-i\omega_{q}t+i{\bf q}_{\perp}\cdot{\bf
x}}\sin\left[q(x-x_{\rm m})\right]$, $\omega_{q}=\sqrt{q^{2}+q_{\perp}^{2}}$
and ${\bf q}_{\perp}=q_{y}\hat{y}+q_{z}\hat{z}$, all this for a single mirror.
Using
$\partial^{0}=g^{00}\partial_{0}=(g_{00})^{-1}\partial_{0}=(|g|)^{-1}\partial_{0}$
for the Rindler metric and
$(d^{3}x\sqrt{|g|})(\phi\partial^{0}\psi)=(d^{3}x/\sqrt{|g|})(\phi\partial_{0}\psi)$
for the volume element, the orthogonality of the following relativistic inner
product is ensured

$\displaystyle\langle\phi_{q},\phi_{q^{\prime}}\rangle$ $\displaystyle=$
$\displaystyle i\int
d^{3}x\sqrt{|g|}(\phi_{q}^{*}\partial^{0}\phi_{q^{\prime}}-\phi_{q^{\prime}}\partial^{0}\phi_{q}^{*})$
$\displaystyle=$
$\displaystyle-\frac{i}{(2\pi)^{3}}\int\frac{d^{3}x}{\sqrt{|g|}}|g|^{1/4}|g|^{1/4}\times$
$\displaystyle\times$
$\displaystyle(i\omega_{q^{\prime}}+i\omega_{q})e^{i({\bf q}_{\perp}-{\bf
q}^{\prime}_{\perp})\cdot{\bf x}}\times$ $\displaystyle\times$
$\displaystyle\sin\left[q(x-x_{\rm m})\right]\sin\left[q^{\prime}(x-x_{\rm
m})\right]$ $\displaystyle=$ $\displaystyle 2\omega_{q}\delta^{3}({\bf q}-{\bf
q}^{\prime}).$ $\displaystyle\langle\phi_{k},\phi^{*}_{\pm k^{\prime}}\rangle$
$\displaystyle=$ $\displaystyle 0.$ (11)

For quantized frequencies $\omega_{n}$ of the accelerated eigenmodes, we have
$\psi_{n,k_{\perp}}\equiv e^{-i\omega_{n}t+i{\bf k}_{\perp}\cdot{\bf
x}}\sqrt{k_{\perp}x}H^{(1)}_{i\omega/\alpha}(ik_{\perp}x)=e^{-i\omega_{n}t+i{\bf
k}_{\perp}\cdot{\bf x}}\sqrt{\alpha x}\chi_{n}(x)$, where the specific form of
$\omega_{n}(k_{\perp})$ will be discussed later. The inner product yields
$\langle\psi_{n},\psi_{m}\rangle=2\omega_{n}(k_{\perp})\delta_{n,m}\delta^{2}({\bf
k}_{\perp}-{\bf k}^{\prime}_{\perp})$, and for negative energies or
frequencies $\langle\psi_{n},\psi_{m}^{*}\rangle=0$. Finally, the Bogoliubov
transformation from inertial operators $\hat{A}$ to non-inertial operators
$\hat{B}$ involves both creation and annhilation of particles as

$\displaystyle\hat{B}_{n,k_{\perp}}=\int d^{3}q$
$\displaystyle\left[\left(\frac{\langle\psi_{n,k_{\perp}},\phi_{q}\rangle^{*}}{2\sqrt{\omega_{q}\omega_{n}}}\right)\hat{A}_{q}\right.$
(12) $\displaystyle+$
$\displaystyle\left.\left(\frac{\langle\psi^{*}_{n,k_{\perp}},\phi_{q}\rangle}{2\sqrt{\omega_{q}\omega_{n}}}\right)\hat{A}^{\dagger}_{q}\right].$

The non-vanishing contribution of the transformed average number is controlled
by the following overlap

$\displaystyle\langle\phi_{q}^{*},\psi_{n,k_{\perp}}\rangle$ $\displaystyle=$
$\displaystyle
e^{i(\omega_{q}-\omega_{n})t}\left[\frac{\omega_{n}-\omega_{q}}{(2\pi)^{3/2}}\right]\times$
$\displaystyle\times$ $\displaystyle\int d^{3}xe^{i({\bf k}_{\perp}-{\bf
q}_{\perp})\cdot{\bf x}}\sin\left[q(x-x_{\rm m})\right]\chi_{n}(x)$
$\displaystyle=$ $\displaystyle
e^{i(\omega_{q}-\omega_{n})t}\left[\omega_{n}-\omega_{q}\right]\delta^{2}({\bf
k}_{\perp}-{\bf q}_{\perp})F_{n,q},$ $\displaystyle F_{n,q}$
$\displaystyle\equiv$ $\displaystyle\frac{1}{\sqrt{2\pi}}\int_{x_{\rm
m}}^{\infty}\sin\left[q(x-x_{\rm m})\right]\chi_{n}(x),$ (13)

where the complete set $\psi_{n,k_{\perp}}=e^{-i\omega_{n}t+i{\bf
k}_{\perp}\cdot{\bf x}}g^{1/4}\chi_{n}(x)$ solves the 1D problem with
anomalous potential and mirror boundary. One can see that $\chi_{n}$ is
proportional to the Hankel function $H^{(1)}$; therefore, the remaining
integral $F_{n,q}$ is the sine-Fourier transform of the corresponding Hankel
function subject to (8). The expectation value of the number operator between
inertial vacuum states $|0\rangle_{A}$ is not zero, which can be shown if we
employ $\hat{A}_{q^{\prime}}\hat{A}_{q}^{\dagger}=\delta^{3}({\bf q}-{\bf
q}^{\prime})+\hat{A}_{q}^{\dagger}\hat{A}_{q^{\prime}}$ in the expression for
$\langle N\rangle$:

$\displaystyle\langle
0|B_{n,k_{\perp}}^{\dagger}B_{n,k_{\perp}}|0\rangle_{A}=\int d^{3}q\int
d^{3}q^{\prime}\delta^{2}({\bf q}^{\prime}_{\perp}-{\bf k}_{\perp})\times$
(14) $\displaystyle\times$ $\displaystyle\delta^{2}({\bf q}_{\perp}-{\bf
k}_{\perp})\delta^{3}({\bf q}-{\bf
q}^{\prime})\frac{|F_{n,q}|^{2}\left[\omega_{n}-\omega_{q}\right]^{2}}{4\omega_{n}\omega_{q}}$
$\displaystyle=$ $\displaystyle V_{k_{\perp}}\int
dq\frac{|F_{n,q}|^{2}\left[\omega_{n}-\omega_{q}\right]^{2}}{4\omega_{n}\omega_{q}}\bigg{\rvert}_{{\bf
q}_{\perp}={\bf k}_{\perp}},$

where $V_{k_{\perp}}$ is the volume of all available states ${\bf k}_{\perp}$,
thus avoiding the evaluation of $\delta(0)$ as costumary. A variety of
techniques can be employed in the estimation of the integrals involved in this
expression. The Fourier transform of the Hankel function is best estimated by
means of the stationary phase approximation applied to its integral
representation (grad , pg. 915, sec. (8.421) formula 8):

$\displaystyle H_{\nu}^{(1)}(xz)$ $\displaystyle=$
$\displaystyle\frac{z^{\nu}e^{-i\nu\pi/2}}{i\pi}\int_{0}^{\infty}dte^{ix(t+z^{2}/t)/2}t^{-\nu-1}$
(15) $\displaystyle=$
$\displaystyle\frac{1}{i\pi}\int_{-\infty}^{\infty}du\quad e^{-\nu
u+i|z|x\sinh u}$ $\displaystyle\approx$
$\displaystyle\sqrt{\frac{8i\pi}{\sqrt{|\nu|^{2}-(|z|x)^{2}}}}\cos\left[|\nu|{\rm
arccosh}\left(\frac{|\nu|}{|z|x}\right)\right.$ $\displaystyle-$
$\displaystyle\left.\sqrt{|\nu|^{2}-(|z|x)^{2}}-\pi/4\right]\Theta\left(\left|\frac{\nu}{z}\right|-x\right).$

In passing, we note that this estimate resembles closely the WKB approximation
of the 1D problem with anomalous $-1/x^{2}$ potential, but the singularity at
the turning point
$x=\sqrt{(\omega/\alpha)^{2}+1/4}/k_{\perp}\approx\omega/(k_{\perp}\alpha)$
can be ’cured’ in the limit of small $\alpha$ and large $x\sim 1/\alpha$. We
have
$\sqrt{|\nu|-|z|^{2}x^{2}}\approx(\sqrt{\omega^{2}-k_{\perp}^{2}})/\alpha=k_{\parallel}/\alpha$.
In the argument of the cosine we have $x\sqrt{(|\nu|/x)^{2}-|z|^{2}}\approx
x\sqrt{\omega^{2}-k_{\perp}^{2}}=k_{\parallel}x$. Also in this limit, we have
${\rm arccosh}(|\nu|/|zx|)\approx\log(2|\nu|/|zx|)$. Regarding the
normalization of bound solutions, we shall preserve the Hankel function’s
Dirac-delta convention satisfied by $\sqrt{zx}H_{\nu}^{(1)}(zx)$ instead of a
Kronecker delta normalization of the square integrable function, because the
limit $\alpha\rightarrow 0$, $x\rightarrow 1/\alpha\rightarrow\infty$ would
make collapse all solutions $\psi_{n}$ to zero otherwise. Also, the argument
of the cosine in (15) must be chosen such that $\psi_{n}(x_{\rm m})=0$. With
these considerations, we have the following wave functions:

$\displaystyle\psi_{n}$ $\displaystyle\approx$ $\displaystyle\sqrt{\frac{8\pi
ik_{\perp}\alpha x}{k_{\parallel}}}\sin\left[k_{\parallel}(x-x_{\rm
m})+\frac{\omega_{n}}{\alpha}\log\left(\frac{x}{x_{\rm
m}}\right)\right]\times$ (16) $\displaystyle\times$
$\displaystyle\Theta\left(\frac{\omega_{n}}{\alpha k_{\perp}}-x\right)$

while the implicit determination of the frequency function
$\omega_{n}(k_{\perp})$ comes from the quantization condition (8) applied to
(15):

$\displaystyle\left(\frac{\omega}{\alpha}\right){\rm arccosh\
}\left(\frac{\omega}{\alpha k_{\perp}x_{\rm
m}}\right)-\sqrt{\left(\frac{\omega}{\alpha}\right)^{2}-k_{\perp}^{2}x_{\rm
m}^{2}}=(n+1/4)\pi,$ $\displaystyle\mbox{with}\quad n=0,1,2,...$ (17)

A reasonable approximation for the frequency emerging from this relation is
$\omega_{n}\approx\sqrt{(\alpha x_{\rm
m}\pi(n+1/4)k_{\perp})^{2}+k_{\perp}^{2}}$, thus quantizing
$k_{\parallel}\approx k_{\perp}\alpha x_{\rm m}\pi(n+1/4)$, see fig. 1(d).
With these approximations, the overlap $F_{n,q}$ can be put in terms of Gamma
functions explicitly. For our purposes, it is better to estimate $F_{n,q}$
using a stationary phase approximation:

$\displaystyle|F_{n,q}|^{2}=\frac{k_{\perp}\omega_{n}}{k_{\parallel}|q-k_{\parallel}|\alpha}\times$
$\displaystyle\times$
$\displaystyle\Theta(q-k_{\parallel}-k_{\perp})\Theta\left(k_{\parallel}+\frac{\omega}{\alpha
x_{\rm m}}-q\right)\times$ $\displaystyle\times$
$\displaystyle\cos^{2}\left[\frac{\omega_{n}}{\alpha}\left\\{1+\log\left(\frac{\omega_{n}}{\alpha
x_{\rm m}|q-k_{\parallel}|}\right)\right\\}-|q-k_{\parallel}|x_{\rm
m}\right].$

The computation of the limit $\alpha\rightarrow 0$ in this expression can be
tricky, as $F_{n,q}\rightarrow\delta(q-k_{\parallel})$. This is important in
the determination of (14) for low accelerations: The average photon number per
mode satisfies $\langle N\rangle\rightarrow 0$ as $\alpha\rightarrow 0$, which
is verified by substituting
$|F_{n,q}\times(\omega_{n}-\omega_{q})|^{2}\rightarrow|\delta(q-k_{\parallel})(\omega_{n}-\omega_{q})|^{2}=0$.
We also note that the location of the stationary point in the integral
defining $F_{n,q}$ must satisfy

$\displaystyle x^{*}=\frac{\omega_{n}}{\alpha(q-k_{\parallel})},\quad x_{\rm
m}<x^{*}<\frac{\omega_{n}}{\alpha k_{\perp}}$ (19)

and it produces a natural infrared cutoff through the second inequality above.
This is expressed in the first Heaviside Theta function in (LABEL:e15), since
$q_{\rm min}=k_{\parallel}+k_{\perp}$ is the minimal parallel momentum of the
photon. This also renders a nonsingular prefactor $1/(q-k_{\parallel})$ in
(LABEL:e15). The second step function in (LABEL:e15) shows a natural
ultraviolet cutoff $q_{\rm max}=k_{\parallel}+\omega_{n}/\alpha x_{\rm m}$,
coming from the first inequality in (19). These cuts in the wavenumber also
produce minimal and maximal photon frequencies $\omega_{\rm min}=\sqrt{q_{\rm
min}^{2}+k_{\perp}^{2}}$, $\omega_{\rm max}=\sqrt{q_{\rm
max}^{2}+k_{\perp}^{2}}$ respectively.

The integral over $q$ for the total average number per excited mode converges,
and it is composed of an oscillatory (trigonometric) part and a monotonic
function, as seen from the simple trick $\cos^{2}\vartheta=(1+\cos
2\vartheta)/2$. We retain the strongest contribution from the non-oscillatory
part of the integrand and find the following estimate:

$\displaystyle\langle N\rangle$ $\displaystyle\approx$
$\displaystyle\frac{V_{k_{\perp}}k_{\perp}}{8k_{\parallel}}\left[\omega_{\rm
max}-\omega_{\rm min}+k_{\parallel}\log\left(\frac{q_{\rm max}+\omega_{\rm
max}}{q_{\rm min}+\omega_{\rm min}}\right)\right.$ (20) $\displaystyle-$
$\displaystyle\left.2\omega_{n}\log\left(\frac{k_{\perp}^{2}+k_{\parallel}q_{\rm
max}+\omega_{n}\omega_{\rm max}}{k_{\perp}^{2}+k_{\parallel}q_{\rm
min}+\omega_{n}\omega_{\rm min}}\right)\right].$

Some important conclusions can be drawn from our results. Firstly, the typical
Boltzmann factor or Bose-Einstein distribution in the occupation numbers as
functions of the energy are absent here; usually one finds such exponentials
(and the corresponding effective temperatures in terms of $\alpha$) by
computing the Bogoliubov coefficients through complex integration techniques
in the absence of boundary conditions. Here, on the other hand, we have an
important perturbation produced by the mirror. For very weak accelerations, we
must consider mirrors located at $x\sim 1/\alpha$ and frontal photons such
that $k_{\parallel}\gg k_{\perp}$. Then $\langle N\rangle\sim
ak_{\perp}/\alpha+bk_{\perp}^{2}/\alpha k_{\parallel}$, with $a,b$ numerical
constants. This means that the number decreases as $1/k_{\parallel}$, i.e.
inversely proportional to photon’s momentum. However, there is a remnant that
vanishes linearly as $k_{\perp}\rightarrow 0$. No threshold mass or critical
energy appeared in our calculations, as expected. As to the mirror’s motion,
it is clear now that the momentum imprinted on this object by the field is
proportional to $\langle\hat{x}\cdot{\bf E}\times{\bf B}\rangle$, which yields
a superposition of all occupation numbers weighed by $k_{\parallel}$. Since
the number per mode depends on $\alpha$ via $\omega_{n}$, the change of the
mirror’s linear momentum depends on the acceleration. Damping sets in, and
there appears an effect similar to the Abraham-Lorentz force for classically
accelerated charges, except that the mirror can be electrically neutral! A
concave mirror could act as a parachute, and help to manoeuvre a spaceship, in
the event that it were ever possible to travel near the vicinity of a black
hole. Plus the energy employed in the sudden acceleration of a photonic bottle
can be converted into trapped radiation, a most attractive possibility.

###### Acknowledgements.

Financial support from VIEP project 100518931 BUAP-CA-289 is acknowledged.

## References

  * (1) W. G. Unruh, Phys. Rev.D 14, 870 (1976)
  * (2) L. C. B. Crispino, Atsushi Higuchi, and George E. A. Matsas, Rev. Mod. Phys. 80, 787 (2008)
  * (3) A. de Gill, American Journal of Physics 78, 685 (2010)
  * (4) N. B. Narozhny, A. M. Fedotov, B. M. Karnakov, V. D. Mur, and V. A. Belinskii, Phys. Rev. D 65, 025004 (2001)
  * (5) Z. Tulipant, The European Physical Journal C 81, 288 (2021)
  * (6) M. O. Scully, S. Fulling, D. M. Lee, D. N. Page, W. P. Schleich, and A. A. Svidzinsky, PNAS 115 (32) 8131–8136 (2018)
  * (7) M. H. Lynch, E. Cohen, Y. Hadad, and I. Kaminer, Experimental Evidence for the Unruh Effect, CLEO: QELS Fundamental Science 2020, Washington, DC United States, 10–15 May 2020. Page FM2C.5
  * (8) C. A. Uliana Lima, F. Brito, J. A. Hoyos and D. A. Turolla Vanzella, Nature Communications 10, 3030 (2019)
  * (9) Jiazhong Hu et al. Nature Physics (2019). DOI: 10.1038/s41567-019-0537-1
  * (10) C. Gooding, S. Biermann, S. Erne, J. Louko, W. G. Unruh, J. Schmiedmayer, and S. Weinfurtner, Phys. Rev. Lett 125 (21): 213603 (2020); doi: 10.1103/PhysRevLett.125.213603
  * (11) E. Martín-Martínez, I. Fuentes and R. B. Mann, Phys. Rev. Lett. 107, 131301 (2011)
  * (12) G. Cozzella, A. G. S. Landulfo, G. E. A. Matsas, and D. A. T. Vanzella, Phys. Rev. Lett. 118, 161102 (2017)
  * (13) R. Carballo-Rubio, L. J. Garay, E. Martín-Martínez, and J. de Ramón, Phys. Rev. Lett. 123, 041601 (2019)
  * (14) G.W. Ford, R.F. O’Connell, Physics Letters A 350, Issues 1–2 (2006)
  * (15) D. Condado, J. L. Díaz-Cruz, A. Rosado and E. Sadurní, Phys. Rev. A 98, 043618 (2018)
  * (16) V. M. Vasyutaa and V. M. Tkachuk, Eur. Phys. J. D (2016) 70: 267
  * (17) S. Moroz, and R. Schmidt, Annals of Physics 325, 2 (2010) 491–513
  * (18) V. N. Efimov, Sov. J. Nucl. Phys. 12, 589 (1971)
  * (19) T. R. Govindarajan, V. Suneeta, S. Vaidya, Nuclear Physics B 583, Issues 1–2, (2000) 291–303
  * (20) I. S. Gradshteyn and I. M. Ryzhik, Tables of Series, Integrals and Products, seventh ed. Academic Press, Amsterdam 2007


# Non-Existence and Existence of Shock Profiles
in the Bemfica-Disconzi-Noronha Model

Heinrich Freistühler  Department of Mathematics, University of Konstanz, 78457
Konstanz, Germany. Supported by DFG Grant No. FR 822/10-1.

(March 27, 2021)

###### Abstract

This note studies a four-field hyperbolic PDE model that was recently
introduced by Bemfica, Disconzi, and Noronha for the pure radiation fluid with
viscosity, and asks whether shock waves admit continuous profiles in this
description. The model containing two free parameters $\mu,\nu$ and being
causal whenever one chooses $(\mu,\nu)$ from a certain range
$\mathcal{C}\subset\mathbb{R}^{2}$, this paper shows that for any choice of
$(\mu,\nu)$ in the interior of $\mathcal{C}$, there is a dichotomy in so far
as (i) shocks of sufficiently small amplitude admit profiles and (ii) certain
other, thus necessarily non-small, shocks do not. This finding does not
preclude the possibility that if one chooses $(\mu,\nu)$ from a specific part
$\mathcal{S}$ of the boundary of $\mathcal{C}$, the dichotomy disappears and
_all_ shocks have profiles; the parameter set $\mathcal{S}$ corresponds to the
“sharply causal” case, in which one of the characteristic speeds of the
dissipation operator is the speed of light.

##  1 Introduction

In their admirable 2018 paper [1], Bemfica, Disconzi, and Noronha have
proposed a four-field PDE formulation

$\displaystyle\partial_{\beta}(T^{\alpha\beta}+\Delta T^{\alpha\beta})=0$
(1.1)

for the dynamics of the pure radiation fluid, ideally

$\displaystyle
T^{\alpha\beta}=\frac{\partial(p(\theta)\psi^{\alpha})}{\partial\psi_{\beta}}=\theta^{3}p^{\prime}(\theta)\psi^{\alpha}\psi^{\beta}+p(\theta)g^{\alpha\beta},\quad\psi^{\gamma}=\frac{u^{\gamma}}{\theta},\quad
p(\theta)=\frac{1}{3}\theta^{4},$ (1.2)

with a dissipation tensor

$\displaystyle\Delta
T^{\alpha\beta}=-B^{\alpha\beta\gamma\delta}(\psi)\displaystyle{\frac{\partial\psi_{\gamma}}{\partial
x^{\delta}}}$ (1.3)

in which the classical Eckart viscosity ansatz [5]

$\displaystyle\begin{aligned}
B^{\alpha\beta\gamma\delta}_{visc}&=\Pi^{\alpha\gamma}\Pi^{\beta\delta}+\Pi^{\alpha\delta}\Pi^{\beta\gamma}-\frac{2}{3}\Pi^{\alpha\beta}\Pi^{\gamma\delta}\quad\text{with}\quad\Pi^{\alpha\gamma}=u^{\alpha}u^{\beta}+g^{\alpha\beta}\end{aligned}$
(1.4)

is augmented as111Our notation is slightly different from theirs. We write
$B_{ther},B_{velo}$ for the two characteristic pieces of the augmentation in
order to make their relation to first-order equivalence transformations [12]
apparent, and use the thermodynamically natural Godunov variable
$\psi^{\gamma}$ as in [18, 11, 6].

$\displaystyle B^{\alpha\beta\gamma\delta}=\eta
B^{\alpha\beta\gamma\delta}_{visc}-\mu B^{\alpha\beta\gamma\delta}_{ther}-\nu
B^{\alpha\beta\gamma\delta}_{velo}$ (1.5)

with

$\displaystyle
B^{\alpha\beta\gamma\delta}_{ther}=(3u^{\alpha}u^{\beta}+\Pi^{\alpha\beta})(3u^{\gamma}u^{\delta}+\Pi^{\gamma\delta}),\quad
B^{\alpha\beta\gamma\delta}_{velo}=(u^{\alpha}{\Pi^{\beta}}_{\epsilon}+u^{\beta}{\Pi^{\alpha}}_{\epsilon})(u^{\gamma}\Pi^{\delta\epsilon}+u^{\delta}\Pi^{\gamma\epsilon}).$
(1.6)

Like that of a different augmentation introduced by B. Temple and the author
in [11, 12], the purpose of the combination (1.5) is to make the dissipation
causal and thus remedy a well-known deficiency of the Eckart theory.222Recall
Israel’s work [14] for an early definitive statement of this deficiency and
the first stringent attempt to handle it. Many others followed, cf. the
introductions of [11] and [6]. It was shown in [1] that this formulation,
which we will here briefly refer to as ‘the BDN model’, is indeed causal if
and only if, relative to the classical coefficient $\eta$ of viscosity, the
coefficients $\mu$ and $\nu$ of the “regulators”
$B^{\alpha\beta\gamma\delta}_{ther},B^{\alpha\beta\gamma\delta}_{velo}$
satisfy

$\displaystyle\mu\geq\frac{4}{3}\eta\quad\text{and}\quad\nu\leq\left(\frac{1}{3\eta}-\frac{1}{9\mu}\right)^{-1}.$
(1.7)

The present note focusses on shock waves, whose ideal version is given by
discontinuous solutions to the corresponding inviscid, Euler, equations

$\displaystyle\partial_{\beta}T^{\alpha\beta}=0$ (1.8)

of the (prototypical) form

$\displaystyle\psi(x)=\begin{cases}\psi_{-},&x^{\beta}\xi_{\beta}<0,\\\
\psi_{+},&x^{\beta}\xi_{\beta}>0,\end{cases}$ (1.9)

and asks whether they can be properly represented in the viscous setting
(1.1). A standard way to achieve such representation of a ‘viscous shock wave’
is a ‘dissipation profile’, i.e., a regular solution of (1.1) that depends
also only on $x^{\beta}\xi_{\beta}$ and connects the two states forming the
shock, in other words, a solution of the ODE

$\displaystyle\xi_{\beta}\xi_{\delta}B^{\alpha\beta\gamma\delta}(\psi)\psi_{\gamma}^{\prime}=\xi_{\beta}T^{\alpha\beta}(\psi)-q^{\alpha},\quad
q^{\alpha}:=\xi_{\beta}T^{\alpha\beta}(\psi_{\pm}),$ (1.10)

on $\mathbb{R}$ which is heteroclinic to them,

$\displaystyle\hat{\psi}(-\infty)=\psi_{-},\quad\hat{\psi}(+\infty)=\psi_{+}.$
(1.11)

The technical main purpose of this note is to show the following.

######  Theorem 1.

For any choice of the coefficients $\eta,\mu,\nu>0$,
(i) every Lax shock of sufficiently small amplitude possesses a dissipation
profile, and
(ii) if the coefficients satisfy the strict causality condition

$\displaystyle\mu\geq\frac{4}{3}\eta\quad\text{and}\quad\nu<\left(\frac{1}{3\eta}-\frac{1}{9\mu}\right)^{-1},$
(1.12)

then there always exist other shock waves that do not admit a dissipation
profile.

Besides the dichotomy that Theorem 1 states so distinctly, we draw the
reader’s attention to the slight difference between the two causality
conditions (1.7) and (1.12). To understand its meaning note that all
characteristic speeds $\sigma$ of the operator

$B^{\alpha\beta\gamma\delta}\displaystyle{\frac{\partial^{2}\psi_{\gamma}}{\partial
x^{\beta}\partial x^{\delta}}}$

satisfy $0\leq\sigma^{2}\leq 1$ if and only if both inequalities in (1.7)
hold, and one of the speeds is luminal, $\sigma^{2}=1$, if and only if

$\displaystyle\mu\geq\frac{4}{3}\eta\quad\text{and}\quad\nu=\left(\frac{1}{3\eta}-\frac{1}{9\mu}\right)^{-1},$
(1.13)

a condition that we therefore call _sharp causality_.

Obviously, Theorem 1 has the following two implications :

(A) In the strictly causal case (1.12), upon variation from small to large
amplitudes, necessarily some kind of transition occurs in the phase portrait
of the dynamical system (1.10).

(B) As it is well possible (though not formally proved in this paper, cf.
Remark 1 at the end of Section 3) that _all_ shocks have profiles when

$\displaystyle\nu=\nu_{*}(\eta,\mu)=\left(\frac{1}{3\eta}-\frac{1}{9\mu}\right)^{-1},$
(1.14)

it seems that for the modelling one might prefer this sharply causal tuning
over others.

Both (A) – cf. Figure 1 – and (B) are under further investigation [8].

Figure 1: Phase portrait of (1.10) for strictly causal parameter values
$(\eta,\mu,\nu)=(1,7,20)$. Left: $\tilde{q}=31/40$ . Right: $\tilde{q}=34/40$.
Gray lines indicate det $B=0$. Plots by V. Pellhammer.

Dissipation profiles for fluid dynamical shock waves have been widely studied.
For a selection of aspects, including the interesting (though quite different)
heteroclinic bifurcation in standard magnetohydrodynamics, the interested
reader is referred to [20, 13, 17, 11, 3, 10, 7].

As regards the present paper, we parametrize the ideal shock waves of the pure
radiation fluid in Section 2, Section 3 contains the proof of the theorem, and
a brief appendix concisely reconsiders the causality conditions (1.7), (1.13).

##  2 Rankine-Hugoniot conditions of pure radiation

Due to Lorentz invariance, we can restrict attention to shocks of speed $s=0$,
and because of the system’s natural isotropy, we may assume that the spatial
direction of propagation is $(1,0,0)$. Correspondingly, we fix the
spatiotemporal normal to the shock as

$(\xi_{0},\xi_{1},\xi_{2},\xi_{3})=(0,1,0,0),$

henceforth consider only states $\psi$ with $\psi_{2}=\psi_{3}=0$, and let the
indices run over 0 and 1 instead from 0 to 3; thus, for instance, the metric
$g^{\alpha\gamma}$ and the projection
$\Pi^{\alpha\gamma}=g^{\alpha\gamma}+u^{\alpha}u^{\gamma}$ on the orthogonal
complement of the velocity $u^{\alpha}$ are given by the matrices

$\left(g^{\alpha\gamma}\right)=\begin{pmatrix}-1&0\\\
0&1\end{pmatrix},\quad\left(\Pi^{\alpha\gamma}\right)=\begin{pmatrix}u_{1}^{2}&u^{0}u^{1}\\\
u^{0}u^{1}&u_{0}^{2}\end{pmatrix}.$

On the state space
$\Psi\equiv\\{\psi=(\psi^{0},\psi^{1})\in\mathbb{R}^{2}:\psi^{0}>|\psi^{1}|\\}$,
we use the temperature and the velocity,

$\theta=(-\psi_{\gamma}\psi^{\gamma})^{-1/2},\quad(u,v):=(u^{0},u^{1})=\theta(\psi^{0},\psi^{1})\
\text{with }u=(1+v^{2})^{1/2},$

as coordinates. All possible shocks can be identified by screening the
preimage set of the $\xi_{\beta}$ component of the ideal stress
$T^{\alpha\beta}$ [4, 19]:

######  Lemma 1.

On $\Psi$, the equation

$\displaystyle T^{\alpha 1}(\psi)=q^{\alpha}$ (2.1)

has more than one solution if and only if

$q^{1}>0\quad\text{and}\quad(q^{1})^{2}<(q^{0})^{2}<\frac{2}{\sqrt{3}}(q^{1})^{2}.$

In that case, it has exactly two, and these two states form a standing Lax
shock in right-moving or left-moving flow if $q^{0}>0$ or $q^{0}<0$,
respectively.

###### Proof.

As the relation (2.1)1, $T^{11}(\psi)=q^{1}$, is equivalent to

$\frac{4}{3}v^{2}+\frac{1}{3}=q^{1}\theta^{-4},$

its solution set is empty if $q^{1}\leq 0$. If $q^{1}>0$, it is the curve

$H^{1}(q^{1})=\\{h^{1}(v,q^{1}):v\in\mathbb{R}\\}$

with

$h^{1}(v,q^{1})=\frac{1}{\theta^{1}(v,q^{1})}\left(\sqrt{1+v^{2}},v\right),\quad\text{where
}\theta^{1}(v,q^{1}):=\left(\frac{1}{q^{1}}\left(\frac{4}{3}v^{2}+\frac{1}{3}\right)\right)^{-1/4}.$

As the relation (2.1)0, $T^{01}(\psi)=q^{0}$, is equivalent to

$\left(\frac{4}{3}\right)^{2}(1+v^{2})v^{2}=(q^{0})^{2}\theta^{-8}\quad\text{and}\quad
q^{0}v\geq 0,$

its solution set is $H^{0}(0)=(0,\infty)\times\\{0\\}$ if $q^{0}=0$, but
otherwise it is the curve

$H^{0}(q^{0})=\\{h^{0}(v,q^{0}):q^{0}v>0\\}$

with

$h^{0}(v,q^{0})=\frac{1}{\theta^{0}(v,q^{0})}\left(\sqrt{1+v^{2}},v\right),\quad\text{where
}\theta^{0}(v,q^{0}):=\left(\left(\frac{4}{3}\right)^{2}\frac{(1+v^{2})v^{2}}{(q^{0})^{2}}\right)^{-1/8}.$

Being the intersection of the two curves $H^{0}(q^{0})$ and $H^{1}(q^{1})$,
the solution set $H(q)$ of (2.1) can have more than one element only if we
have $q^{0}\neq 0$. As (2.1) is invariant under the reflection
$q^{0}\mapsto-q^{0},u^{1}\mapsto-u^{1}$, it is w. l. o. g. that we henceforth
assume

$\displaystyle q^{1}>0\quad\text{and}\quad q^{0}>0.$ (2.2)

Observing that

$\displaystyle h^{0}(v,q^{0})=h^{1}(v,q^{1})$
$\displaystyle\Leftrightarrow\theta^{0}(v,q^{0})=\theta^{1}(v,q^{1})$
$\displaystyle\Leftrightarrow(4v^{2}+1)^{2}=16\tilde{q}^{2}v^{2}(1+v^{2})$
$\displaystyle\Leftrightarrow(16(1-\tilde{q})v^{4}+8(1-2\tilde{q})v^{2}+1$

with

$\tilde{q}=(q^{1}/q^{0})^{2},$

we find that the equation $h^{0}(v,q^{0})=h^{1}(v,q^{1})$ has two different
positive solutions $v_{-}>v_{+}$ if and only if

$3/4<\tilde{q}<1$

and in that case, these solutions are given through

$v_{\pm}^{2}(\tilde{q})=\frac{(2\tilde{q}-1)\mp\sqrt{\tilde{q}(4\tilde{q}-3)}}{4(1-\tilde{q})}.$

We note that

$\lim_{\tilde{q}\searrow 3/4}v_{\pm}(\tilde{q})=:v_{*}=\frac{1}{\sqrt{2}}$

while

$\lim_{\tilde{q}\nearrow
1}v_{+}(\tilde{q})=\frac{1}{2\sqrt{2}}\quad\text{and}\quad\lim_{\tilde{q}\nearrow
1}v_{-}(\tilde{q})=\infty.$

I. e., the case $\tilde{q}\searrow 3/4$ corresponds to the zero amplitude and
$\tilde{q}\nearrow 1$ to the infinite amplitude limit.

As (2.1) is also invariant under homotheties
$(\psi,q)\mapsto(a\psi,a^{-2}q),a>0,$ it again causes no loss of generality
when we from now on also assume that $q^{1}=1.$ This means that we consider a
$\tilde{q}$-dependent pair of states

$(\psi_{-}(\tilde{q}),\psi_{+}(\tilde{q}))\quad\text{with
}\psi_{\pm}=\left(((4/3)v^{2}+1/3)^{1/4}((1+v^{2})^{1/2},v)\right)|_{v=v_{\pm}(\tilde{q})},\quad
3/4<\tilde{q}<1,$

that originates from the bifurcation point

$\psi_{*}=\frac{1}{\sqrt{2}}\left(\sqrt{3},1\right)$

when $\tilde{q}$ grows starting from its lower limiting value $3/4$. Since
$v_{-}>v_{+}>0$, $\psi_{-}$ as a left-hand (upstream) state and $\psi_{+}$ as
a right-hand (downstream) state form a shock that stands in a right-moving
decelerating flow.

To confirm that it is a Lax shock, we observe that up to a positive scalar
factor, the Jacobian

$\frac{\partial T^{\alpha
1}}{\partial\psi_{\gamma}}=\theta^{2}p_{\theta}(g^{\alpha
1}u^{\gamma}+g^{\alpha\gamma}u^{1}+g^{1\gamma}u^{\alpha})+(\theta^{3}p_{\theta})_{\theta}u^{\alpha}u^{1}u^{\gamma}=\frac{4}{3}\theta^{5}\left((g^{\alpha
1}u^{\gamma}+g^{\alpha\gamma}u^{1}+g^{1\gamma}u^{\alpha})+6u^{\alpha}u^{1}u^{\gamma}\right)$

is given by the matrix

$\displaystyle\begin{pmatrix}-v+6u^{2}v&u+6uv^{2}\\\
u+6uv^{2}&3v+6v^{3}\end{pmatrix}=\begin{pmatrix}v(6v^{2}+5)&\sqrt{1+v^{2}}(6v^{2}+1)\\\
\sqrt{1+v^{2}}(6v^{2}+1)&v(6v^{2}+3)\end{pmatrix},$ (2.3)

whose determinant

$\displaystyle v^{2}(6v^{2}+5)(6v^{2}+3)-(1+v^{2})(6v^{2}+1)^{2}=2v^{2}-1$
(2.4)

is positive for $v^{2}>1/2$ and negative for $v<1/2$. Thus both characteristic
speeds are positive at $\psi_{-}$, while for $\psi_{+}$ one of them is
positive and the other negative: this defines a 1-shock [16]. ∎

##  3 Shock profiles in the BDN picture

We write the profile equation

$B^{\alpha 1\gamma 1}(\psi)\psi_{\gamma}^{\prime}=T^{\alpha
1}(\psi)-q^{\alpha}.$

and a _scaled_ version of its linearization at a given state $\bar{\psi}$,

$B^{\alpha 1\gamma 1}(\bar{\psi})\psi_{\gamma}^{\prime}=\frac{\partial
T^{\alpha 1}}{\partial\psi_{\gamma}}(\bar{\psi})\psi_{\gamma},$

as

$\displaystyle B(\psi)\psi^{\prime}=F(\psi,\tilde{q})$ (3.1)

and

$\displaystyle B(\bar{\psi})\psi^{\prime}=A(\bar{\psi})\psi,$ (3.2)

respectively, with

$\displaystyle B=\tilde{\eta}\tilde{B}_{visc}-\mu B_{ther}-\nu B_{velo}$ (3.3)

where

$\displaystyle\tilde{B}_{visc}(\psi)=\begin{pmatrix}u^{2}v^{2}&-u^{3}v\\\
-u^{3}v&u^{4}\end{pmatrix},\qquad\tilde{\eta}=\frac{4}{3}\eta,\qquad$ (3.4)
$\displaystyle B_{ther}(\psi)=\begin{pmatrix}16u^{2}v^{2}&-4uv(4v^{2}+1)\\\
-4uv(4v^{2}+1)&(4v^{2}+1)^{2}\end{pmatrix},$ (3.5)
$\displaystyle\begin{aligned}
B_{velo}(\psi)&=\begin{pmatrix}(u^{2}+v^{2})^{2}&-2(u^{2}+v^{2})uv\\\
-2(u^{2}+v^{2})uv&4u^{2}v^{2}\end{pmatrix},\end{aligned}$ (3.6)

and

$A(\psi)=\begin{pmatrix}v(6u^{2}-1)&-u(6v^{2}+1)\\\
-u(6v^{2}+1)&v(6v^{2}+3)\end{pmatrix};$

according to (1.4),(1.6), and (2.3). (The minus signs in the offdiagonal
entries of the above matrices are induced by the distinction between contra-
and covariant indices, $\psi^{\epsilon}=g^{\epsilon\gamma}\psi_{\gamma}.$)

For any $\tilde{q}\in(3/4,1)$, the two states
$\psi_{-}(\tilde{q}),\psi_{+}(\tilde{q})$ we have identified above, obviously
are (the only) rest points of the ODE system (3.1).

We first consider the small-amplitude case. At the bifurcation point
$\psi_{*}$, all four matrices

$A(\psi_{*})=\frac{1}{\sqrt{2}}\begin{pmatrix}8&-4\sqrt{3}\\\
-4\sqrt{3}&6\end{pmatrix},\quad\tilde{B}_{visc}(\psi_{*})=\begin{pmatrix}u^{2}v^{2}&-u^{3}v\\\
-u^{3}v&u^{4}\end{pmatrix}=\frac{3}{4}\begin{pmatrix}1&-\sqrt{3}\\\
-\sqrt{3}&3\end{pmatrix},$
$B_{ther}(\psi)=\begin{pmatrix}16u^{2}v^{2}&-4uv(4v^{2}+1)\\\
-4uv(4v^{2}+1)&(4v^{2}+1)^{2}\end{pmatrix}=\begin{pmatrix}12&-6\sqrt{3}\\\
-6\sqrt{3}&9\end{pmatrix},$ $\displaystyle B_{velo}(\psi)$
$\displaystyle=\begin{pmatrix}(u^{2}+v^{2})^{2}&-2(u^{2}+v^{2})uv\\\
-2(u^{2}+v^{2})uv&4u^{2}v^{2}\end{pmatrix}=\begin{pmatrix}4&-2\sqrt{3}\\\
-2\sqrt{3}&3\end{pmatrix}.$

are positive multiples of orthogonal projectors, with

$\text{ ker $A(\psi_{*})=\ $ker $B_{ther}(\psi_{*})=\ $ker
$B_{velo}(\psi_{*})$ spanned by }r=\begin{pmatrix}3\\\ 2\sqrt{3}\end{pmatrix}$

and

$\text{ker }\tilde{B}_{visc}(\psi_{*})\text{ spanned by }\begin{pmatrix}3\\\
\sqrt{3}\end{pmatrix}\not\parallel r.$

As thus both $\eta\tilde{B}_{visc}(\psi_{*})$ and $\mu B_{ther}(\psi_{*})+\nu
B_{velo}(\psi_{*})$ are nontrivial multiples of orthogonal projectors with
different one-dimensional images, the matrix $B(\psi_{*})$, and so $B(\psi)$
for any $\psi$ near $\psi_{*}$, is invertible. Consequently, the augmented
profile equation

$\displaystyle\psi^{\prime}=(B(\psi))^{-1}F(\psi,\tilde{q}),\quad\tilde{q}^{\prime}=0$
(3.7)

has a 1+1-dimensional center manifold $\mathcal{C}$ at
$(\psi_{*},\tilde{q}_{*}=\frac{3}{4})$. For any value of $\tilde{q}$ above and
sufficiently close to $\tilde{q}_{*}$, the 1-dimensional fibre

$\mathcal{C}_{\tilde{q}}=\\{\psi\in\Psi:(\psi,\tilde{q})\in\mathcal{C}\\}$

contains (cf. [17], p. 242) the nearby rest points
$\psi_{-}(\tilde{q}),\psi_{+}(\tilde{q})$, and as these are the only rest
points, the segment of the curve $\mathcal{C}_{\tilde{q}}$ between them is a
single orbit $\hat{\psi}(\mathbb{R})$ which is heteroclinic to them – this
orbit is the sought after traveling wave! A straightforward center manifold
reduction (cf.[17], pp. 245, 246) confirms that it has
$\hat{\psi}(-\infty)=\psi_{-}(\tilde{q}),\hat{\psi}(+\infty)=\psi_{+}(\tilde{q})$
and not vice versa. This finishes the proof of Assertion (i).

Regarding arbitrary shocks, we first show the following properties of $B$.

######  Lemma 2.

Assume $\mu$ and $\nu$ satisfy the causality condition (1.7). Then
(i) $\psi_{-}(\tilde{q})$ is an attractor at least for certain values of
$\tilde{q}\in(3/4,1)$, unless (1.14) holds.
(ii) If (1.14) holds, $\psi_{-}(\tilde{q})$ is a hyperbolic saddle for all
values of $\tilde{q}\in(3/4,1)$.

###### Proof.

Relations (3.3)–(3.6) readily yield

$B(\psi)=\begin{pmatrix}\tilde{\eta}u^{2}v^{2}-16\mu
u^{2}v^{2}-\nu(u^{2}+v^{2})^{2}&-uv[\tilde{\eta}u^{2}-4\mu(4v^{2}+1)-2\nu(u^{2}+v^{2})]\\\
-uv[\tilde{\eta}u^{2}-4\mu(4v^{2}+1)-2\nu(u^{2}+v^{2})]&\tilde{\eta}u^{4}-\mu(4v^{2}+1)^{2}-4\nu
u^{2}v^{2}\end{pmatrix},$

which entails

$\displaystyle\det B(\psi)$
$\displaystyle=-9\tilde{\eta}\mu(1+v^{2})v^{2}-\tilde{\eta}\nu(1+v^{2})^{2}+\mu\nu(2v^{2}-1)^{2}$
$\displaystyle=c_{4}(\tilde{\eta},\mu,\nu)v^{4}+c_{2}(\tilde{\eta},\mu,\nu)v^{2}+c_{0}(\tilde{\eta},\mu,\nu)$
(3.8)

with

$\displaystyle c_{4}(\tilde{\eta},\mu,\nu)$
$\displaystyle=-9\tilde{\eta}\mu-\tilde{\eta}\nu+4\mu\nu$ $\displaystyle
c_{2}(\tilde{\eta},\mu,\nu)$
$\displaystyle={\color[rgb]{1,0,0}-9\tilde{\eta}\mu-2\tilde{\eta}\nu-4\mu\nu}$
$\displaystyle c_{0}(\tilde{\eta},\mu,\nu)$
$\displaystyle={\color[rgb]{1,0,0}-\tilde{\eta}\nu+\mu\nu}.$

Assume first that (1.14) holds. In this case,

$\displaystyle c_{4}(\tilde{\eta},\mu,\nu)$ $\displaystyle=0$ $\displaystyle
c_{2}(\tilde{\eta},\mu,\nu)$
$\displaystyle={\color[rgb]{1,0,0}-\nu(\tilde{\eta}+8\mu)}$ (3.9)
$\displaystyle c_{0}(\tilde{\eta},\mu,\nu)$
$\displaystyle={\color[rgb]{1,0,0}\nu(\mu-\tilde{\eta})}.$ (3.10)

thus $\det B(\psi^{\pm})=\nu(-(\tilde{\eta}+8\mu)v^{2}+(\mu-\tilde{\eta}))<0$
and, as according to (2.4) $\det A(\psi_{-})>0$, also

$\det(B(\psi_{-})^{-1}A(\psi_{-}))<0,$

which proves (ii).

As $c_{4}(\tilde{\eta},\mu,\nu_{*}(\eta,\mu))=0$ and $\partial
c_{4}(\tilde{\eta},\mu,\nu)/\partial\nu=4\mu-\tilde{\eta}>0$, the leading-
order coefficient in (3.8) is positive whenever

$\displaystyle\nu>\nu_{*}(\tilde{\eta},\mu).$ (3.11)

Assuming (3.11), we thus have $\det B(\psi_{-})>0$ for sufficiently large $v$,
which implies that

$\displaystyle\det B(\psi_{-})>0\text{ for any $\tilde{q}<1$ sufficiently
close to 1.}$ (3.12)

Since $B(\psi)$ has negative trace and $A(\psi_{-})>0$, this shows that the
eigenvalues of $B(\psi_{-})^{-1}$ and thus those of
$B(\psi_{-})^{-1}A(\psi_{-})$ are negative and thus

$\displaystyle\psi_{-}\text{ is an attractor for these same }\tilde{q};$
(3.13)

this verifies (i). ∎

_Proof of Assertion (ii) of Theorem 1_ : This is a direct corollary of
Assertion (i) of Lemma 2, as an attractor cannot be the $\alpha$-limit of any
other state; in particular, no heteroclinic orbit can originate in $\psi_{-}$
for any of the $\tilde{q}$ in (3.12),(3.13). Cf. the second plot in Figure 1.
$\Box$

Remark 1. Assertion (ii) of Lemma 2 leaves the possibility that in the sharply
causal case (1.13), $\psi_{-}(\tilde{q})$ is connected to
$\psi_{+}(\tilde{q})$ by a heteroclinic orbit for _all_ choices of
$\tilde{q}\in(3/4,1)$.

Remark 2. We conjecture that the phenomena established in this note have
counterparts in generalizations of the original BDN model to other fluids as
have meanhwile been considered in [15, 2, 6].

##  Appendix A Strict and sharp causality

For the reader’s convenience, we briefly discuss the causality condition (1.7)
following [1].333In their notation, (1.7) says that cofficients $\chi$ and
$\lambda$, which correspond to what here is $3\mu$ and $\nu$, must satisfy
$\chi=a_{1}\eta,\lambda\geq 3a_{1}\eta/(a_{1}-1)$ for some $a_{1}\geq 4$. See
[1], p. 10, and [6]. Fix $\tilde{\eta},\mu,\nu>0$. The dispersion relation of
$B$ being444We evaluate it in the fluid’s rest frame. This means no loss of
generality as these modes are neutral, cf. the proof of Proposition 1 in [9].

$\displaystyle 0$ $\displaystyle=\det\left(\lambda^{2}B^{\alpha 0\gamma
0}+i\xi\lambda(B^{\alpha 0\gamma 1}+B^{\alpha 1\gamma 0})-\xi^{2}B^{\alpha
1\gamma
1}\right)=\left|\begin{matrix}9\mu\lambda^{2}-\nu\xi^{2}&i\xi\lambda(3\mu+\nu)\\\
i\xi\lambda(3\mu+\nu)&\nu\lambda^{2}+(\tilde{\eta}-\mu)\xi^{2}\end{matrix}\right|,$

the speeds $\sigma=-i\lambda/\xi$ satisfy

$0=\pi(\sigma^{2})=9\mu\nu\sigma^{4}-3\mu(3\tilde{\eta}+2\nu)\sigma^{2}-\nu(\tilde{\eta}-\mu).$

Looking at the signs of the coefficients and the discriminant of the
polynomial $\pi$,

$\Delta=\tilde{\eta}\\{81\mu^{2}\tilde{\eta}+36\mu\nu(3\mu+\nu)\\},$

one sees that the speeds have $\sigma^{2}\geq 0$ if and only if
$\mu\geq\tilde{\eta}$; we assume this now.

Observing that

$\pi(1)=4\mu\nu-(9\mu+\nu)\tilde{\eta},$

one easily concludes that the speeds satisfy $\sigma^{2}\leq 1$ if and only if

$\nu\leq\left(\frac{1}{3\eta}-\frac{1}{9\mu}\right)^{-1}$

holds, with $\sigma^{2}=1$ occurring exactly in case the last inequality holds
with “=”.

## References

  * [1] F. S. Bemfica, M. M. Disconzi, J. Noronha: Causality and existence of solutions of relativistic viscous fluid dynamics with gravity, Phys. Rev. D 98 (2018), 104064.
  * [2] F. S. Bemfica, M. M. Disconzi, J. Noronha: General relativistic viscous fluid dynamics, 2020, arxiv 2009.11388.
  * [3] C. Conley, J. Smoller: On the structure of magnetohydrodynamic shock waves, Comm. Pure Appl. Math. 27 (1974), 367–375.
  * [4] C. M. Dafermos: Hyperbolic Conservation Laws in Continuum Physics. Fourth edition. Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], 325. Springer-Verlag, Berlin, 2016.
  * [5] C. Eckart: The thermodynamics of irreversible processes. 3: Relativistic theory of the simple fluid, Phys. Rev. 58 (1940), 919-924.
  * [6] H. Freistühler: A class of Hadamard well-posed five-field theories of dissipative relativistic fluid dynamics. J. Math. Phys. 61 (2020), 033101, 17 pp.
  * [7] H. Freistühler, V. Pellhammer: Dependence on the background viscosity of solutions to a prototypical nonstrictly hyperbolic system of conservation laws, SIAM J. Math. Anal. 52 (2020), 5658–5674.
  * [8] H. Freistühler, V. Pellhammer: Existence and bifurcation of shock profiles in dissipative relativistic fluid dynamics, in preparation.
  * [9] H. Freistühler, M. Reintjes, B. Temple: Decay and subluminality of modes of all wave numbers in the relativistic dynamics of viscous and heat conductive fluids, arxiv 2012.01220.
  * [10] H. Freistühler, P. Szmolyan: Existence and bifurcation of viscous profiles for all intermediate magnetohydrodynamic shock waves. SIAM J. Math. Anal. 26 (1995), 112–128.
  * [11] H. Freistühler, B. Temple: Causal dissipation and shock profiles in the relativistic fluid dynamics of pure radiation, Proc. R. Soc. A 470 (2014), 20140055.
  * [12] H. Freistühler, B. Temple: Causal dissipation in the relativistic dynamics of barotropic fluids, J. Math. Phys. 59 (2018), 063101.
  * [13] D. Gilbarg: The existence and limit behavior of the one-dimensional shock layer, Amer. J. Math. 73 (1951). 256-274.
  * [14] W. Israel: Nonstationary irreversible thermodynamics, a causal relativistic theory, Ann. Phys. 100 (1976), 310–331.
  * [15] P. Kovtun: First-order relativistic hydrodynamics is stable, J. High Energy Phys. 2019, 034, 25 pp.
  * [16] P.D. Lax: Hyperbolic systems of conservation laws II, Comm. Pure Appl. Math. 10 (1957), 537–566.
  * [17] A. Majda, R. L. Pego: Stable viscosity matrices for systems of conservation laws, J. Differ. Eqs. 56 (1985), 229 –262.
  * [18] T. Ruggeri, A. Strumia: Main field and convex covariant density for quasilinear hyperbolic systems. Relativistic fluid dynamics, Ann. Inst. H. Poincaré Sect. A (N.S.) 34 (1981), 65–84.
  * [19] A. H. Taub: Relativistic Rankine-Hugoniot equations, Phys. Rev. 74 (1948), 328–334.
  * [20] H. Weyl: Shock waves in arbitrary fluids, Comm. Pure Appl. Math. 2 (1949), 103–122.