# On the Practical Design of Tube-enhanced Multi-stage Nonlinear Model
Predictive Control

Sankaranarayanan Subramanian  Yehia Abdelsalam  Sergio Lucia  Sebastian Engell
Process Dynamics and Operations Group, TU Dortmund, Germany (e-mail:
{sankaranarayanan.subramanian, yehia.abdelsalam, sebastian.engell}@tu-
dortmund.de).  Process Automation Systems, TU Dortmund, Germany (e-mail:
<EMAIL_ADDRESS>
###### Abstract

Tube-enhanced multi-stage nonlinear model predictive control is a robust
control scheme that can handle a wide range of uncertainties with reduced
conservatism and manageable computational complexity. In this paper, we
elaborate on the flexibility of the approach from an application point of
view. We discuss the path to making design decisions to implement the novel
scheme systematically. We illustrate the critical steps in the design and
implementation of the scheme for an industrial example.

###### keywords:

Process control, model predictive control, robust control, constrained
systems, batch process

††thanks: ©2022 Sankaranarayanan Subramanian, Yehia Abdelsalam, Sergio Lucia,
Sebastian Engell. This work has been accepted to IFAC for publication under a
Creative Commons Licence CC-BY-NC-ND.

##  1 Introduction

Nonlinear Model Predictive Control (NMPC) is an optimization based advanced
control strategy that can be used to control nonlinear systems with
constraints efficiently (Findeisen and Allgöwer, 2002; Rawlings and Mayne,
2009). NMPC uses a nonlinear model to forecast the system’s future behavior
until a finite horizon and obtains a sequence of control inputs that optimizes
a given objective and satisfies the state and input constraints. The control
input obtained at the first prediction step is applied to the system, and the
optimization problem is solved at the next step in a receding horizon fashion.
The performance of the NMPC scheme depends critically on the accuracy of the
model. The uncertainties present in the model may lead to erroneous forecasts
resulting in constraint violation, loss of recursive feasibility, or even
closed-loop instability. Hence, to systematically address the presence of
uncertainties in the model, robust NMPC schemes are necessary.

Open-loop min-max model predictive control (MPC) is one of the earliest
proposed robust strategies in the context of linear systems (Campo and Morari,
1987). The scheme predicts a sequence of control inputs that minimize the
worst-case objective while satisfying the constraints for all realizations of
the uncertainty. However, predicting one control input per prediction step may
lead to conservative behavior. The closed-loop min-max MPC (Scokaert and
Mayne, 1998) and the multi-stage MPC formulations (Muñoz de la Peña et al.,
2005; Lucia et al., 2013, 2020) were proposed to address this conservatism by
adapting the control inputs in the predictions for different possible
realizations of the uncertainties. This models recourse in the problem
formulation and achieves less conservative results. Because of the reduced
conservatism, the multi-stage approach has been applied in different case
studies (Puschke and Mitsos, 2018; Jang et al., 2016; Krishnamoorthy et al.,
2018).

Unlike the min-max approaches, multi-stage NMPC optimizes the weighted mean of
all the scenarios considered in the prediction while accounting for the
feedback in the predictions. The drawback of the scheme is that the complexity
of the optimization problem grows rapidly when the number of uncertainties and
the prediction horizon increase. Tube-based NMPC is another class of robust
NMPC approach that provides computationally simpler alternative that generates
nominal trajectories optimally and then regulates the processes to converge to
or stay close to these (see (Mayne et al., 2011; Rubagotti et al., 2011;
Villanueva et al., 2017)). In (Mayne et al., 2011), a nominal and an ancillary
controller are employed hierarchically to achieve robustness. In the presence
of significant uncertainties, tube-based NMPC can be conservative while the
complexity remains close to the nominal NMPC controller.

To improve the trade-off between optimality and complexity compared to the
existing robust MPC strategies, Tube-enhanced multi-stage (TEMS) MPC was
proposed for the linear case in (Subramanian et al., 2021) and for the
nonlinear case in (Subramanian et al., 2018). The sufficient conditions to
achieve closed-loop stability using the TEMS NMPC scheme were derived in
(Subramanian et al., 2022). While it is possible to obtain prior safety and
stability guarantees of the closed-loop, it is however very difficult to
verify in practice whether the required assumptions are satisfied.

The TEMS NMPC approach, in its basic form, employs two hierarchical robust
multi-stage NMPC controllers. In this contribution, we elaborate on the
flexibility and the applicability of the TEMS NMPC approach from an
application point of view. The following aspects of the design process are
explained:

  1. 1.

Construction of the scenario tree for handling different kinds of
uncertainties with varying significance.

  2. 2.

Different formulations of the ancillary controller to achieve robust
constraint satisfaction.

  3. 3.

Constraint tightening of the primary controller.

  4. 4.

Initialization of the scenario tree for the primary controller.

This is followed by the application of the method to an industrially relevant
fed-batch polymerization process.

##  2 System description and control goals

The dynamics of the system is described as follows:

$\displaystyle x^{+}=f(x,u,d),$ (1)

where
$x\in\mathbb{R}^{n_{x}},u\in\mathbb{R}^{n_{u}},\,d\in\mathbb{D}\subset\mathbb{R}^{n_{d}}$
denote the states, inputs and model uncertainties (both parametric and
additive uncertainties). The nonlinear mapping
$f:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{d}}\to\mathbb{R}^{n_{x}}$
denotes the system dynamics. The uncertainty bound $\mathbb{D}$ is assumed to
be compact and can be represented as $\mathbb{D}=\\{d\mid d^{L}\leq d\leq
d^{U}\\}$, where $d^{L}$ and $d^{U}$ denote the lower and the upper bounds of
the uncertainty. The state constraints are denoted by $\mathbb{X}$, and the
input constraints are denoted by $\mathbb{U}$. The constraint set $\mathbb{X}$
is assumed to be closed and $\mathbb{U}$ is assumed to be compact. The control
goal can be either to track a setpoint or to economically optimize the
performance of the plant while satisfying the state and input constraints for
all realizations of the uncertainties at all times. In addition, the
controller should not be overly conservative and should be real-time
implementable. Here, we assume for simplicity that full-state information is
available.

##  3 Tube-enhanced multi-stage NMPC

Tube-enhanced multi-stage NMPC employs two controllers hierarchically to
achieve the control goals: a primary multi-stage NMPC controller and an
ancillary multi-stage NMPC controller. The primary controller achieves
robustness of the primary system, for which the dynamics are described by

$\displaystyle z^{+}=f(z,v,\overline{d})$ (2)

where $z\in\mathbb{R}^{n_{x}}$ denotes the state, $v\in\mathbb{R}^{n_{u}}$
denotes the control input, and
$\overline{d}\in\overline{\mathbb{D}}\subset\mathbb{R}^{n_{{d}}}$ denotes the
uncertainties of the primary system. The nonlinear mapping
$f:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{d}}\to\mathbb{R}^{n_{x}}$
is the same as in (1) for all $\overline{d}\in\overline{\mathbb{D}}$. The
uncertainty set $\overline{\mathbb{D}}$ is a finite set and is generated by
sampling the true uncertainty set $\mathbb{D}$. The predicted trajectories
generated by the primary controller are then tracked by the ancillary
controller to mitigate the uncertainties that are not considered in the
primary system.

###  3.1 Primary controller

To achieve robustness against the uncertainties that are considered in the
primary system described in (2), multi-stage NMPC is employed. The
optimization problem formulation for the primary controller is given below:

$\displaystyle\min_{v_{k}^{j}\forall(j,k)\in
I_{[0,N-1]}}\;\;\;\sum_{k=0}^{N-1}\sum_{j=1}^{s^{k}}\omega_{k}^{j}\ell(z_{k}^{j},v_{k}^{j})+\sum_{j=1}^{s^{N}}V_{f}(z_{N}^{j}),$
(3a) subject to: $\displaystyle
z_{k+1}^{c(j,r)}=f(z_{k}^{j},v_{k}^{j},d^{r}),$
$\displaystyle\,\forall\,(j,k)\in I_{[0,N-1]},\,r\in\\{1,\,\dots,\,s\\},$ (3b)
$\displaystyle\qquad z_{k}^{j}\in\mathbb{Z},v_{k}^{j}\in\mathbb{V},$
$\displaystyle\,\forall\;(j,k)\in I_{[0,N-1]},$ (3c) $\displaystyle\qquad
z_{N}^{j}\in\mathbb{Z}_{f},$ $\displaystyle\,\forall\;(j,N)\in I_{N},$ (3d)
$\displaystyle v_{k}^{j}=v_{k}^{l}\;\text{if }z_{k}^{j}=z_{k}^{l}$
$\displaystyle\,\forall\;(j,k),(l,k)\in I_{[0,N-1]},$ (3e)

where $N$ denotes the prediction horizon, $s$ denotes the number of branches
in the scenario tree, $I$ denotes the set of indices $(j,k)$ in the scenario
tree, $I_{[k_{1},k_{2}]}$ denotes the set of indices from the time step
$k_{1}$ until the step $k_{2}$, $I_{N}$ denotes the indices at the terminal
nodes, $\ell(z_{k}^{j},v_{k}^{j})$ denotes the stage cost, and
$V_{f}(z_{N}^{j})$ denotes the terminal cost. The weights associated with each
node are denoted by $\omega_{k}^{j},\forall(j,k)\in I$. (3b) denotes the model
of the system where the predicted state $z_{k+1}^{c(j,r)}$ is the child node
obtained from its parent state $z_{k}^{j}$ for the input $v_{k}^{j}$ and the
realization of the uncertainty $d_{k}^{r(j)}\in\overline{\mathbb{D}}$. The
state and input constraints are defined in (3c). The terminal constraints are
defined in (3d). The state constraint set is denoted by $\mathbb{Z}$, the
input constraint set is denoted by $\mathbb{V}$, and the terminal constraint
set is denoted by $\mathbb{Z}_{f}$. Both $\mathbb{Z}$ and $\mathbb{V}$ are
tightened compared to original constraint sets, and
$\mathbb{Z}_{f}\subseteq\mathbb{Z}\subseteq\mathbb{X}$ and
$\mathbb{V}\subseteq\mathbb{U}$ hold. The non-anticipativity constraints are
defined in (3e). The realized uncertainty of the primary system
$\overline{d}(t)$ for all $t\geq 0$ can be estimated by solving the
optimization problem:

$\displaystyle\overline{d}(t){\in}\arg\min_{\overline{d}\in\mathbb{\overline{D}}}|x(t+1)-f(x(t),u(t),\overline{d})|,$
(4)

where $x(t)$ and $x(t+1)$ denote the state measurements at time steps $t$ and
$t+1$. The resulting $\overline{d}(t)$ is used to obtain $z(t+1)$ using (2),
where $\overline{d}=\overline{d}(t)$.

###  3.2 Ancillary controller

Robustness against the uncertainties that are not explicitly part of the
primary system (2) is achieved by the ancillary controller. The ancillary
controller tracks the predictions of the primary controller, and the
optimization problem formulation is

$\displaystyle\min_{u_{k}^{j}\forall(j,k)\in
I_{[0,N-1]}}\;\;\;\sum_{k=0}^{N-1}\sum_{j=1}^{s^{k}}\omega_{a,k}^{j}\ell_{a}(x_{k}^{j}-{z}_{k}^{j*},u_{k}^{j}-{v}_{k}^{j*})+\sum_{j=1}^{s^{N}}V_{fa}(x_{N}^{j}),$
(5a) subject to: $\displaystyle
x_{k+1}^{c(j,r)}=f(x_{k}^{j},u_{k}^{j},d^{r}),\,$
$\displaystyle\forall(j,k)\in I_{[0,N-1]},\,r\in\\{1,\dots,s\\},$ (5b)
$\displaystyle u_{k}^{j}\in\mathbb{U},\,$ $\displaystyle\forall(j,k)\in
I_{[0,N-1]},$ (5c) $\displaystyle u_{k}^{j}=u_{k}^{l}\;\text{if
}x_{k}^{j}=x_{k}^{l},\,$ $\displaystyle\forall(j,k),(l,k)\in I_{[0,N-1]}.$
(5d)

where the stage cost of the ancillary controller is defined as
$\ell_{a}(x_{k}^{j}-{z}_{k}^{j*},u_{k}^{j}-{v}_{k}^{j*})$ and is tuned to
track the optimal state and input trajectories that were obtained from the
primary controller ${z}_{k}^{j*}$ and ${v}_{k}^{j*}$ for all $(j,k)\in
I_{[0,N-1]}$. The terminal cost is denoted as $V_{fa}(x_{N}^{j})$. The weights
associated with the nodes of the scenario tree of the ancillary controller are
denoted by $\omega_{a,k}^{j}$ for all $(j,k)\in I$. The state constraints are
not part of the ancillary controller problem. This enables the recursive
feasibility of the ancillary controller by construction. The robust constraint
satisfaction is achieved by tuning the stage cost of the ancillary controller.
At the initial time step $t=0$, the primary controller is initialized with
$x(0)$ and $u=v_{0}^{1*}$ is applied to the system. For all $t\geq 1$, the
implementation steps are summarized in Algorithm 1.

Algorithm 1 Tube-enhanced multi-stage NMPC algorithm Require |  All the ingredients of (3) and (5).
---|---
Step 1 |  Measure the state $x(t)$.
Step 2 |  Solve the optimization problem (4) and obtain $\overline{d}(t-1)$ and simulate (2) to obtain $z(t)$.
Step 3 |  Solve (3) using $z_{0}^{1}=z(t)$ to get $\\{{z}_{k}^{j*}\\},\,\forall(j,k)\in I_{[0,N]}$ and $\\{{v}_{k}^{j*}\\},\,\forall(j,k)\in I_{[0,N-1]}$.
Step 4 |  Initialize the ancillary controller with $x(t)$, pass the reference trajectories $\\{{z}_{k}^{j*}\\},\,\forall(j,k)\in I_{[0,N]}$ and $\\{{v}_{k}^{j*}\\},\,\forall(j,k)\in I_{[0,N-1]}$ and solve (5).
Step 5 |  Apply $u_{0}^{1*}$ as the input to the plant, set $t=t+1$ and go to step 1 at the next time step.

The closed-loop dynamics of the composite system is given by

$\displaystyle x^{+}$ $\displaystyle=f(x,\kappa_{N}(x,z),d),\,d\in\mathbb{D},$
(6a) $\displaystyle z^{+}$
$\displaystyle=f(z,\kappa^{p}_{N}(z),\overline{d}),\,\overline{d}\in\mathbb{\overline{D}},$
(6b)

where $\kappa_{N}(x,z)=u_{0}^{1*}$ denotes the control law of the ancillary
controller obtained by solving the optimization problem (5),
$\kappa_{N}^{p}(z,t)=\overline{v}_{0}^{1*}$ denotes the control law for the
primary system (2), $\overline{d}$ is obtained by solving (4), and
$d\in\mathbb{D}$ denotes the realized uncertainty in the plant.

The differences among the robust NMPC schemes are illustrated in Figure 1.
Multi-stage NMPC (Lucia et al., 2013) typically requires more branches to
handle the same set of uncertainties compared to TEMS NMPC and employs a
single controller which is initialized with the true state of the system.
Tube-based NMPC (Mayne et al., 2011) uses two controllers as in the TEMS NMPC
scheme but employs only a single scenario. This can lead to stringent
constraint tightening in the presence of significant uncertainties and
therefore negatively affect the performance. Considering only the significant
uncertainties in the TEMS NMPC controllers results in an improved trade-off
between optimality and complexity.

Figure 1: Illustration of multi-stage NMPC, tube-based NMPC, and TEMS NMPC.

##  4  Design details

Before discussing the design steps, we first discuss ways to simplify the
controllers.

###  4.1 Controller simplifications

The scenario tree of the primary and ancillary controller problems grows
exponentially with respect to the prediction horizon. To avoid the growth in
problem complexity, the uncertainty realizations can be assumed to remain
constant or a disturbance model can be employed after a certain prediction
step called robust horizon $N_{R}$. Beyond $N_{R}$, there exists only one
child per parent node. Though there exists no prior guarantees, robust
constraint satisfaction and recursive feasibility are often achieved in
practice even with $N_{R}=1$ as demonstrated in (Lucia et al., 2013, 2014).

The second possible simplification is the choice of the ancillary controller.
Instead of tracking all the state and input trajectories of the primary
controller, only a subset of the trajectories that are critical from control
perspective can be tracked. In the simplest possible case, a single trajectory
can be tracked. In this case, the ancillary control problem simplifies as
follows:

$\displaystyle\min_{{u_{k}},\forall
k\in\\{0,\dots,N-1\\}}\;\;\;\sum_{k=0}^{N-1}\ell_{a}(x_{k}-{z}_{k}^{*},u_{k}-{v}_{k}^{*})+V_{fa}(x_{N}),$
(7a) subject to: $\displaystyle x_{k+1}=f(x_{k},u_{k},d_{\mathrm{nom}}),$
$\displaystyle\quad\forall\,k\in\\{0,\dots,N-1\\},$ (7b) $\displaystyle
u_{k}^{j}\in\mathbb{U},$ $\displaystyle\quad\forall\,k\in\\{0,\dots,N-1\\},$
(7c)

where only the trajectories corresponding to the nominal uncertainty
$d_{\mathrm{nom}}$ are tracked. This reduces the complexity of the online
implementation of the ancillary controller. Instead of tracking a nominal
trajectory, one can even adapt the trajectories to be tracked depending on the
realization of the uncertainty (Abdelsalam et al., 2020).

###  4.2  Construction of the controllers

The scenario tree is constructed for all the realizations of the uncertainties
in the set $\overline{\mathbb{D}}$. The choice of elements of the set
$\overline{\mathbb{D}}$ (defined in (2)) is the first design decision to make.
A heuristic for scenario tree construction is to consider extreme values of
the uncertainties along with the nominal values of the parameter. In Lucia et
al. (2014), it was shown that robust constraint satisfaction was achieved for
all $d\in\mathbb{D}$ even for nonlinear systems by considering only the
extreme realizations and nominal values of the uncertainties for industrially
relevant case studies. However, the problem complexity increases rapidly with
the dimension of the uncertainty $n_{d}$.

In TEMS NMPC, the uncertainties shall be ordered in terms of significance
(e.g., based on sensitivity analysis) and only the most significant
uncertainties should be used for the construction of the scenario trees. Let
the number of significant uncertainties be denoted as $n_{\overline{d}}$ that
is smaller than $n_{d}$. Let the uncertainty $d$ be defined as
$d\triangleq\\{d_{1},d_{2},\dots,d_{n_{d}}\\}$, the elements of which are
ordered such that the first $n_{\overline{d}}$ elements are significant. Let
$\overline{\mathbb{D}}_{c}$ be defined as

$\displaystyle\overline{\mathbb{D}}_{c}=\bigg{\\{}d\mid$ $\displaystyle
d_{i}^{L}\leq d_{i}\leq d_{i}^{U},\forall i\in\\{1,\dots,n_{\overline{d}}\\},$
$\displaystyle d_{i}=d_{\mathrm{nom}},\forall
i\in\\{n_{\overline{d}}+1,\dots,n_{d}\\}\bigg{\\}}.$ (8)

Then $\overline{\mathbb{D}}$ can be constructed by sampling the extreme and
nominal values of all $d_{i},\,i\in\\{1,\dots,n_{\overline{d}}\\}$. Various
works in the literature concerning scenario tree generation can also be
employed in the context of TEMS NMPC: quadrature-based approximation
(Leidereiter et al., 2014), vector quantization approach (Goodwin et al.,
2009) or sensitivity-based construction (Thombre et al., 2021). The scenario
tree of the ancillary has the same structure used in the primary controller.
However, there exist computationally simpler alternatives as discussed in
Section 4.1.

###  4.3 Constraint tightening of the primary controller and tuning the
ancillary controller

The constraints of the primary controller should be tightened to provide room
for the ancillary controller to satisfy the constraints of the original
system. There is an interplay between the stage costs for the ancillary
controller and the constraint tightening of the primary controller. Hence
during the simulation study, it is recommended to study the closed-loop
response as the first step. The stage cost of the ancillary controller can be
chosen as a quadratic function as follows:

$\displaystyle\ell_{a}(x-z,u-v)=(x-z)^{T}Q(x-z)+(u-v)^{T}R(u-v),$ (9)

where $Q$ and $R$ are positive semi-definite matrices. For simpler
implementations, $Q$ and $R$ can be chosen as diagonal matrices with non-
negative elements. The gains can be chosen to track the critical states and
inputs. If the tracking of the ancillary controller is satisfactory, the
primary controller constraints can be tightened to account for the magnitude
of the constraint violations.

A simple way to tighten the constraints is by simulation studies as proposed
in Mayne et al. (2011). Let the state and input constraints be represented as
$g_{i}(x,u)\leq 0$, $i=\\{1,2,\cdots,n_{c}\\}$, where $n_{c}$ denotes the
number of constraints. The constraints of the primary controller can be
represented as $g_{i}(x,u)\leq-\delta_{i}$, $i=\\{1,2,\dots,n_{c}\\}$. The
parameters $\delta_{i}\geq 0,\,\forall i=\\{1,2,\dots,n_{c}\\}$ associated
with each constraint can be chosen based on the maximum violations observed
for the final tuning of the ancillary controller. Once the constraints are
tightened, rigorous constraint satisfaction can be verified by performing
extensive simulation studies for different values of uncertainties and for
different initial conditions.

###  4.4 Reinitialization and solution of the controllers

Both (3) and (5) are solved recursively. At the initial step $(t=0)$, the
primary controller is initialized with the true state (or the initial
estimate) of the system. For all $t>0$, the primary controller is always
initialized with one of its predicted states. This is because (4) finds the
closest uncertainty realization considered in the scenario-tree branches to
the measured state. The optimization problem (3) can be solved in previous
time steps for all uncertainty realizations considered in the scenario tree.
As soon as the measurement is received, the trajectories corresponding to the
realized uncertainty can be passed to the ancillary controller. In some cases,
it is possible that the optimization problem (4) results in fluctuations of
the uncertainty because of various disturbances. If the disturbance dynamics
are available, they can be employed to improve the accuracy of the resulting
estimates of $\overline{d}$. An alternative is to include regularization term
for all $t\geq 1$ in (4) as follows:

$\displaystyle\overline{d}(t){\in}\arg\min_{\overline{d}\in\mathbb{\overline{D}}}|x(t+1)-f(x(t),u(t),\overline{d})|+|W(\overline{d}(t-1)-\overline{d})|,$

were $W$ is a weighting matrix. The weighting matrix can be chosen as a
diagonal matrix with large values for time-invariant or slowly time varying
parameters and small values for rapidly time-varying parameters. To improve
the performance of the closed-loop, the parameters of the primary system can
be estimated for all $t>0$ as follows:

$\displaystyle\overline{d}(t){\in}\arg\min_{\overline{d}\in\mathbb{\overline{D}}_{c}}|x(t+1)-f(x(t),u(t),\overline{d})|+|W(\overline{d}(t-1)-\overline{d})|,$
(10)

where $\mathbb{\overline{D}}_{c}$ contains the closed intervals of the
parameters considered in the scenario tree as defined in (4.2). The primary
state at the next time step ($z(t+1)$) can be computed from the current values
of the state and input of the primary controller and the resulting uncertainty
of the problem (10), using (2). If (10) is employed, the primary dynamics (2)
are defined for all $\overline{d}\in\overline{\mathbb{D}}_{c}$. When the
initialization based on (10) is employed, it is assumed that the primary
controller can handle all uncertainties
$\overline{d}\in\mathbb{\overline{D}}_{c}$. The ancillary controller is always
initialized with the true state of the system and the solution is computed
online.

###  4.5 Tracking vs. economic objectives

For tracking objectives, a guaranteed stabilizing control law using TEMS NMPC
is possible. See (Subramanian et al., 2022) for the necessary assumptions for
prior guarantees. If the system under investigation is too complex to verify
the assumptions, the design steps provided here can help achieve required
robust constraint satisfaction and closed-loop stability. Offset-free tracking
requires additional steps which are beyond the scope of the paper. The scheme
implemented based on the discussed methods may lead to offsets at steady-
state.

In the case of economic objectives, the system often operates at the system
boundary. Employing multi-stage NMPC, ignoring minor disturbances, may lead to
oscillatory input moves (as a response to minor state constraint violations).
This oscillation issue is resolved in the TEMS NMPC because the primary
controller is affected only by the uncertainties considered in the tree and
the ancillary controller does not have state constraints. This avoids
aggressive moves from the ancillary controller if tuned properly. At the same
time, implementing tube-based NMPC in the presence of significant
uncertainties may lead to conservatism. TEMS NMPC can improve the economic
performance because of a less conservative primary controller.

###  4.6 Further extensions

The TEMS NMPC scheme can be easily extended to the output feedback case as
well. See (Subramanian et al., 2018) for details on the same. One can use any
efficient estimator to estimate the states of the system and initialize the
ancillary controller using the state estimates. The primary controller is
implemented the same way as in the full-state feedback case with a more
stringent constraint tightening to account also for the estimation errors.

In the case of time-invariant parameters, suitable parameter estimation
algorithms can be used to adapt the scenario tree online. The adaptive TEMS
NMPC scheme can yield major performance advantages as shown in (Abdelsalam et
al., 2021).

##  5 Example

We present the details on the design of TEMS NMPC for the example of a
polymerization process that was provided by BASF SE (Lucia et al., 2014).
There are mass balances of water $m_{\mathrm{W}}$, monomer $m_{\mathrm{A}}$
and polymer $m_{\mathrm{P}}$. There are five energy balances providing the
dynamics of temperatures of reactor ${T}_{\text{R}}$, vessel ${T}_{\text{S}}$,
jacket ${T}_{\text{M}}$, heat-exchanger mixture ${T}_{\text{EK}}$ and coolant
at the outlet of the heat exchanger ${T}_{\text{AWT}}$. The reaction is
exothermic, and the jacket and the heat exchanger are employed to regulate the
temperature of the reactor. The parameter values and the further details of
the model can be found in (Lucia et al., 2014). The constraints on the reactor
temperature $T_{\text{R}}$ help in achieving the desired quality of the end
product, and the adiabatic safety constraint on ${T}_{\text{ad}}$ ensures safe
operation of the plant. The manipulated variables are feed rate of the monomer
$\dot{m}_{\text{F}}$, the jacket inlet temperature $T_{\text{M}}^{\text{IN}}$,
and the inlet coolant temperature of the heat exchanger
$T_{\text{AWT}}^{\text{IN}}$. The bounds on the state and the control inputs
are given in Table 1.

The reaction rate $k_{0}$ and the enthalpy of the reaction $\Delta
H_{\text{R}}$ are uncertain by $\pm 30\%$ relative to the nominal values. The
nominal values of $\Delta H_{\text{R}}$ and $k_{0}$ are $-950$kJ kg-1 and $7$,
respectively. The bounds are given by $\Delta H_{\text{R}}\in[-665,-1235]$kJ
kg-1 and $k_{0}\in[4.9,9.7]$. In addition, the dynamics of all the states are
affected by additive disturbances whose bounds at each discrete time-steps are
given by $\pm 0.5$ kg/h for the mass balance of water, $\pm 5$ kg/h for the
mass balances of monomer and polymer, and $\pm 0.1^{\circ}$C for the energy
balances at discrete time steps. The magnitudes of the uncertainties are
different when compared to the simulation study presented in (Subramanian et
al., 2018).

Table 1. State and input constraints  State/ input Min. Max. Primary
controller bounds Unit $T_{\text{R}}$ 88.0 92.0 [88.3, 91.7]
${}^{\circ}\text{C}$ $T_{\text{ad}}$ 0 109 [1, 108] ${}^{\circ}\text{C}$
$\dot{m}_{\text{F}}$ 0 30000 [0, 29990] $\frac{\text{kg}}{\text{h}}$
$T_{\text{M}}^{\text{IN}}$ 60 100 [61, 99] ${}^{\circ}\text{C}$
$T_{\text{AWT}}^{\text{IN}}$ 60 100 [61, 99] ${}^{\circ}\text{C}$

###  5.1 Scenario tree construction

For the construction of the scenario tree, only the parametric uncertainties
were considered because they are the most significant ones. Using the
heuristics discussed in Section 4.2, we built the scenario for the extreme and
nominal realizations of $\Delta H_{\text{R}}$ and $k_{0}$, where $\Delta
H_{\text{R}}\in\\{-950,-655,-1235\\}$ and $k_{0}\in\\{7,4.9,9.1\\}$. It can be
seen that $s=3^{2}=9$. For the additive uncertainties, the nominal values are
considered. The approach was implemented for a robust horizon $N_{R}=1$. After
the first prediction step, all the uncertainties are assumed to be constant.
This results in nine scenarios for both the primary controller and the
ancillary controller. If we do not classify the uncertainty and employ multi-
stage NMPC, it would require $3^{10}=59,049$ scenarios (including the additive
uncertainties) for $N_{R}=1$. The advantage of TEMS in the complexity
reduction is obvious in comparison to multi-stage NMPC.

###  5.2 Cost function of the primary controller

The goal is to produce a required amount of polymer as fast as possible. To
achieve this goal, the following stage cost was chosen as in (Lucia et al.,
2014) for the primary controller.

$\displaystyle\ell(z_{k}^{j},u_{k}^{j})$
$\displaystyle=-{m}_{\text{P},k}^{j}+r_{1}(\Delta\dot{m}_{\text{F},k}^{j})^{2}+r_{2}(\Delta
T_{\text{M},k}^{\text{IN},j})^{2}+r_{3}(\Delta
T_{\text{AWT},k}^{\text{IN},j})^{2},$

where $r_{1}$, $r_{2}$ and $r_{3}$ are tuning parameters that penalize the
control moves. The values of the tuning parameters were chosen as
$r_{1}=0.125$, $r_{2}=4$ and $r_{3}=0.25$. The terminal cost was chosen as
$V_{f}(z)=0,\,\forall z\in\mathbb{R}^{n_{x}}$ and the terminal set is chosen
as $\mathbb{Z}_{f}=\mathbb{Z}$. The sampling time was chosen as fifty seconds
and the prediction horizon of $N=20$ was chosen. Extensive simulation studies
were performed for all uncertainties in $\overline{\mathbb{D}}$ without
constraint tightening.

###  5.3 Cost function of the ancillary controller

The stage cost of the ancillary controller is chosen as a quadratic function
defined in (9). The input weight matrix was fixed to identity matrix
$R=\mathtt{I}^{3\times 3}$. The tuning of $Q$ was done considering two
aspects: Good tracking of $m_{\text{P}}$ achieves the required goal and
tracking of $T_{\text{R}}$ is required to guarantee robust constraint
satisfaction. The constraint $T_{\text{R}}$ remains active throughout the
batch and hence, poses a challenge in tuning of $Q$. After several tries, $Q$
was chosen as $\text{diag}(0,0,1,500,0,0,0,0)$. The gain of $500$ was required
because of the differences in the absolute values of $m_{\text{P}}$ and
$T_{\text{R}}$. If smaller gains were chosen for tracking $T_{R}$, the
constraint violations were relatively large at the end of the batch. The
terminal cost of the ancillary controller was chosen as $V_{fa}(x)=0,\,\forall
x\in\mathbb{R}^{n_{x}}$.

###  5.4 Constraint tightening of the primary controller

After choosing the ancillary stage cost, simulation studies were performed
mainly for the extreme realizations of the uncertainties and the worst-case
constraint violations were recorded. The primary controller constraints were
tightened based on the procedure given in Section 4.3. The resulting
constraint tightening for the primary controller can be found in Table 1.

###  5.5 Computation of the solution

Orthogonal collocation on finite elements was employed for the discretization
of the nonlinear dynamics, and CasADi (Andersson et al., 2012) was employed
for the automatic generation of first and second-order exact derivatives. The
resulting nonlinear programming problem was solved using IPOPT (Wächter and
Biegler, 2006). For the primary controller, the re-initialization strategy
(10) was employed.

###  5.6 Results and comparison

After all the aspects of the controller had been tuned, a uniform grid of
$100$ realizations was chosen in the uncertain parametric space. Each
realization was selected for a batch run and $100$ batch runs were simulated.
The additive disturbances were varied randomly at every time step during the
batch. Selected closed-loop trajectories obtained by applying TEMS NMPC scheme
for the industrial polymerization process are illustrated in Figure 2. It can
be seen that both the quality and safety constraints are satisfied at all
times. Also no constraint violations were noted for all $100$ batch runs.

The method is compared with the tube-based NMPC and the multi-stage NMPC
approaches. As explained above, considering all the uncertainties in the
multi-stage NMPC controller leads to a high computational complexity.
Completely ignoring the additive uncertainties resulted in feasibility issues
because of violations of the constraint on $T_{\text{R}}$. The additive
disturbance on the dynamics of $T_{\text{R}}$ was therefore considered
together with the parametric uncertainties in the scenario tree. For the
implementation of multi-stage NMPC, three uncertainties were considered
leading to $27$ scenarios are employed. It was very difficult to tune the
tube-based NMPC because the nominal trajectory varied considerably compared to
the uncertain trajectory. The tube-based NMPC was implemented with very
stringent constraint tightening. Even then, robust constraint satisfaction was
not possible for all the batches.

Figure 2: Trajectories of reactor temperature $T_{\text{R}}$, adiabatic safety
temperature $T_{\text{ad}}$, feed rate $\dot{m}_{\text{F}}$, jacket inlet
temperature $T_{\text{M}}^{\text{IN}}$ for various simulation runs of the TEMS
NMPC scheme.

Table 2. Performance comparison between multi-stage NMPC, tube-based NMPC
(Tube) and TEMS NMPC scheme concerning batch times, constraint satisfaction
and computation times

Parameters |  MS |  Tube |  TEMS
---|---|---|---
Average batch time [h] |  2.1 |  3.1 |  1.96
No. of violation of $T_{\text{R}}$ |  0 |  0 |  0
No. of. violation of $T_{\text{ad}}$ |  0 |  14 |  0
Avg. comp. time per iter. [s] |  3.07 |  0.07 |  0.98

The results are given in Table 2. It can be seen that the TEMS NMPC scheme
yields the best average batch time. The average computational complexity was
less than one second. It can also be seen that no constraints were violated.
Multi-stage NMPC had a higher computational cost. However, it does not result
in an improved performance. This is because there were small oscillations in
the control moves due to ignored disturbances. Tube-based NMPC on the average
leads to 1 hour longer batch times when compared to both multi-stage NMPC and
TEMS NMPC. Tube-based NMPC requires the smallest computation time, but it also
resulted in violations of safety constraints in $14$ batches. Overall, it can
be concluded that the TEMS NMPC scheme provides a better trade-off between
optimality and complexity compared to the other two approaches.

##  6 Conclusion

Tube-enhanced multi-stage NMPC offers flexibility to achieve the desired
trade-off between optimality and complexity. We presented various design
details, simplifications, improvements and extensions of the tube-enhanced
multi-stage NMPC approach. The scheme was applied to an industrially relevant
example and all the design steps were discussed. The resulting scheme showed
superior performance when compared to both the multi-stage NMPC and the tube-
based NMPC schemes both in terms of performance and computational complexity.

## References

  * Abdelsalam et al. (2021) Abdelsalam, Y., Subramanian, S., Aboelnour, M., and Engell, S. (2021).  Adaptive tube-enhanced multi-stage nonlinear model predictive control.  _IFAC-PapersOnLine_ , 54(3), 212–218.
  * Abdelsalam et al. (2020) Abdelsalam, Y., Subramanian, S., and Engell, S. (2020).  A Simplified Implementation of Tube-Enhanced Multi-Stage NMPC.  _IFAC-PapersOnLine_ , 53(2), 11563–11569.  21st IFAC World Congress.
  * Andersson et al. (2012) Andersson, J., Åkesson, J., and Diehl, M. (2012).  Dynamic optimization with CasADi.  In _Proc. of the 51st IEEE Conference on Decision and Control_ , 681–686.
  * Campo and Morari (1987) Campo, P. and Morari, M. (1987).  Robust model predictive control.  In _Proc. of the American Control Conference_ , 1021–1026.
  * Findeisen and Allgöwer (2002) Findeisen, R. and Allgöwer, F. (2002).  An introduction to nonlinear model predictive control.  In _21st Benelux meeting on systems and control_ , volume 11, 119–141. Citeseer.
  * Goodwin et al. (2009) Goodwin, G.C., Østergaard, J., Quevedo, D.E., and Feuer, A. (2009).  A vector quantization approach to scenario generation for stochastic nmpc.  In _Nonlinear Model Predictive Control_ , 235–248. Springer.
  * Jang et al. (2016) Jang, H., Lee, J.H., and Biegler, L.T. (2016).  A robust nmpc scheme for semi-batch polymerization reactors.  _IFAC-PapersOnLine_ , 49(7), 37–42.
  * Krishnamoorthy et al. (2018) Krishnamoorthy, D., Suwartadi, E., Foss, B., Skogestad, S., and Jäschke, J. (2018).  Improving scenario decomposition for multistage mpc using a sensitivity-based path-following algorithm.  _IEEE Control Systems Letters_ , 2(4), 581–586.
  * Leidereiter et al. (2014) Leidereiter, C., Potschka, A., and Bock, H.G. (2014).  Quadrature-based scenario tree generation for nonlinear model predictive control.  _IFAC Proceedings Volumes_ , 47(3), 11087–11092.
  * Lucia et al. (2014) Lucia, S., Andersson, J., Brandt, H., Diehl, M., and Engell, S. (2014).  Handling uncertainty in economic nonlinear model predictive control: a comparative case-study.  _Journal of Process Control_ , 24, 1247–1259.
  * Lucia et al. (2013) Lucia, S., Finkler, T., and Engell, S. (2013).  Multi-stage nonlinear model predictive control applied to a semi-batch polymerization reactor under uncertainty.  _Journal of Process Control_ , 23, 1306–1319.
  * Lucia et al. (2020) Lucia, S., Subramanian, S., Limon, D., and Engell, S. (2020).  Stability properties of multi-stage nonlinear model predictive control.  _Systems & Control Letters_, 143, 104743.
  * Mayne et al. (2011) Mayne, D.Q., Kerrigan, E.C., van Wyk, E.J., and Falugi, P. (2011).  Tube-based robust nonlinear model predictive control.  _International Journal of Robust and Nonlinear Control_ , 21(11), 1341–1353.
  * Muñoz de la Peña et al. (2005) Muñoz de la Peña, D., Bemporad, A., and Alamo, T. (2005).  Stochastic programming applied to model predictive control.  In _Proc. of the 44th IEEE Conference on Decision and Control_ , 1361–1366.
  * Puschke and Mitsos (2018) Puschke, J. and Mitsos, A. (2018).  Robust feasible control based on multi-stage enmpc considering worst-case scenarios.  _Journal of Process Control_ , 69, 8–15.
  * Rawlings and Mayne (2009) Rawlings, J. and Mayne, D. (2009).  _Model Predictive Control Theory and Design_.  Nob Hill Pub.
  * Rubagotti et al. (2011) Rubagotti, M., Raimondo, D.M., Ferrara, A., and Magni, L. (2011).  Robust model predictive control with integral sliding mode in continuous-time sampled-data nonlinear systems.  _IEEE Transactions on Automatic Control_ , 56(3), 556–570.
  * Scokaert and Mayne (1998) Scokaert, P. and Mayne, D. (1998).  Min-max feedback model predictive control for constrained linear systems.  _IEEE Transactions on Automatic Control_ , 43(8), 1136–1142.
  * Subramanian et al. (2018) Subramanian, S., S. Lucia, and Engell, S. (2018).  A Synergistic Approach to Robust Output Feedback Control: Tube-based Multi-stage NMPC.  In _Proc. of the International Symposium on Advanced Control of Chemical Processes. IFAC-PapersOnLine_ , 500–505.
  * Subramanian et al. (2022) Subramanian, S., Abdelsalam, Y., Lucia, S., and Engell, S. (2022).  Robust tube-enhanced multi-stage nmpc with stability guarantees.  _IEEE Control Systems Letters_ , 6, 1112–1117.
  * Subramanian et al. (2021) Subramanian, S., Lucia, S., Paulen, R., and Engell, S. (2021).  Tube-enhanced multi-stage model predictive control for flexible robust control of constrained linear systems with additive and parametric uncertainties.  _International Journal of Robust and Nonlinear Control_ , 31(9), 4458–4487.
  * Thombre et al. (2021) Thombre, M., Yu, Z.J., Jäschke, J., and Biegler, L.T. (2021).  Sensitivity-assisted multistage nonlinear model predictive control: Robustness, stability and computational efficiency.  _Computers & Chemical Engineering_, 148, 107269.
  * Villanueva et al. (2017) Villanueva, M.E., Li, J.C., Feng, X., Chachuat, B., and Houska, B. (2017).  Computing ellipsoidal robust forward invariant tubes for nonlinear mpc.  _IFAC-PapersOnLine_ , 50(1), 7175 – 7180.  20th IFAC World Congress.
  * Wächter and Biegler (2006) Wächter, A. and Biegler, L. (2006).  On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming.  _Mathematical Programming_ , 106, 25–57.


# Twistor Space Origins of the Newman-Penrose Map

Kara Farnsworth  Michael L. Graesser  and Gabriel Herczeg

###### Abstract

Recently, we introduced the “Newman-Penrose map,” a novel correspondence
between a certain class of solutions of Einstein’s equations and self-dual
solutions of the vacuum Maxwell equations, which we showed was closely related
to the classical double copy. Here, we give an alternative definition of this
correspondence in terms of quantities that are defined naturally on twistor
space, and a shear-free null geodesic congruence on Minkowski space whose
twistorial character is articulated by the Kerr theorem. The advantage of this
reformulation is that it is purely geometrical in nature, being manifestly
invariant under both spacetime diffeomorphisms and projective transformations
on twistor space. While the original formulation of the map may be more
convenient for most explicit calculations, the twistorial formulation we
present here may be of greater theoretical utility.

##  1 Introduction

Some of our deepest insights into gauge theories and gravity have emerged from
the study of correspondences relating the two. One such promising relationship
is the double copy, which can be summarized schematically by the equation

$\textrm{gravity}=\textrm{gauge}\otimes\textrm{gauge}.$ (1.1)

This general structure has been revealed in many contexts over the past few
decades. The prehistory of the double copy began in 1986, when Kawai, Lewellen
and Tye [1] realized that there was a close relationship between the tree-
level amplitude of a closed string and the square of the tree-level amplitude
of an open string. The subject began in earnest a little more than twenty
years later, when Bern, Carrasco and Johansson [2] observed that a similar
relationship existed between perturbative amplitudes in gauge theories and
gravity. Their results have now been proven at tree-level [3], and a growing
body of evidence suggests that they can be extended beyond tree-level in
perturbation theory [4, 5, 6, 7, 8, 9, 10, 11]. With such promising evidence
at the amplitude level, it is natural to look for a classical formulation of
the double copy in order to understand how fundamental this relationship
between gauge theory and gravity is.

Although there have been attempts at a Lagrangian formulation [12, 13, 14],
much recent success has been found realizing double copies in a classical
setting, as maps between exact solutions of gauge theories and gravity.
Several such maps have been proposed and explored [15, 16, 17, 18, 19, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40],
and many of them share the common feature that null geodesic congruences play
a prominent role. For example, in the Kerr-Schild double copy of [15], the
gauge field is itself the dual of the tangent vector to a geodesic null
congruence. The Newman-Penrose map introduced in [16] also has this flavor,
where the gravity solutions one considers are of Kerr-Schild type, and the
Kerr-Schild vector is assumed to be tangent to a shear-free null geodesic
congruence (SNGC) with non-vanishing expansion. In the Weyl double copy of
[17], the existence of an SNGC is guaranteed by virtue of the Goldberg-Sachs
theorem for algebraically special spacetimes [41]. In light of the Kerr
Theorem [42, 43, 44, 45], which realizes any SNGC in Minkowski space as the
solution of an equation defined by the vanishing of a holomorphic, homogeneous
function of twistor variables, it is perhaps not surprising that the Weyl
double copy was recently given an elegant twistorial formulation [46, 47] that
can also accommodate algebraically special _linearized_ solutions of _any_
Petrov type.

One might wonder whether the Newman-Penrose map can also be given a twistorial
formulation. The primary purpose of this article is to answer this question in
the affirmative: We provide a twistorial formulation of the Newman-Penrose map
that is manifestly invariant under both spacetime diffeomorphisms and
projective transformations on twistor space.

This paper is organized as follows: In Section 2 we review the Newman-Penrose
map, originally defined in [16]. In Section 3 we give an overview of the
relevant aspects of the two-component spinor and twistor formalisms, and in
Section 4 we present the twistor space version of the Newman-Penrose map.
Finally, we conclude in Section 5 and discuss the consequences of the choice
of SNGC normalization in Appendix A.

##  2 The Newman-Penrose map

###  2.1 Kerr-Schild spacetimes

Here we review the salient features of the Newman-Penrose map as defined in
[16]. We begin by recalling some details about Kerr-Schild spacetimes and
their construction in the Newman-Penrose formalism. A Kerr-Schild solution of
Einstein’s equations is one that can be written in the form

$g_{\mu\nu}=\eta_{\mu\nu}+Vl_{\mu}l_{\nu},$ (2.1)

where $\eta_{\mu\nu}$ is a flat metric, $l_{\mu}$ is a null vector111Note that
if $l_{\mu}$ is null with respect to either $g_{\mu\nu}$ or $\eta_{\mu\nu}$
then it is null with respect to both., and $V$ is a scalar function. If we
assume that $l^{\mu}$ is geodesic, shear-free, and expanding (see section 3.2
of [16] for a concise discussion of shear and expansion) then we can write
simple expressions defining a null tetrad for the metric (2.1) [48, 49]:

$\displaystyle l$ $\displaystyle=$
$\displaystyle\partial_{v}-\Phi\partial_{\zeta}-\bar{\Phi}\partial_{\bar{\zeta}}+\Phi\bar{\Phi}\partial_{u}\,$
$\displaystyle n$ $\displaystyle=$
$\displaystyle\partial_{u}-\tfrac{1}{2}Vl\,$ $\displaystyle m$
$\displaystyle=$ $\displaystyle\partial_{\bar{\zeta}}-\Phi\partial_{u}\,$
$\displaystyle\bar{m}$ $\displaystyle=$
$\displaystyle\partial_{\zeta}-\bar{\Phi}\partial_{u},$ (2.2)

where $u,v$ are real light-cone coordinates, and $\zeta,\bar{\zeta}$ are
complex conjugate coordinates related to the usual Cartesian coordinates
$(t,x,y,z)$ by

$u=\tfrac{1}{\sqrt{2}}(t-z),\quad
v=\tfrac{1}{\sqrt{2}}(t+z),\quad\zeta=\tfrac{1}{\sqrt{2}}(x+iy),$ (2.3)

and $\Phi(u,v,\zeta,\bar{\zeta})$ is a complex scalar. In these coordinates,
the flat metric takes the form

$\eta_{\mu\nu}dx^{\mu}dx^{\nu}=2(dudv-d\zeta d\bar{\zeta}).$ (2.4)

In terms of the null tetrad (2.2), we can write the inverse metric as

$g^{\mu\nu}=l^{\mu}n^{\nu}+n^{\mu}l^{\nu}-m^{\mu}\bar{m}^{\nu}-\bar{m}^{\mu}m^{\nu}=\eta^{\mu\nu}-Vl^{\mu}l^{\nu}.$
(2.5)

The shear-free and geodesic conditions on $l^{\mu}$ are equivalent to the
following nonlinear partial differential equations for $\Phi$:

$\Phi_{,v}=\Phi\Phi_{,\zeta}\,,\qquad\Phi_{,\bar{\zeta}}=\Phi\Phi_{,u}.$ (2.6)

The equations (2.6) together imply that the complex scalar $\Phi$ is harmonic
with respect to the flat metric $\eta_{\mu\nu}$, so that

$\Box_{0}\Phi=2(\Phi_{,uv}-\Phi_{,\zeta\bar{\zeta}})=0.$ (2.7)

When the spacetime solves the vacuum Einstein equations, the function $V$ can
be solved for in terms of $\Phi$, so that any $\Phi$ satisfying (2.6)
completely specifies the solution up to constants of integration.

###  2.2 Newman-Penrose map

The spin-raising operator

$\hat{k}=-\frac{Q}{2\pi\epsilon_{0}}(dv\,\partial_{\zeta}+d\bar{\zeta}\,\partial_{u})\,$
(2.8)

was used in [15, 22] to define the self-dual double copy. Here we have
included the constants $Q$, representing the charge, and $\epsilon_{0}$,
representing the vacuum permittivity, to directly map to solutions of
Maxwell’s equations. In that work, the authors justified the form of the
operator by observing that it satisfied certain properties that made it
analogous to a “momentum-space” version of the Kerr-Schild vector $l^{\mu}$.
However, this is not the unique operator possessing those properties, and the
particular choice (2.8) is not justified over any other similar choice of
operator. In fact, other possible choices of operator can easily be found,
e.g., by substituting either $u\leftrightarrow v$, or
$\zeta\leftrightarrow\bar{\zeta}$, or both, in (2.8).222 In [40], it was
pointed out that this operator, which they used to study the symmetries of the
self-dual sectors of Yang-Mills theory and gravity, is closely related to the
area form on certain two-dimensional (complex) subspaces of (complex)
Minkowski space.

It is not a new observation that spin-raising operators such as the one
appearing in (2.8) map solutions of the wave equation (often called
‘Hertz’-type potentials in this context) to solutions of Maxwell’s equations.
In addition to their appearance in the self-dual double copy [15, 22], spin-
raising operators of a very similar form have also appeared in [50], in which
scattering of plane wave solutions was studied in the context of the double
copy. Much earlier, a rather general spinorial version of the spin-raising
operators that takes a solution of the zero rest-mass equations with spin $s$
to a solution of the zero rest-mass equations with spin $s+1/2$ was described
by Penrose in [51] (see also [52]); in fact, this construction was used to
give a spinorial description of the Newman-Penrose map in appendix B of [16].
It has also been known for some time that spin-raising and spin-lowering
operators that map solutions of the zero rest-mass equations to solutions with
arbitrarily raised or lowered spins can be associated with a choice of twistor
of an appropriate rank and type (see section 6.4 of [53], for example). While
all of these structures are related to the Newman-Penrose map, we emphasize
that the identification of a complex solution of the wave equation associated
with a given _real_ solution of the Einstein equations of Kerr-Schild type,
upon which such a spin-raising operator may act to generate a solution of the
vacuum Maxwell equations, is a novel feature of the Newman-Penrose map. The
connection with twistor theory that we present here is similarly novel: while
the twistorial spin-raising operators of [53] depend on an arbitrary choice of
twistor, the construction we present here is defined for a subset of twistors,
and is _independent_ of the choice of twistor within that subset.

In [16] we used this operator to define the Newman-Penrose map as follows:
Given a Kerr-Schild spacetime with an expanding, shear-free Kerr-Schild
vector, we can fix the null tetrad in the form (2.2), and read off the complex
function $\Phi$, which satisfies the non-linear partial differential equations
(2.6) and consequently is harmonic. Then the gauge field defined by

$A=\hat{k}\Phi$ (2.9)

is necessarily a self-dual solution of the vacuum Maxwell equations, since we
have

$A=-\frac{Q}{2\pi\epsilon_{0}}\left(\Phi_{,\zeta}dv+\Phi_{,u}d\bar{\zeta}\right),$
(2.10)

from which we can compute the field strength two-form $F=dA$:

$\displaystyle F=-\frac{Q}{2\pi\epsilon_{0}}\left(\Phi_{,u\zeta}du\wedge
dv-\Phi_{,\zeta\zeta}dv\wedge d\zeta+\Phi_{,uu}du\wedge
d\bar{\zeta}+\Phi_{,u\zeta}d\zeta\wedge d\bar{\zeta}\right),$ (2.11)

where we used
$\Phi_{,uv}-\Phi_{,\zeta\bar{\zeta}}=\tfrac{1}{2}\Box_{0}\Phi=0$. The field
strength in (2.11) is self-dual with respect to $\star_{0}$, the Hodge star
operator associated with the background metric $\eta_{\mu\nu}$—that is, $F$
satisfies

$F=i\star_{0}\\!F,$ (2.12)

which, together with the fact that $F$ is exact, implies the vacuum flat-space
Maxwell equation

$d\star_{0}\\!F=0.$ (2.13)

The primary purpose of this article is to give a purely geometric
interpretation of the origin of the spin-raising operator (2.8) and $\Phi$ in
terms of structures naturally defined on projective twistor space. This gives
a geometric foundation to the coordinate-dependent description which we have
reviewed in this section, and we hope that it will provide some additional
insight on the self-dual double copy of [15, 22], and the recent discovery of
the twistorial origins of the Weyl double copy based on the Penrose transform
[46, 47].

##  3 Aspects of twistor theory

Twistors are natural objects for studying null geodesics in complexified
Minkowski space $\mathbb{CM}$. Here we briefly review the relevant aspects of
twistor theory, focusing on the relationships between geometric structures on
$\mathbb{CM}$, and the corresponding geometric strucures on _twistor space_
$\mathbb{T}$. Our discussion follows closely the approach taken in chapter 7
of [54]. See also [53, 55, 56] for more background and [57] for a more recent
review of twistor theory.

###  3.1 Spin space

In four spacetime dimensions, Minkowski space—or more generally, the tangent
space $T_{p}\mathbb{N}$ to any point $p$ of a pseudo-Riemannian manifold
$\mathbb{N}$ of signature $(3,1)$—is isomorphic to the set of $2\times 2$
Hermitian matrices. One can make this isomorphism explicit by specifying a
soldering form333For the considerations of the present article, it will
suffice to consider a particular soldering form that is adapted to Minkowski
space. Solderings for a general spacetime can be specified by contracting the
flat-space soldering form with an orthonormal frame for the curved spacetime
metric, but we will not need such structures here. $\bm{\sigma}_{\mu}$ with
components $\sigma_{\mu}^{AA^{\prime}}$, which assigns a Hermitian matrix
$\bm{V}\\!=\bm{\sigma}_{\mu}V^{\mu}$ with components
$V^{AA^{\prime}}=\sigma_{\mu}^{AA^{\prime}}V^{\mu}$ to every real spacetime
vector $V^{\mu}$. Using standard Cartesian coordinates on Minkowski space, a
common choice is $\bm{\sigma}_{\mu}=(\bf{I}\\!$ , $\\!\bm{\sigma}_{i})$, where
$\bm{\sigma}_{i}$ are the Pauli matrices normalized so that
$\\{\bm{\sigma}_{i},\bm{\sigma}_{j}\\}=\delta_{ij}\bf{I}$. Such a soldering
form maps an arbitrary vector $V^{\mu}=(V^{0},V^{i})$ to the Hermitian matrix

$\bm{V}=\frac{1}{\sqrt{2}}\begin{pmatrix}V^{0}+V^{3}&V^{1}+iV^{2}\\\
V^{1}-iV^{2}&V^{0}-V^{3}\end{pmatrix},$ (3.1)

with $\textrm{det}(\bf{V})$ $=\frac{1}{2}|V|^{2}$. Then for any $\bf{U}$ $\in
SL(2,\mathbb{C})$, the Hermitian matrix
$\bf{V}^{\prime}=\bf{U}^{\dagger}\bf{V}\bf{U}$, has determinant
$\textrm{det}(\bf{V}^{\prime})$ $=\textrm{det}(\bf{V})$ $=\frac{1}{2}|V|^{2}$,
meaning $\bf{U}$ corresponds to a linear transformation preserving the norm of
vectors on Minkowski space, i.e. a Lorentz transformation.

We define _spin space_ , $\mathbb{S}$, as the space of complex-valued, two-
component column vectors on which the elements $\mathbf{U}$ act by left
multiplication. A vector $\alpha\in\mathbb{S}$ is called a _spinor_ , and it’s
components are denoted by $\alpha^{A}$. Similarly, the complex conjugate space
is denoted by $\mathbb{\bar{S}}$, and the components of a conjugate spinor are
written with primed indices, e.g. $\beta^{A^{\prime}}$.

Given any spinor $\alpha^{A}$, one can consider the Hermitian matrix with
components $\alpha^{A}\bar{\alpha}^{A^{\prime}}$, which has vanishing
determinant, so that $\alpha^{A}\bar{\alpha}^{A^{\prime}}$ corresponds to a
real, null vector. Thus, a spinor can be viewed as the “square root of a null
vector.” More generally, given $\lambda^{A}\in\mathbb{S}$,
$\mu^{A^{\prime}}\in\bar{\mathbb{S}}$, the non-Hermitian matrix
$\lambda^{A}\mu^{A^{\prime}}$ corresponds to a complex null vector. Spinor
indices are raised and lowered using the antisymmetric, rank two spinor tensor
$\epsilon_{AB}$ via

$\displaystyle\alpha_{B}=\alpha^{A}\epsilon_{AB},\qquad\alpha^{A}=\epsilon^{AB}\alpha_{B}$
(3.2)

and every spinor satisfies $\alpha_{A}\alpha^{A}=0$. It is often convenient to
fix a _normalized dyad_ for spin space $\mathbb{S}$, i.e. a pair of spinors
$(o^{A},\iota^{A})$ satisfying $o_{A}\iota^{A}=1$. A normalized dyad on
$\mathbb{S}$ naturally defines a null tetrad on $\mathbb{M}$ via

$l_{0}^{AA^{\prime}}=o^{A}\bar{o}^{A^{\prime}},\quad
n_{0}^{AA^{\prime}}=\iota^{A}\bar{\iota}^{A^{\prime}},\quad
m_{0}^{AA^{\prime}}=o^{A}\bar{\iota}^{A^{\prime}},$ (3.3)

where we have decorated the tetrad elements with subscripts to emphasize that
this is a tetrad for flat Minkowski space $\mathbb{M}$, in contrast to the
tetrad (2.2) for the full, Kerr-Schild spacetime.

###  3.2 Twistor space

So far, the discussion has focused only on spinors at a point, but all of
these structures can be generalized to smoothly varying functions of
$\mathbb{CM}$, or more general four-dimensional manifolds. In particular, we
may consider _spinor fields_ $\alpha^{A}(x)$ over $\mathbb{CM}$. The covariant
derivative $\nabla_{\mu}$ on $\mathbb{CM}$ can be extended to a covariant
derivative on spinor fields satisfying
$\nabla_{AA^{\prime}}\epsilon_{BC}=0=\nabla_{AA^{\prime}}\epsilon_{B^{\prime}C^{\prime}}$.

A _twistor_ is a spinor field $\Omega^{A}(x)$ satisfying the twistor equation

$\nabla_{A^{\prime}}{}^{(A}\Omega^{B)}=0.$ (3.4)

Over complexified Minkowski space $\mathbb{CM}$, (3.4) can be solved exactly
to give

$\Omega^{A}=\omega^{A}-ix^{AA^{\prime}}\pi_{A^{\prime}},$ (3.5)

where $\omega^{A}$ and $\pi_{A^{\prime}}$ are constant spinors. Thus, the
space of solutions to the twistor equation—i.e., the twistor space,
$\mathbb{T}$—is coordinatized by a pair of spinors
$Z=(\omega^{A},\pi_{A^{\prime}})$ and we can regard $\mathbb{T}$ as a four-
dimensional complex vector space. Twistor space can be equipped with a natural
Hermitian inner product, and in particular, each twistor $Z\in\mathbb{T}$ can
be assigned a norm

$|Z|^{2}:=\omega^{A}\bar{\pi}_{A}+\bar{\omega}^{A^{\prime}}\pi_{A^{\prime}}.$
(3.6)

The (squared) norm of a twistor can be positive, negative or zero, so we
obtain a decomposition of the twistor space
${\mathbb{T}=\mathbb{T}_{+}\cup\,\mathbb{T}_{0}\,\cup\,\mathbb{T}_{-}}$, where
$\mathbb{T}_{+}$ is the set of twistors with positive norm, $\mathbb{T}_{0}$
is the set of twistors with zero norm, which are called _null twistors_ , and
$\mathbb{T}_{-}$ is the set of twistors with negative norm. In what follows,
we will be particularly interested in null twistors, which we will see
determine real, null geodesics in $\mathbb{CM}$.

In order to understand the relationship between null twistors and real, null
geodesics, consider the subset of $\mathbb{CM}\,$ defined by the zero locus of
an arbitrary twistor

$\Omega^{A}(x)=0.$ (3.7)

The general solution (3.5) then specifies this subset as the points
$x^{AA^{\prime}}\in\mathbb{CM}$ that satisfy the incidence relation

$\omega^{A}=ix^{AA^{\prime}}\pi_{A^{\prime}}.$ (3.8)

Given a particular solution $x_{0}^{AA^{\prime}}$ to (3.8), the general
solution can be written as

$x^{AA^{\prime}}=x_{0}^{AA^{\prime}}+\lambda^{A}\pi^{A^{\prime}},$ (3.9)

where $\lambda^{A}$ is an arbitrary spinor. Equation (3.9) describes a totally
null two-plane in $\mathbb{CM}$, called an _$\alpha$ -plane_. Every tangent
vector to an $\alpha$-plane is null and orthogonal to every other tangent
vector. Such a plane cannot exist in real Minkowski space, so the points
contained in the plane will generically be complex.

Since (3.9) is the general solution to (3.8), the twistor
$tZ=(t\omega^{A},t\pi_{A^{\prime}})$ defines the same $\alpha$-plane as $Z$
for any $t\in\mathbb{C}$. The equivalence relation $\,\sim\,$ on $\mathbb{T}$,
defined by $tZ\sim Z$ for some $t\in\mathbb{C}$, gives rise to the notion of
_projective twistor space_ $\mathbb{PT}:=\mathbb{T}/\\!\sim$. We may then
regard the points of $\mathbb{PT}$ as being labeled by $\alpha$-planes
444$\mathbb{PT}$ is a 3-dimensional complex manifold and can be realized as an
open subset of $\mathbb{CP}^{3}$. The projective space $\mathbb{PT}_{0}$ of
null twistors on $\mathbb{PT}$ is a five-dimensional real manifold.. The
intersection of the $\alpha$-plane of a null twistor $Z\in\mathbb{T}_{0}$ and
the real subset $\mathbb{M}$ of $\mathbb{CM}$ is described by a real, null
geodesic, as we illustrate in Fig. 1 and now show.

Figure 1: $\mathbb{M}$ and an $\alpha$-plane $\\{\Omega^{A}(x)=0\\}$ as
subsets of $\mathbb{CM}$. Their intersection is a real light ray $\gamma$.

We have already observed that the points on an $\alpha$-plane are generically
complex. Suppose that some twistor $Z$ defines an $\alpha$-plane that contains
at least one real point. Then, without loss of generality, let us assume
$x_{0}^{AA^{\prime}}$ is real. Contracting (3.8) with $\bar{\pi}_{A}$ then
gives

$\omega^{A}\bar{\pi}_{A}=ix_{0}^{AA^{\prime}}\bar{\pi}_{A}\pi_{A^{\prime}}.$
(3.10)

Since $x_{0}^{AA^{\prime}}$ is real, the right hand side of (3.10) is clearly
imaginary:

$\omega^{A}\bar{\pi}_{A}=ia,\qquad a\in\mathbb{R}.$ (3.11)

Then we immediately find that

$|Z|^{2}=\omega^{A}\bar{\pi}_{A}+\bar{\omega}^{A^{\prime}}\pi_{A^{\prime}}=0.$
(3.12)

Thus, whenever the $\alpha$-plane defined by a twistor $Z$ contains a real
point, $Z$ is null. The converse can be proven as well (see chapter 7 of [54]
for details) so $Z$ is null if and only if its $\alpha$-plane contains at
least one real point. Moreover, if $x_{0}^{AA^{\prime}}$ is real, then the
$\alpha$-plane contains the whole real null geodesic

$x^{AA^{\prime}}=x_{0}^{AA^{\prime}}+r\bar{\pi}^{A}\pi^{A^{\prime}},\quad
r\in\mathbb{R},$ (3.13)

which establishes that null twistors determine real null geodesics in
$\mathbb{CM}$, as claimed. In the following section we introduce a twistorial
definition of the Newman-Penrose map that is projectively invariant, and so
may be defined directly on $\mathbb{PT}$.

##  4 Newman-Penrose map from twistor space geometry

Consider a null twistor $Z=(\omega^{A},\pi_{A^{\prime}})$ with the property
that

$l_{AA^{\prime}}\bar{\pi}^{A}\pi^{A^{\prime}}=1,\qquad$ (4.1)

where $l_{AA^{\prime}}$ is the expanding, null, shear-free and geodesic Kerr-
Schild vector. Note that it is always possible to choose a scale for
$l_{AA^{\prime}}$ so that such a twistor exists, and we suppose that such a
scale has been chosen. In fact, as we shall see explicitly below, this
condition defines a 4-dimensional real subspace of $\mathbb{T}_{0}$. Moreover
the condition (4.1) has a clear geometric interpretation. The vector
$\bar{\pi}^{A}\pi^{A^{\prime}}$ is a tangent vector to the real null geodesic
contained in the $\alpha$-plane corresponding to $Z$. This can be extended to
all of spacetime via parallel transport with respect to the flat metric
$\eta_{\mu\nu}$. The normalization condition then requires that this tangent
vector be everywhere transverse to $l_{AA^{\prime}}$ and scaled so that their
inner product is normalized to unity. Kerr’s Theorem [42, 43, 44, 45] implies
that, without loss of generality, we can write 555This corresponds to a shift
in the dyad (3.3) of
$o_{A}\rightarrow\tilde{o}_{A}=o_{A}-\bar{\Phi}\iota_{A}$, as explained in
[58] and appendix A.2 of [16].

$l_{AA^{\prime}}dx^{AA^{\prime}}=du+\bar{\Phi}d\zeta+\Phi
d\bar{\zeta}+\Phi\bar{\Phi}dv,$ (4.2)

in which case the normalization condition (4.1) implies that
$\bar{\pi}^{A}\pi^{A^{\prime}}\partial_{AA^{\prime}}=\partial_{u}$. It is also
possible to choose the normalization so that the coefficient of $dv$ is equal
to one in (4.2), however we show in the Appendix that this choice leads to the
same field strength we obtain using the current normalization.

In order to write explicit expressions in terms of $\omega^{A},$ and
$\pi_{A^{\prime}}$, we introduce a constant, normalized dyad
$o_{A}\iota^{A}=1$ with

$o^{A}\bar{o}^{A^{\prime}}\partial_{AA^{\prime}}=\partial_{v},\quad
o^{A}\bar{\iota}^{A^{\prime}}\partial_{AA^{\prime}}=\partial_{\zeta},\quad\iota^{A}\bar{\iota}^{A^{\prime}}\partial_{AA^{\prime}}=\partial_{u}.$
(4.3)

Note that the dyad $(o^{A},\iota^{A})$ is tied to a tetrad for flat Minkowski
space, in contrast to the tetrad for the full Kerr-Schild spacetime (2.2). In
this basis, the null Kerr-Schild vector (4.2) can be written as
$l_{AA^{\prime}}=\tilde{o}_{A}\overline{\tilde{o}}_{A^{\prime}}=\left(o_{A}-\bar{\Phi}\iota_{A}\right)\left(\bar{o}_{A^{\prime}}-\Phi\bar{\iota}_{A^{\prime}}\right)$.

Now we can expand the constant spinors $\omega^{A},$ and $\pi_{A^{\prime}}$ as

$\omega^{A}=\alpha
o^{A}+\beta\iota^{A},\qquad\pi_{A^{\prime}}=\gamma\bar{o}_{A^{\prime}}+\delta{\bar{\iota}}_{A^{\prime}},$
(4.4)

where $\alpha,\beta,\gamma,\delta\in\mathbb{C}$.666These parameters are not
the spin coefficients of the tetrad. The normalization condition (4.1) then
tells us that

$\gamma=0,\qquad|\delta|^{2}=1.$ (4.5)

Moreover, since $Z$ is null, we have

$\omega^{A}\bar{\pi}_{A}=ia,\quad a\in\mathbb{R},$ (4.6)

which implies that $\alpha=-ia\delta$. Now our expessions for the spinor
components of $Z$ become

$\omega^{A}=-i\delta
ao^{A}+\beta\iota^{A},\qquad\pi_{A^{\prime}}=\delta\bar{\iota}_{A^{\prime}}.$
(4.7)

In the generic case when $a\neq 0$, the $\alpha$-plane associated with $Z$ can
be described by the equation

$\displaystyle
x^{AA^{\prime}}=\frac{1}{a}\omega^{A}\bar{\omega}^{A^{\prime}}+\lambda^{A}\pi^{A^{\prime}}$
(4.8)

for arbitrary $\lambda^{A}$. The special case $a=0$ requires some care, as it
includes both ordinary $\alpha$-planes and $\alpha$-planes that contain real
geodesics that are generators of null infinity rather than geodesics on the
interior of Minkowski space. Given these subtleties, we will only consider
twistors for which $a\neq 0$ from now on. For reasons that will soon be clear,
we expand $\lambda^{A}$ in terms of the spinor dyad as

$\lambda^{A}=\bar{\delta}\left(\zeta
o^{A}+\left(u+u_{0}-\tfrac{|\beta|^{2}}{a}\right)\iota^{A}\right),$ (4.9)

where $u_{0}$ is a real constant, and $u$ and $\zeta$ vary over all complex
values. Using this expansion for $\lambda^{A}$, the equation for the
$\alpha$-plane becomes

$x^{AA^{\prime}}=v_{0}o^{A}\bar{o}^{A^{\prime}}+(\zeta+\zeta_{0})o^{A}\bar{\iota}^{A^{\prime}}+\bar{\zeta}_{0}\iota^{A}\bar{o}^{A^{\prime}}+(u+u_{0})\iota^{A}\bar{\iota}^{A^{\prime}},$
(4.10)

where we have identified $a=v_{0}$, $\beta=-i\delta\bar{\zeta}_{0}$. This is
then the $\alpha$-plane defined by $v=v_{0}$, $\bar{\zeta}=\bar{\zeta}_{0}$,
and the tangent bivector $\tau$ to any such plane satisfies the
proportionality

$\tau\propto\partial_{u}\wedge\partial_{\zeta}.$ (4.11)

Lowering an index with the flat metric leads to

$\tau^{\mu\rho}\eta_{\rho\nu}dx^{\nu}\partial_{\mu}\propto-(d\bar{\zeta}\partial_{u}+dv\partial_{\zeta})\propto\hat{k}.$
(4.12)

The above discussion shows that we are able to motivate the form of the
operator $\hat{k}$ geometrically from twistor space. However, we would like to
describe _all_ quantities in the NP map in terms of objects naturally
associated to the null twistor $Z$ and the dual to the Kerr-Schild SNGC
$l_{AA^{\prime}}$ without reference to any particular coordinate system. For
this, we can introduce the vectors

$p^{AA^{\prime}}=\bar{\pi}^{A}\pi^{A^{\prime}},\quad
q^{AA^{\prime}}=\omega^{A}\pi^{A^{\prime}}$ (4.13)

which span the tangent space to the $\alpha$-plane. Now we can take

$\tau=\frac{1}{\omega^{A}\bar{\pi}_{A}}p\wedge
q=-\delta^{2}\partial_{u}\wedge\partial_{\zeta},$ (4.14)

and define

$\hat{\kappa}:=-\frac{Q}{2\pi\epsilon_{0}}\tau^{\mu\rho}\eta_{\rho\nu}dx^{\nu}\partial_{\mu}=-\frac{Q}{2\pi\epsilon_{0}}\delta^{2}(dv\,\partial_{\zeta}+d\bar{\zeta}\,\partial_{u})$
(4.15)

which matches the operator (2.8) up to an overall, arbitrary phase
$\delta^{2}$. Next, we need an invariant definition of $\Phi$. For this, we
note that

$\Psi:=\frac{1}{\omega^{A}\bar{\pi}_{A}}l_{BB^{\prime}}\bar{q}^{BB^{\prime}}=\bar{\delta}^{2}\Phi+\frac{\bar{\delta}\bar{\beta}}{ia}.$
(4.16)

The first term in $\Psi$ differs from $\Phi$ by the constant phase
$\bar{\delta}^{2}$, complementary to the phase of $\hat{\kappa}$. The second
term is a constant, so it is annihilated by $\hat{\kappa}$. Hence, we can
define the NP map in a manifestly invariant manner, depending only on the SNGC
$l_{AA^{\prime}}$ and the twistor $Z$ by

$A=\hat{\kappa}\Psi.$ (4.17)

In particular, note that this definition of $A$ is independent of the phase
$\delta$, the individual dependence in $\Psi$ and $\hat{\kappa}$ having
canceled, and is equivalent to the original definition in (2.9). In fact, if
we insist that the normalization condition (4.1) be preserved under rescalings
$Z\to tZ$, $t\in\mathbb{C}$, so that
$l_{AA^{\prime}}\to|t|^{-2}l_{AA^{\prime}}$, then it can be shown that the
Newman-Penrose map depends only on the equivalence class
$[Z]=[tZ],t\in\mathbb{C}$, and so is projectively invariant.

We end with a comment about the uniqueness of $\hat{\kappa}$. In general one
may consider other operators similar to $\hat{\kappa}$, as described in
Section 2.2, and it isn’t clear why this particular form of $\hat{\kappa}$ is
preferred. The twistor space construction of the NP map described here offers
some answers. At the outset, a normalization for $l$ is made to set one of its
coefficients equal to a constant, but which coefficient to choose is seemingly
arbitrary. To keep $l_{AA^{\prime}}$ real, the only choices are to set the
coefficient of $dv$ or $du$ equal to a constant. Choosing $du$ leads to the
operator $\hat{\kappa}$, through the tangent bivector of the $\alpha$-plane
associated to projective twistor $[Z]$ described above. If we instead choose
$dv$, the tangent bivector is proportional to
$\partial_{v}\wedge\partial_{\bar{\zeta}}$, leading to a new differential
operator proportional to $du\,\partial_{\bar{\zeta}}+d\zeta\,\partial_{v}$.
However, the resulting field strength is the same as the one obtained above,
as shown in the Appendix.

##  5 Discussion

The twistorial formalism we have presented here firmly establishes the
geometric foundations of the Newman-Penrose map. Implicit in this construction
is the fact that any SNGC on Minkowski space can be given a purely twistorial
definition in terms of the vanishing set of a homogeneous and holomorphic
function defined on twistor space [42, 43]. In particular, for any Kerr-Schild
spacetime with an expanding SNGC, we have defined a coordinate-independent
construction for both the operator $\hat{\kappa}$ and the complex scalar
$\Psi$ in terms of twistor variables.

These quantities are then used to define the self-dual gauge field in the
Newman-Penrose map, $A=\hat{\kappa}\Psi$, which we showed in previous work
[16] coincides with other double copy prescriptions. The close relationship
between twistor theory and the geometry of SNGCs is at the heart of this
construction, and given the distinguished role that shear-free rays play in
several of the approaches to the classical double copy, it is no wonder that
twistors are increasingly being recognized as a useful tool in this setting.
The work we presented here was partially motivated by applications of twistor
theory to the Weyl double copy [46, 47], where the Penrose transform was of
central importance; however our investigations seem to run somewhat orthogonal
to those references, and it would be very desirable to understand more clearly
the relationship between the two approaches.

Given the work that has been done on extending the classical double copy on
maximally symmetric backgrounds [20, 21], it could also be interesting to
consider maximally symmetric extensions of the twistorial formulation of the
Newman-Penrose map. Since (A)dS shares the same twistor space as flat
spacetime [59], and since the property of being an SNGC is conformally
invariant, it seems reasonable to expect that much of the formalism we have
developed here can be adapted to the maximally symmetric case. It would also
be interesting to extend this formalism to arbitrary asymptotically flat
spacetimes using the generalized Kerr Theorem of [60].

## Acknowledgements

We thank Gilly Elor, Sylvia Nagy, and Paul Tod for discussions and Jake
Herczeg for his help designing figure 1. KF acknowledges support from Simons
Foundation Award Number 658908. The work of MG is supported by the LDRD
program at Los Alamos National Laboratory and by the U.S. Department of
Energy, Office of High Energy Physics, under Contract No. DE-AC52-06NA25396.

##  Appendix A SNGC Normalization

The point of this Appendix is to show that no loss of generality occurs in
normalizing the coefficient of $du$ in the expression (4.2) for
$l_{AA^{\prime}}dx^{AA^{\prime}}$ to be equal to one. To show this, we rescale
$l_{AA^{\prime}}dx^{AA^{\prime}}$ so that its coefficient of $dv$ is
normalized to one. Next, we show that the twistorial Newman-Penrose map
defined in terms of this rescaled vector gives rise to the same field-strength
as in (2.11).

After rescaling, we have

$l_{AA^{\prime}}dx^{AA^{\prime}}=dv+\Phi^{-1}d\zeta+\bar{\Phi}^{-1}d\bar{\zeta}+(\Phi\bar{\Phi})^{-1}du,$
(A.1)

where $\Phi$ satisfies (2.6). Then for a null twistor
$Z=(\omega^{A},\pi_{A^{\prime}})$, the normalization condition

$l_{AA^{\prime}}\bar{\pi}^{A}\pi^{A^{\prime}}=1$ (A.2)

implies $\bar{\pi}^{A}\pi^{A^{\prime}}\partial_{AA^{\prime}}=\partial_{v}$.
Expanding $\omega^{A}$ and $\pi_{A^{\prime}}$ as in (4.4), the condition (A.2)
leads to

$|\gamma|^{2}=1,\qquad\delta=0.$ (A.3)

Since $Z$ is null, $\omega^{A}\bar{\pi}_{A}=ia$ implies $\beta=i\gamma a$,
giving

$\omega^{A}=\alpha o^{A}+i\gamma
a\iota^{A},\qquad\pi_{A^{\prime}}=\gamma\bar{o}_{A^{\prime}}.$ (A.4)

The $\alpha$-plane associated to $Z$ has the same form as (4.8), and here is
defined by $u=u_{0}$ and $\zeta=\zeta_{0}$, where we have identified
$a=u_{0}$, $\alpha=i\gamma\zeta_{0}$. This leads to the tangent bivector

$\tau=\frac{1}{\omega^{A}\bar{\pi}_{A}}p\wedge
q=\gamma^{2}\partial_{v}\wedge\partial_{\bar{\zeta}}\,$ (A.5)

which in turn leads to

$\hat{\kappa}=\frac{Q}{2\pi\epsilon_{0}}\gamma^{2}(d\zeta\partial_{v}+du\partial_{\bar{\zeta}})\,.$
(A.6)

Once expanded in a coordinate spinor basis, the operator $\hat{\kappa}$ has a
different expression than the one in (2.8). One might then expect the
resulting field strength $F$ to be different too. However, the associated
scalar $\Psi$ is also modified, and these two effects compensate such that the
field strength is unchanged.

The definition of $\Psi$ (4.16) in terms of the twistor (A.4) is

$\Psi=-\bar{\gamma}^{2}\Phi^{-1}+\frac{\bar{\alpha}\bar{\gamma}}{ia},$ (A.7)

which gives for the field strength associated with the gauge field
$A=\hat{\kappa}\Psi$,

$F=-\frac{Q}{2\pi\epsilon_{0}}\left(\Phi_{,u\zeta}(du\wedge dv+d\zeta\wedge
d\bar{\zeta})+\Phi_{,uu}du\wedge d\bar{\zeta}-\Phi_{\zeta\zeta}dv\wedge
d\zeta\right)$ (A.8)

which is identical to (2.11).

## References

  * [1] H. Kawai, D.C. Lewellen and S.-H. Tye, _A relation between tree amplitudes of closed and open strings_ , _Nuclear Physics B_ 269 (1986) 1\.
  * [2] Z. Bern, J.J.M. Carrasco and H. Johansson, _New Relations for Gauge-Theory Amplitudes_ , _Phys. Rev._ D78 (2008) 085011 [0805.3993].
  * [3] Z. Bern, T. Dennen, Y.-t. Huang and M. Kiermaier, _Gravity as the Square of Gauge Theory_ , _Phys. Rev._ D82 (2010) 065003 [1004.0693].
  * [4] Z. Bern, J.J.M. Carrasco and H. Johansson, _Perturbative Quantum Gravity as a Double Copy of Gauge Theory_ , _Phys. Rev. Lett._ 105 (2010) 061602 [1004.0476].
  * [5] N.E.J. Bjerrum-Bohr, P.H. Damgaard, R. Monteiro and D. O’Connell, _Algebras for Amplitudes_ , _JHEP_ 06 (2012) 061 [1203.0944].
  * [6] J. Nohle, _Color-Kinematics Duality in One-Loop Four-Gluon Amplitudes with Matter_ , _Phys. Rev._ D90 (2014) 025020 [1309.7416].
  * [7] R.H. Boels, R.S. Isermann, R. Monteiro and D. O’Connell, _Colour-Kinematics Duality for One-Loop Rational Amplitudes_ , _JHEP_ 04 (2013) 107 [1301.4165].
  * [8] Z. Bern, S. Davies, T. Dennen, Y.-t. Huang and J. Nohle, _Color-Kinematics Duality for Pure Yang-Mills and Gravity at One and Two Loops_ , _Phys. Rev._ D92 (2015) 045041 [1303.6605].
  * [9] Y.-J. Du, B. Feng and C.-H. Fu, _Dual-color decompositions at one-loop level in Yang-Mills theory_ , _JHEP_ 06 (2014) 157 [1402.6805].
  * [10] S. Oxburgh and C.D. White, _BCJ duality and the double copy in the soft limit_ , _JHEP_ 02 (2013) 127 [1210.1110].
  * [11] S.G. Naculich, H. Nastase and H.J. Schnitzer, _All-loop infrared-divergent behavior of most-subleading-color gauge-theory amplitudes_ , _JHEP_ 04 (2013) 114 [1301.2234].
  * [12] R. Monteiro and D. O’Connell, _The Kinematic Algebra From the Self-Dual Sector_ , _JHEP_ 07 (2011) 007 [1105.2565].
  * [13] C. Cheung and G.N. Remmen, _Hidden Simplicity of the Gravity Action_ , _JHEP_ 09 (2017) 002 [1705.00626].
  * [14] M. Tolotti and S. Weinzierl, _Construction of an effective Yang-Mills Lagrangian with manifest BCJ duality_ , _JHEP_ 07 (2013) 111 [1306.2975].
  * [15] R. Monteiro, D. O’Connell and C.D. White, _Black holes and the double copy_ , _JHEP_ 12 (2014) 056 [1410.0239].
  * [16] G. Elor, K. Farnsworth, M.L. Graesser and G. Herczeg, _The Newman-Penrose Map and the Classical Double Copy_ , _JHEP_ 12 (2020) 121 [2006.08630].
  * [17] A. Luna, R. Monteiro, I. Nicholson and D. O’Connell, _Type D Spacetimes and the Weyl Double Copy_ , _Class. Quant. Grav._ 36 (2019) 065003 [1810.08183].
  * [18] A. Luna, R. Monteiro, D. O’Connell and C.D. White, _The classical double copy for Taub-NUT spacetime_ , _Phys. Lett._ B750 (2015) 272 [1507.01869].
  * [19] A. Luna, R. Monteiro, I. Nicholson, D. O’Connell and C.D. White, _The double copy: Bremsstrahlung and accelerating black holes_ , _JHEP_ 06 (2016) 023 [1603.05737].
  * [20] M. Carrillo-González, R. Penco and M. Trodden, _The classical double copy in maximally symmetric spacetimes_ , _JHEP_ 04 (2018) 028 [1711.01296].
  * [21] N. Bahjat-Abbas, A. Luna and C.D. White, _The Kerr-Schild double copy in curved spacetime_ , _JHEP_ 12 (2017) 004 [1710.01953].
  * [22] D.S. Berman, E. Chacón, A. Luna and C.D. White, _The self-dual classical double copy, and the Eguchi-Hanson instanton_ , _JHEP_ 01 (2019) 107 [1809.04063].
  * [23] I. Bah, R. Dempsey and P. Weck, _Kerr-Schild Double Copy and Complex Worldlines_ , _JHEP_ 02 (2020) 180 [1910.04197].
  * [24] M. Carrillo González, B. Melcher, K. Ratliff, S. Watson and C.D. White, _The classical double copy in three spacetime dimensions_ , _JHEP_ 07 (2019) 167 [1904.11001].
  * [25] Y.-T. Huang, U. Kol and D. O’Connell, _The Double Copy of Electric-Magnetic Duality_ , 1911.06318.
  * [26] R. Alawadhi, D.S. Berman, B. Spence and D. Peinador Veiga, _S-duality and the double copy_ , _JHEP_ 03 (2020) 059 [1911.06797].
  * [27] N. Bahjat-Abbas, R. Stark-Muchão and C.D. White, _Monopoles, shockwaves and the classical double copy_ , _JHEP_ 04 (2020) 102 [2001.09918].
  * [28] A. Luna, S. Nagy and C. White, _The convolutional double copy: a case study with a point_ , 2004.11254.
  * [29] K. Lee, _Kerr-Schild Double Field Theory and Classical Double Copy_ , _JHEP_ 10 (2018) 027 [1807.08443].
  * [30] W. Cho and K. Lee, _Heterotic Kerr-Schild Double Field Theory and Classical Double Copy_ , _JHEP_ 07 (2019) 030 [1904.11650].
  * [31] K. Kim, K. Lee, R. Monteiro, I. Nicholson and D. Peinador Veiga, _The Classical Double Copy of a Point Charge_ , _JHEP_ 02 (2020) 046 [1912.02177].
  * [32] A. Banerjee, E.O. Colgáin, J.A. Rosabal and H. Yavartanoo, _Ehlers as EM duality in the double copy_ , _Phys. Rev. D_ 102 (2020) 126017 [1912.02597].
  * [33] N. Arkani-Hamed, Y.-t. Huang and D. O’Connell, _Kerr black holes as elementary particles_ , _JHEP_ 01 (2020) 046 [1906.10100].
  * [34] D.A. Easson, C. Keeler and T. Manton, _Classical double copy of nonsingular black holes_ , _Phys. Rev. D_ 102 (2020) 086015 [2007.16186].
  * [35] R. Monteiro, D. O’Connell, D.P. Veiga and M. Sergola, _Classical Solutions and their Double Copy in Split Signature_ , 2012.11190.
  * [36] H. Godazgar, M. Godazgar, R. Monteiro, D.P. Veiga and C.N. Pope, _Weyl Double Copy for Gravitational Waves_ , _Phys. Rev. Lett._ 126 (2021) 101103 [2010.02925].
  * [37] D.S. Berman, K. Kim and K. Lee, _The classical double copy for M-theory from a Kerr-Schild ansatz for exceptional field theory_ , _JHEP_ 04 (2021) 071 [2010.08255].
  * [38] L. Alfonsi, C.D. White and S. Wikeley, _Topology and Wilson lines: global aspects of the double copy_ , 2004.07181.
  * [39] G. Alkac, M.K. Gumus and M. Tek, _The Classical Double Copy in Curved Spacetime_ , 2103.06986.
  * [40] M. Campiglia and S. Nagy, _A double copy for asymptotic symmetries in the self-dual sector_ , _JHEP_ 03 (2021) 262 [2102.01680].
  * [41] J. Goldberg and R. Sachs, _Republication of: A theorem on petrov types_ , _General Relativity and Gravitation_ 41 (2009) 433.
  * [42] R. Penrose, _Twistor algebra_ , _J. Math. Phys._ 8 (1967) 345.
  * [43] R. Penrose, _Twistor quantization and curved space-time_ , _Int. J. Theor. Phys._ 1 (1968) 61.
  * [44] D. Cox, _Kerr’s theorem and the kerr–schild congruences_ , _Journal of Mathematical Physics_ 18 (1977) 1188 [https://doi.org/10.1063/1.523388].
  * [45] D. Cox and E.J. Flaherty, _A conventional proof of Kerr’s Theorem_ , _Communications in Mathematical Physics_ 47 (1976) 75.
  * [46] C.D. White, _Twistorial Foundation for the Classical Double Copy_ , _Phys. Rev. Lett._ 126 (2021) 061602 [2012.02479].
  * [47] E. Chacón, S. Nagy and C.D. White, _The Weyl double copy from twistor space_ , 2103.16441.
  * [48] C. McIntosh et al., _Kerr-schild spacetimes revisited_ , in _Conference on Mathematical Relativity_ , pp. 201–206, Centre for Mathematics and its Applications, Mathematical Sciences Institute, 1989.
  * [49] C. McIntosh and M. Hickman, _Single kerr-schild metrics: a double view_ , _General relativity and gravitation_ 20 (1988) 793.
  * [50] T. Adamo, E. Casali, L. Mason and S. Nekovar, _Scattering on plane waves and the double copy_ , _Classical and Quantum Gravity_ 35 (2017) 015004.
  * [51] R. Penrose, _Zero rest mass fields including gravitation: Asymptotic behavior_ , _Proc. Roy. Soc. Lond. A_ 284 (1965) 159.
  * [52] L.J. Mason, _On ward’s integral formula for the wave equation in plane wave space-times_ , _Twistor Newsl. 28 17–9, 1989_ .
  * [53] R. Penrose and W. Rindler, _SPINORS AND SPACE-TIME. VOL. 2: SPINOR AND TWISTOR METHODS IN SPACE-TIME GEOMETRY_ , Cambridge Monographs on Mathematical Physics, Cambridge University Press (4, 1988), 10.1017/CBO9780511524486.
  * [54] S. Huggett and K. Tod, _AN INTRODUCTION TO TWISTOR THEORY_ , London Mathematical Society Student Texts (9, 1986).
  * [55] R. Penrose and W. Rindler, _Spinors and Space-Time_ , Cambridge Monographs on Mathematical Physics, Cambridge Univ. Press, Cambridge, UK (4, 2011), 10.1017/CBO9780511564048.
  * [56] R.S. Ward and R.O. Wells, _Twistor geometry and field theory_ , Cambridge Monographs on Mathematical Physics, Cambridge University Press (8, 1991), 10.1017/CBO9780511524493.
  * [57] T. Adamo, _Lectures on twistor theory_ , _PoS_ Modave2017 (2018) 003 [1712.02196].
  * [58] R.P. Kerr and W.B. Wilson, _Singularities in the Kerr-Schild metrics._ , _General Relativity and Gravitation_ 10 (1979) 273.
  * [59] P. Tod, _Some geometry of de Sitter space_ , 1505.06123.
  * [60] E.T. Newman, _Asymptotic twistor theory and the Kerr theorem_ , _Class. Quant. Grav._ 23 (2006) 3385 [gr-qc/0512079].