# Data-driven Heuristics for DC optimal transmission switching problem

Juncheng Li Email<EMAIL_ADDRESS>Department of Management Science,
Lancaster University, United Kingdom  Trivikram Dokka Email<EMAIL_ADDRESS>Department of Management, Queens University Belfast, United Kingdom  Guglielmo
Lulli Email<EMAIL_ADDRESS>Department of Management Science,
Lancaster University, United Kingdom  Fabrizio Lacalandra Email:
<EMAIL_ADDRESS>The Italian Regulatory Authority for Energy, Networks and
Environment (ARERA), Italy

###### Abstract

The goal of Optimal Transmission Switching (OTS) problem for power systems is
to identify a topology of the power grid that minimizes the cost of the system
operation while satisfying the operational and physical constraints. Among the
most popular methods to solve OTS is to construct approximation via integer
linear programming formulations, which often come with big-M inequalities.
These big-M inequalities increase, considerably, the difficulty of solving the
resulting formulations. Moreover, choosing big-M values optimally is as hard
as solving OTS itself. In this paper, we devise two data-driven big-M bound
strengthening methods which take network structure, power demands and
generation costs into account. We illustrate the robustness of our methods to
load changes and impressive runtime improvements of mixed-integer solvers
achieved by our methods with extensive experiments on benchmark instances. The
speedup by one of the proposed methods is almost 13 times with respect to the
exact method.

##  1 Introduction

The economical dispatch of electric power on grids has been attracting the
attention of both academia and Operations Research (OR) practitioners in the
energy sector for decades. The optimal power flow (OPF) problem is to
determine a minimum-cost delivery scheme of electric power to meet power
demand subject to physical laws - such as Kirchhoff’s Laws and Ohm’s Laws - as
well as other operational limits imposed by the grid. The OPF problem is a
non-convex quadratic constrained quadratic programming problem first proposed
by Carpentier [4] and studied since then. While the OPF problem treats the
configuration of the grid as fixed, it is now a common practice to temporarily
switch off some branches in specific operating hours [13].

Economic benefits from switching off branches of a grid can be tremendous, as
it enables system operators to lower significantly the overall production
costs and/or balance the load on the grid [8, 12, 14]. Indeed, the optimal
unit commitment schedule can change by modifying the grid topology [11].
Equally important, is that the use of switching in grids may facilitate the
integration of renewable power into the system [22]. The problem of finding
the optimal grid topology, i.e., the topology that allows meeting the energy
demand at the lowest generation costs, is known in the literature as Optimal
Transmission Switching (OTS). The name derives from the fact that the optimal
topology is obtained by switching off some branches. From the mathematical
modelling point of view, the OTS problem is formulated by adding the switching
component of branches to the OPF model. This leads to a mixed integer non-
convex nonlinear programming formulation that poses an extraordinary
computational challenge. Under some realistic operational assumptions, the OTS
formulation is often approximated as a mixed integer linear program (MILP)
with big-M parameters [8], to take advantage of mature MIP solvers. Such a
MILP approximation - which is still NP-hard [18][17] \- is referred in the
literature as the DC OTS problem. Given the computational complexity of the
problem, a number of heuristic methods have been proposed in the literature.
Fuller et al. [9] presented two greedy algorithms that iteratively switch off
branches guided by linear programming duality information. Ruiz et al. [21]
analyzed three different policies to remove braches in a greedy type
algorithm. Prescreening heuristics to reduce the set of switchable arcs
according to properties of the electric power system were proposed in [2, 19,
24]. Johnson et al. [16] studied a k-nearest-neighbor clustering approach to
compute DC OTS solution making use of historical information. On the other
side of the methodological spectrum, i.e, exact methods, Kocuk et al.[17] gave
a polyhedral study of the problem and proposed cycle-based valid inequalities
to strengthen the MILP formulation.

In this work, we propose to compute solutions of the OTS problem using MIP
solvers. However, to exploit the capabilities of these solvers, the values of
the big M parameters have to be carefully set; and this is the objective of
this paper. In fact, while small big-M parameters result in strong linear
relaxation that speeds up solution procedure in MIP solvers, undersized ones
can cut off the optimal solution. Binato et al. [3] observed that elementary
weighted longest paths between the endpoints of each branch of the grid can be
used to define the big-M parameters for the grid expansion planning problem -
a problem that shares similar characteristics with DC OTS. Indeed, we show in
§2 that big-M parameters determined by elementary weighted longest paths
reserves at least one optimal solution of the DC OTS problem. However,
longest-path big-M parameters may be too conservative (big) because they
preserve all the feasible region without taking any power demand and
production cost into account. This motivates the development of methods to
compute smaller big-M parameters. However, the problem of computing the
smallest big-M parameters that do not cut off any feasible solution of the DC
OTS problem is NP-hard and it is as difficult as solving the DC OTS problem,
see [6]. This prompts the development of heuristic methods to compute big-M
parameters for DC OTS problem, which is absent in the existing literature. In
this work, we propose two polynomial-time heuristic methods to compute big-M
parameters for DC OTS problem by solving a sequence of linear programming
problems. The rest of the paper is organized as follows: §2 provides a brief
description of DC OTS problem as well as an example motivating the use of
input data to compute small big-M values. §3 illustrates the 2 proposed
heuristic methods to compute big-M values. §4 shows computational experiments
comparing the proposed heuristic methods and the exact method. §5 draws the
conclusion of this work and mentions possible future research.

##  2 Big-M values and motivating example

For the sake of completeness, we first present the mathematical programming
formulation (2.1) of the DC OTS problem. The interested reader may refer to
[8] for further details. The grid is represented by a graph $D=(N,E)$, where
$N$ denotes the set of buses (nodes) and $E$ denotes the set of branches
(edges) linking the buses. $G\subset N$ is the set of buses connected to power
generators. Let the positive $f_{ij}$ denote the power flow from bus $i$ to
bus $j$. The negative value of $f_{ij}$ denotes the power flow in the opposite
direction (from bus $j$ to bus $i$).

$\displaystyle\min_{p,f,\theta,x}\quad$ $\displaystyle\sum_{i\in G}p_{i}\cdot
c_{i}$ (2.1a) s.t. $\displaystyle\sum_{j\in V}f_{ji}-\sum_{j\in
V}f_{ij}=d_{i},\;\forall i\in N\setminus G$ (2.1b) $\displaystyle\sum_{j\in
V}f_{ji}-\sum_{j\in V}f_{ij}+p_{i}=d_{i},\;\forall i\in G$ (2.1c)
$\displaystyle f_{ij}=B_{ij}(\theta_{i}-\theta_{j})\cdot
x_{ij},\;\forall(i,j)\in E$ (2.1d) $\displaystyle-\bar{f_{ij}}\cdot x_{ij}\leq
f_{ij}\leq\bar{f_{ij}}\cdot x_{ij},\;\forall(i,j)\in E$ (2.1e) $\displaystyle
p_{i}^{min}\leq p_{i}\leq p_{i}^{max},\;\forall i\in G$ (2.1f) $\displaystyle
x_{ij}\in\\{0,1\\},\;\forall(i,j)\in E.$ (2.1g)

The objective function (2.1a) minimizes the power production costs.
Constraints (2.1b) and (2.1c) ensure that net power flow into each bus equals
its power demand ($d_{i}$), where $p_{i}$ is the power provided by generator
i. (2.1d) imposes that power flow across branch $(i,j)$ equals the voltage
angle difference ($\theta_{i}-\theta_{j}$) multiplied by its susceptance
($B_{ij}$), if branch $(i,j)$ is on service (i.e., $x_{ij}=1$. $x_{ij}$ is
binary variable representing on/off status of the branch). Constraints (2.1e)
and (2.1f) impose limits on power flow ($\bar{f_{ij}}$) on each branch and on
power production ($p_{i}^{min},p_{i}^{max}$) at each bus i of $G$. The
bilinear constraint (2.1d) can be linearized with big-M parameters using the
following inequalities:

$-(1-x_{ij})M_{ij}\leq
f_{ij}-B_{ij}(\theta_{i}-\theta_{j})\leq(1-x_{ij})M_{ij},\\\ \forall(i,j)\in
E$ (2.2)

$M_{ij}$ is the big M value, i.e., a "large enough" number such that
constraint (2.2) is redundant when $x_{ij}=0$.(i.e. implied by other
constraints), thus equivalent to non-existent. We here show that $M_{ij}$
computed by elementary weighted longest path from node $i$ to node $j$ does
not cut off any feasible solution under realistic assumption.

######  Proposition 1.

[17] In the DC OTS problem, there exists an optimal solution in which the
lines switched on form a connected network.

It is also desirable for electric grid to connected, see [15]. In the
remaining of this paper, we assume solution of DC OTS problem is always
connected.

######  Proposition 2.

In the DC OTS problem, if a elementary path $P$ from node $i$ to node $j$ is
on service, then weight of $P$ times $B_{ij}$ is a valid $M_{ij}$.

###### Proof.

As the path $P$ is on service, $x_{uv}=1$ for every edge $(u,v)\in P$. Thus,
for each edge $(u,v)\in P$, $-\bar{f}_{uv}\leq
B_{uv}(\theta_{u}-\theta_{v})\leq\bar{f}_{uv}$. Therefore, $-\sum_{(u,v)\in
P}\frac{\bar{f}_{uv}}{B_{uv}}\leq\sum_{(u,v)\in
P}(\theta_{u}-\theta_{v})\leq\sum_{(u,v)\in P}\frac{\bar{f}_{uv}}{B_{uv}}$. As
$\theta_{i}-\theta_{j}=\sum_{(u,v)\in P}(\theta_{u}-\theta_{v})$, this implies
that $-B_{ij}\sum_{(u,v)\in P}\frac{\bar{f}_{uv}}{B_{uv}}\leq
B_{ij}(\theta_{i}-\theta_{j})\leq B_{ij}\sum_{(u,v)\in
P}\frac{\bar{f}_{uv}}{B_{uv}}$. Thus, constraint (2.8) is redundant when
$x_{ij}=0$, with $M_{ij}:=B_{ij}\sum_{(u,v)\in P}\frac{\bar{f}_{uv}}{B_{uv}}$.
∎

Proposition 2 shows that $M_{ij}$ equalling $B_{ij}$ times the weight of the
longest elementary path from node $i$ to node $j$ does not cut off any
feasible solution, where weight of a path $P$ is defined as $\sum_{(u,v)\in
P}\frac{\bar{f}_{uv}}{B_{uv}}$. We refer to this big M value as longest-path
big-M value.

The choice of $M_{ij}$ parameters influences significantly the performance of
MIP solvers based on LP branch-and-cut procedure. Large values of $M_{ij}$’s
often result in “weak" LP relaxations thus hampering the performance of MIP
solvers. Indeed, in the formulation it is desirable to have the smallest big M
values that do not cut off the optimal solution. Fattahi et al. [6] define the
theoretically smallest $M_{ij}$ as

$M_{ij}^{opt}=B_{ij}\times\max_{(p,f,\theta,x)\in
F}\\{|\theta_{i}-\theta_{j}|\\}$ (2.3)

where
$F:=\\{(p,f,\theta,x):\eqref{eq:DC_OTS_power_balance}-\eqref{eq:x_bin}\\}$.
They also showed that i) computing $M_{ij}^{opt}$ is NP-hard and there is no
polynomial-time constant-factor approximation algorithm for $M_{ij}^{opt}$;
and ii) $M_{ij}^{opt}$ can be arbitrarily large on networks with special
topology. But even using the weighted longest path between the endpoints of
each edge of the network as big M value, as suggested by Binato et al. [3]
among others, the computation of such values is still demanding at least from
the theoretical point of view. Moreover, the weighted longest past is based
exclusively only on the network topology and ignores completely the power
demand at each bus and power production costs at each power plant, which may
lead to too conservative values (larger than $M_{ij}^{opt}$). Power demand,
generation capacity and production costs are all relevant information that can
be used to compute smaller big M values, as shown in the example below.

Consider the example depicted in Figure 1. Suppose that the links marked with
a thicker line are switchable, while the others are not. The number next to
the link is the transmission capacity $\bar{f}$ and the susceptance $B$ is 1
for all the links. In this case, the big-M parameters suggested by (weighted)
longest path are equal to 6 for both the switchable links. Suppose now that
the generator $g1$ is cheaper and there is a power demand of 2.5 MW at bus t
of the network. The power demand at the other buses is null. In this case, a
big-M value of 1 for link $(f,g2)$ is sufficient. To see this note that in the
event $(b,g1)$ is switched on a maximum of 1.5 units can be supplied from
$g1$, this implies a maximum of 1 units is supplied from $g2$. In fact even
$M_{f,g2}=0.5$ is sufficient as no unit is required from $g2$ if $(b,g1)$ is
switched off. This example suggests taking network structures, power demand
and generation costs into consideration to get smaller big M parameter.

$t$$f$$d$$b$$h$$g1$$g2$$3$$3$$3$$3$$3$$3$$1$$1$ Figure 1: Big-M Example
(Switchable links - highlighted)

##  3 Heuristic methods for big M values

In this section, we present two data driven polynomial-time heuristic methods
to compute big M parameters. Both the methods take power demands and
production costs into account. We call the first method k-shortest-path (kSP)
and the second k-nearest-neighborhood simulation method (kNN).

###  3.1 k-shortest-path method

The k-shortest-path method is based on the observation that follows directly
from Proposition 2. Indeed, the phase angle difference between any pair of
nodes -and hence the big M value - is limited above by the weighted shortest
path between the two nodes, where the weight of each edge is given by
$\frac{\bar{f}}{B}$. If one edge of the shortest path is removed from the
graph (i.e., it is switched off), than the second shortest path becomes the
one limiting the phase angle difference between the two nodes. We can iterate
the procedure, i.e., removing one edge from the current shortest path and
selecting the next shortest path, until we either reach the k-shortest-path
(for a given parameter k) or we obtain an infeasible problem, meaning that no
edge of the path can be switched off. In other words, this path has to be on
service, and its weight determines the big M value. Given a path, we use the
approach proposed by Fuller et al. [9] to select the edge to be switched off.
Fueller et al. used DC OPF duality information and KKT optimality conditions
to compute edges’ shadow price, which provides a (local) measure of objective
function’s improvement if the edge is switched off. More in particular, they
showed that the shadow price $\alpha_{ij}$ of edge (i,j) is given by
$(\pi_{i}-\pi_{j})f_{ij}$, where $\pi_{i}$ and $\pi_{j}$ are the dual variable
of the power balance constraint at buses i and j of the OPF problem
respectively. These constraints are exactly the same as constraints (2.1b) and
(2.1c) of the DC OTS formulation described above.

A negative value of $\alpha$ "infers" that removing edge $(i,j)$ from the
network would results in lower production costs. Therefore, for each shortest
path, we iteratively check whether an edge can be off service without making
DC OPF infeasible from the smallest $\alpha$ to the largest one. We also test
edges with $\alpha>0$, because we are only interested in verifying if the path
is a be-on-service  path or not, i.e., if the underlying OPF problem - after
disconnecting the path - is feasible or not.

Once we find that an edge can be switched off, we remove this edge from the
graph and proceed to the next shortest path. If none of the edges in the path
can be switched off (i.e. removing any of those edges makes the reduced OPF
problem infeasible), we attempt to restore feasibility of the reduced OPF
problem by removing one more edge outside the path that contributes to
infeasibility. If feasibility of the reduced OPF problem is successfully
retrieved in this way, it shows that path under consideration can be off-
service and the algorithm proceed to consider the next shortest path.
Otherwise, we label the path under consideration as "must-on-service" and
outputs the big M value determined by weight of this path. In principle, we
can attempt to restore feasibility of the infeasible reduced DC OPF problem by
removing more than one edge. However, that leads to significantly heavier
computation burden thus we choose not to do so.

$\displaystyle\min_{p,f,\theta}\quad$ $\displaystyle c^{T}\cdot p$ (3.1a) s.t.
$\displaystyle\sum_{j\in V}f_{ji}-\sum_{j\in V}f_{ij}+p_{i}=d_{i},\;\forall
i\in N$ (3.1b) $\displaystyle
f_{ij}=B_{ij}(\theta_{i}-\theta_{j}),\;\forall(i,j)\in E$ (3.1c)
$\displaystyle-\bar{f_{ij}}\leq f_{ij}\leq\bar{f_{ij}},\;\forall(i,j)\in E$
(3.1d) $\displaystyle p_{i}^{min}\leq p_{i}\leq p_{i}^{max},\;\forall i\in G$
(3.1e)

$\displaystyle\min_{p,f,\theta}\quad$ $\displaystyle c^{T}\cdot p$ (3.2a) s.t.
$\displaystyle\sum_{j\in V}f_{ji}-\sum_{j\in V}f_{ij}+p_{i}=d_{i},\;$
$\displaystyle\forall i\in N$ $\displaystyle\qquad[\pi_{i}]$ (3.2b)
$\displaystyle f_{ij}=B_{ij}\cdot x_{ij}\cdot(\theta_{i}-\theta_{j}),\;$
$\displaystyle\forall(i,j)\in E$ $\displaystyle\qquad[\sigma_{ij}]$ (3.2c)
$\displaystyle 1-x_{ij}=0,\;$ $\displaystyle\forall(i,j)\in E$
$\displaystyle\qquad[\alpha_{ij}]$ (3.2d) $\displaystyle-\bar{f_{ij}}\cdot
x_{ij}\leq f_{ij}\leq\bar{f_{ij}}\cdot x_{ij},\;$
$\displaystyle\forall(i,j)\in E$ (3.2e) $\displaystyle p_{i}^{min}\leq
p_{i}\leq p_{i}^{max},\;$ $\displaystyle\forall i\in G$ (3.2f)

In order to find the edges contributing to infeasibility, we use a method
inspired by the method of identifying minimal (or irreducible) infeasible
subsystem (MIS or IIS) of linear inequalities in [5]. Consider an arbitrary
infeasible set of linear equalities and linear inequalities

$\\{x:Ax=b,Hx\geq g\\}$ (3.3)

A minimal infeasible subsystem of (3.3) is a inclusion-minimal set of linear
equalities and linear inequalities selected from $Ax=b$ and $Hx\geq g$,
respectively:

$\\{x:A^{\prime}x=b^{\prime},H^{\prime}x\geq g^{\prime}\\}$ (3.4)

such that it has no feasible solution (i.e. (3.4) is empty). In order to find
a feasible subsystem of (3.3), we must remove at least one linear equality or
inequality from each of its MIS.

In the DC OPF problem (3.1), removing an edge $(i,j)$ from the network is
equivalent to removing the constraint (3.1c) corresponding to the edge $(i,j)$
and imposing an additional constraint $f_{ij}=0$. Let us consider the
constraints (3.1c) in a MIS of the OPF problem (3.1) and the set of edges
$E_{MIS}$ corresponding to those constraints. OPF problem (3.1) on switchable
network can be feasible only if at least one edge in $E_{MIS}$ is switched
off. Thus, we restrict the candidates for $e_{out}$ to $E_{MIS}$. Next we
describe implementation details of this method. Let us use $n$ to represent
the number of nodes, $m$ to represent the number of edges, $v_{n}$ to
represent a vector $v$ with $n$ entries and $V_{n\times m}$ to represent a
matrix $V$ with $n\times m$ entries. For the ease of presentation, in the
remaing of this work, we will present the OPF problem (3.1) in a compact
matrix-vector form
$\min_{x}\\{\hat{c}^{T}x:\hat{A}x=\hat{b},\hat{H}x\geq\hat{g}\\}$, where
$\hat{c}=\begin{pmatrix}c_{n}\\\ 0_{n}\\\ 0_{m}\end{pmatrix}$ ,
$x=\begin{pmatrix}p_{n}\\\ \theta_{n}\\\ f_{m}\end{pmatrix}$,
$\hat{A}=\begin{pmatrix}I_{n\times n}&0_{n\times n}&\Omega_{n\times m}\\\
0_{m\times n}&B_{m\times n}&-I_{m\times m}\end{pmatrix}$,
$\hat{H}=\begin{pmatrix}0_{m\times n}&0_{m\times n}&I_{m\times m}\\\
0_{m\times n}&0_{m\times n}&-I_{m\times m}\\\ I_{n\times n}&0_{n\times
n}&0_{n\times m}\\\ -I_{n\times n}&0_{n\times n}&0_{n\times m}\end{pmatrix}$,
$\hat{b}=\begin{pmatrix}d_{n}\\\ 0_{m}\end{pmatrix}$ and
$\hat{g}=\begin{pmatrix}-\bar{f}_{m}\\\ -\bar{f}_{m}\\\ \underline{p}_{n}\\\
-\bar{p}_{n}\end{pmatrix}$. $I_{n\times n}$ means $n$-dimensional identity
matrix, $\Omega_{n\times m}$ denotes the node-arc incidence matrix, matrix
$B_{m\times n}$ is a sparse matrix with non-zero components: $B_{l,i}=B_{ij}$
and $B_{l,j}=-B_{ij}$ for each edge $l$ starting from node $i$ and ending at
node $j$. To identify $e_{out}$, we first remove the edge with the smallest
$\alpha$ value in the path from the network. Let us denote the remaining
constraints in (3.1) as $\\{\hat{A}x=\hat{b},\hat{H}x\geq\hat{g}\\}$ and
construct the LP $\min_{x}\\{0^{T}x:\hat{A}x=\hat{b},\hat{H}x\geq\hat{g}\\}$
as well as its dual
$\max_{y,u}\\{y^{T}\hat{g}+u^{T}\hat{b}:y^{T}\hat{H}+u^{T}\hat{A}=0^{T},y\geq
0\\}$ with $u$ being the dual variables of equality constraints and $y$ being
the dual variables of inequality constraints. Since dual infeasibility is
excluded by the solution $(y,u)=0$, primal infeasibility corresponds to dual
unboundedness, i.e. there exists a dual solution $(y,u)$ such that
$y^{T}\hat{g}+u^{T}\hat{b}>0$. Therefore, we can normalize the dual solution
by replacing the objective with a constraint $y^{T}\hat{g}+u^{T}\hat{b}=1$. It
is shown in [10] that the non-zero elements of each vertex of the polyhedron

$\\{(y,u):y^{T}\hat{H}+u^{T}\hat{A}=0^{T},y^{T}\hat{g}+u^{T}\hat{b}=1,y\geq
0\\}$ (3.5)

defines a MIS of (3.3) and thus a candidate set for $e_{out}$. In order to
find a vertex of (3.5), we minimize $1^{T}y$ over (3.5). This LP is guaranteed
to be solved to optimality. We then use the set of edges corresponds to non-
zero $u_{i},i\in J$ in the optimal solution as the candidates for $e_{out}$,
where $J$ is the set of indices of constraints (3.1c). We then attempt to
select an edge from the candidates to be $e_{out}$ as shown in Algorithm 1. If
no edge in the candidate set satisfies the requirement of $e_{out}$, we decide
that the current path "must" be on service in an economical network
configuration, which determines the k-shortest-path big-M values. A formal
description of kSP method is given in Algorithm 1, where "dualMIS" denotes the
problem of minimizing $1^{T}y$ over (3.5) and $w(P^{k}_{ij})$ denote weight of
$k$-th shortest path from node $i$ to node $j$.

Algorithm 1 \- kSP

1:Set parameters: $k_{max}$, $e_{max}$, $l$

2:for $(i,j)\in E$ do

3: Solve DC OPF on G

4: $\alpha^{0}\leftarrow\alpha^{0}(G)$

5: $k\leftarrow 1$

6: while $k<k_{max}$ do

7: if all edges in $k$-th shortest path is on service then

8: Order edges in $k$-th shortest path from smallest $\alpha^{k-1}$ to
largest: $\\{j_{1},\ldots,j_{m}\\}$

9: $i\leftarrow 0$

10: while $i<e_{max}$ do

11: $i\leftarrow i+1$

12: Solve DC OPF on $G^{\prime}=G(N;E\setminus j_{i})$

13: if DC OPF feasible then

14: $\alpha^{k}\leftarrow\alpha^{k}(G^{\prime})$

15: $E\leftarrow E\setminus j_{i}$

16: Go to line 34

17: end if

18: end while

19: Solve dualMIS on $G^{\prime}=G(N;E\setminus j_{1})$

20: Compute $E^{\prime}=\\{(p,q)\in E:u_{pq}\neq 0\\}$

21: Compute $E^{\prime\prime}=\\{h_{1},\ldots,h_{s}\\}$ by ordering edges in
$E^{\prime}$ from smallest $\alpha^{k-1}$ to largest

22: $v\leftarrow 0$

23: while $v<s$ do

24: $v\leftarrow v+1$

25: Solve DC OPF on $\hat{G}=G(N;E\setminus\\{j_{1},h_{v}\\})$

26: if DC OPF feasible then

27: $\alpha^{k}\leftarrow\alpha^{k}(\hat{G})$

28: $E\leftarrow E\setminus\\{j_{1},h_{v}\\}$

29: Go to line 34

30: end if

31: end while

32: Go to line 36

33: end if

34: $k\leftarrow k+1$

35: end while

36: $k\leftarrow k+l$

37: $M_{ij}=w(P_{ij}^{k})\times B_{ij}$

38:end for

###  3.2 k-nearest-neighbor Simulation Method

Compared to k-SP, our second heuristic is a randomized method that relies on
the following intuition. As illustrated in Figure 1, the $M_{ij}$ value
depends on which edges remain in the network. More importantly, the value of
$M_{ij}$ is likely to depend on the status of switchable edges that are in
close vicinity rather than those which are spatially far. This is because
switching off nearby edges likely removes shorter paths forcing a larger
$M_{ij}$ value as compared to faraway edges. We use this observation to
explore the impact on $M_{ij}$ by switching off randomly selected edges in a
given neighbourhood. Moreover, this heuristic is inspired by large class of
randomized local search heuristics often used in large scale combinatorial
optimization such as GRASP heuristics [7]. A formal description of the
k-nearest-neighbor (kNN) simulation method is given below in Algorithm 2.

Algorithm 2 \- kNN

1:Set parameters: $k$, $h<1$, $s$, $r$

2:for $(i,j)\in E$ do

3: $inf\leftarrow 0$

4: Compute $N_{ij}=\\{(u,v)\in E:d_{ij}(u)\leq k\land d_{ij}(v)\leq k\\}$

5: for $iter=1,2,\ldots,r$ do

6: Select $N^{\prime}_{ij}\subset N_{ij}:|N^{\prime}_{ij}|=h\cdot|N_{ij}|$

7: Solve DC OPF on $G^{\prime}=G(N;E\setminus N^{\prime}_{ij})$

8: if DC OPF feasible then

9: $\Delta\Theta^{iter}_{ij}=|\theta_{i}-\theta_{j}|$

10: else$\hskip 5.69054pt\Delta\Theta^{iter}_{ij}=0$ and $inf\leftarrow inf+1$

11: end if

12: end for

13: $\Delta\Theta_{ij}\leftarrow max_{iter}\\{\Delta\Theta^{iter}_{ij}\\}$

14: if $inf<r$ then

15: $M_{ij}=\min\\{s\times\Delta\Theta_{ij},w(P_{ij}^{L})\\}\times B_{ij}$

16: else$\hskip 5.69054ptM_{ij}=w(P_{ij}^{L})\times B_{ij}$

17: end if

18:end for

At line 4 of the algorithm, to compute the neighborhood of edge $(i,j)$
($N_{ij}$), i.e., the set of edges that are close to $(i,j)$, we use Dijktra’s
algorithm. Indeed, we consider in the neighborhood all the edges whose
endpoints are at most $k$ edges away from either i or j (in formula,
$d_{ij}(u)\leq k$). If the DC OPF problems are all infeasible - line 7 of the
algorithms -, we use as big M value for edge $(i,j)$ the value provided by the
longest past method $w(P_{ij}^{L})\times B_{ij}$; wehre $w(P_{ij}^{L})$ is the
weight of longest weighted path from node $i$ to node $j$ ($P_{ij}^{L}$).

###  3.3 An illustrative example

In this section we use a small benchmark instance IEEE 14-bus system [1]
displayed in Figure 2, to illustrate with an example on how longest-path,
k-shortest-path and kNN simulation method compute big-M value. For this
purpose, we show the process of computing big M value of edge $(3,4)$ -i.e.,
$M_{3,4}$ \- with the listed three methods. Transfer capacity is set to 27 for
all the branches, as these values were not provided in the original data. For
kSP and kNN simulation methods, we consider the following parameters’ setting:
$k_{max}=5$, $e_{max}=3$ and $l=1$ for kSP and $k=2,h=20,s=10$, and $r=30$ for
kNN. The big-M value of edge $(3,4)$ returned by the three methods is 355.5
(longest-path method), 260.9 (kSP method) and 216.4 (kNN simulation method).

In the longest-path method, $M_{3,4}$ is given by $B_{34}$ times the weight of
the (weighted) longest path, i.e., path $3\to 2\to 1\to 5\to 6\to 12\to 13\to
14\to 9\to 4$.

As far as the kSP method, the first shortest path to consider (first
iteration) is edge $(3,4)$ itself. We remove edge $(3,4)$ from the network and
solve the OPF problem $OPF((3,4))$, which turns out to be infeasible. We
attempt to restore the feasibility of the OPF problem by removing one
additional edge $(e_{out})$ from the network. To identify candidate edges for
$e_{out}$, we solve $dualMIS((3,4))$. We rank this set of edges by linear
programming duality as described in §2:
$\\{(4,5),(2,4),(6,11),(1,2),\ldots\\}$. We check if the OPF problem becomes
feasible by removing one additional edge of the list. The first feasible OPF
problem is $OPF((3,4),(1,2))$, with $e_{out}:=(1,2)$. We then remove $(3,4)$
and $(1,2)$ from the network. At the second iteration, the shortest path (on
the current network) is $3\to 2\to 4$. Using arguments of §2, we rank the two
edges as follows: $\\{(4,2),(2,3)\\}$. We solve $OPF((4,2))$, which turns out
to be infeasible. We then solve $OPF((2,3))$, which is feasible. We then
remove edge $(2,3)$ from the network. At the third iteration, it turns out
that the third, fourth and fifth shortest paths are no longer active as some
edges of those paths have already been removed. As $k_{max}$ is reached when
considering the fifth shortest path, we terminate the algorithm and output
$k_{34}=5+l=6$. Big-M value on edge $(3,4)$ is defined by $B_{34}$ times
weight of the sixth shortest path from node $3$ to node $4$.

In the kNN simulation method, first we remove edge $(3,4)$ from the network
and identify the set of nodes whose distance is at most $k=2$ from either node
$3$ or node $4$, i.e., $\\{1,2,5,7,9,10,14,8\\}$. All the edges with both the
endpoints in the list define the neighborhood, specifically
$\\{(4,2),(3,2),(2,1),(2,5),$
$(1,5),(4,5),(9,10),(4,9),(7,9),(4,7),(7,8),(9,14)\\}$. Select randomly
$h=20\%$ of the edges (i.e., 2 edges in this example) from the neighborhood.
Solve OPF problem on the network without these randomly selected (two) edges
and compute the phase angle difference
$\Delta\theta_{34}=|\theta_{4}-\theta_{3}|$ from the OPF solution. Repeat this
procedure for $r=30$ times and denote the largest $\Delta\theta_{43}$ as
$\Delta^{\prime}\theta_{43}$. We set as big M value for edge $(4,3)$
$M_{3,4}=\min\\{s\times\Delta\theta_{43}^{\prime},w(P_{43}^{L})\\}\times
B_{43}$ , where $P_{43}^{L}$ denote the longest weighted path from node 3 to
node 4.

(a) diagram of 14-bus system

(b) network of 14-bus system

Figure 2: IEEE 14 bus instance

##  4 Computational Experiments

In this section, we present an extensive computational study that demonstrates
the effectiveness of our proposed approaches. We compare the kSP and kNN
methods with the longest-path approach (LWP) on a set of benchmark instances.
More specifically, we consider two IEEE test cases, i.e., case118B and
case300B, widely used in the literature. The former consists of 118 buses, 180
branches and 19 committed generators with different power production costs.
The latter includes 300 buses, 409 branches and 61 committed generators with
different power production costs. In addition to the nominal demand scenario,
which is part of the IEEE instances, we also consider both a low-demand and a
high-demand scenario. These demand scenarios are obtained from the nominal one
multiplying the load at each bus by 0.95 and 1.05 respectively. For each pair
of grid and demand scenario - defining a class of instances -, we generated
twenty instances by multiplying the load at each bus by a random variable with
a uniform distribution in the interval [0.95, 1.05]. Hence, in total we
considered 120 instances of the OTS problem with different characteristics.
"118/300L", "118/300N", "118/300H" denote the low, normal and high demand
class of instances, respectively.
All the computational experiments are run on a laptop with Intel(R) Core(TM)
i7-8750H @2.20 GHz 2.21GHz processor and 16GB RAM. We use Julia v1.5.2 with
JuMP v0.21.5 programming language [20] to implement the model and the
algorithms, and Gurobi 9.0.1 as a solver. To compute the big M values with the
longest weighted path approach, the Miller-Tucker-Zemlin (MTZ) longest path
formulation [23] is implemented and solved by Gurobi 9.0.1, imposing a time
limit of 600 seconds and setting an optimality gap to 1%. k-shortest paths are
computed using the algorithm proposed in [25].

###  4.1 Parameters’ tuning and results

The first phase of the computational study has been devoted to tuning the
parameters of both the kSP and kNN methods. The goal is to identify a set of
parameters which work well across all instances of a given class.

Table 1 lists all the values of the parameters we have tested for each class
of instances (CoI). The “best" parameters’ setting for each class of instances
are highlighted with bold fonts.

| kSP | kNN
---|---|---
CoI | $k_{max}$ | $e_{max}$ | $l$ | $k$ | $h$ | $s$ | $r$
118L | 11,14,17,20 | 3,5,7 | 0,1 | 3,4,5 | 5, 10 | 5,10,15 | 50
118N |  11,14,17,20 | 3,5,7 | 0,1 | 3,4,5 | 5, 10 | 5,10,15 | 50
118H |  11,14,17,20 | 3,5,7 | 0,1 | 3,4,5 |  5, 10 | 5,10,15 | 50
300L | 15,20,25,30 | 5,7,9 |  0,1 | 4,5,6,7 |  5, 10 | 5,10,15 | 50
300N |  15,20,25,30 |  5,7,9 | 0,1,2,3 |  4,5,6,7 |  5, 10 | 5,10,15 | 50
300H |  15,20,25,30 | 5,7,9 | 0,1,2,3,4,5 | 4,5,6,7 | 5, 10 | 5,10,15 | 50
Table 1: Parameters’ values tested for kSP and kNN.

Given the nature of our approach, i.e., compute small enough big M values that
allow computing (near) optimal OTS solutions in short computational time,
quality of the solution and computational time are the two measures of
interest. Therefore, we use both the measures to assess the effectiveness of
proposed methods. To measure the solutions’ quality, we use the following
ratio $\frac{z_{m}-z_{LP}}{z_{LP}}$. This ratio - relative gap - compares the
objective value of the solution computed by the proposed method ($z_{m}$) with
the objective value of the solution computed with big M values obtained by the
longest-path method ($z_{LP}$). The smaller the relative gap, the better the
solution quality. Negative values of the relative gap indicate that the
proposed approach provides better solutions than the longest-path method. It
is important to recall that all the instances are solved by the solver with a
time limit of 600 sec. Therefore, for the sake of accuracy, we are comparing
the best solution computed within the imposed time limit, which is the optimal
one in many instances. Figure 3 display the scatter plots of the kSP method
for each class of instances. Each point of the scatter plot represents the
average performance of kSP with one setting of the parameters on the twenty
instances of the class. It appears that parameter l is key to the good
performance of the method as it influences the trade off between quality of
solution and computational time. Indeed this parameter was introduced to
ensure a certain degree of robustness of the proposed approach by enlarging
the big M values. To larger values of l correspond better solutions’ quality
(smaller relative gaps) at the cost of longer computational times. For the
instances derived from case118B IEEE, $l=1$ provides a very good balance
between solutions’ quality and computational time. The relative gaps are very
small - negative in several cases - independently of the demand scenario and
the setting of the other two parameters. The computational times are also
short, especially in the low and high demand scenarios. However, in this case
the choice of $k_{max}$ and $e_{max}$ settings is relevant for the time. The
settings that guarantee the shortest computational times are the chosen ones
or the "best". It is also important to observe, that these settings are quite
robust. In particular, the triplet $k_{max}=11$, $e_{max}=5$ and $l=1$ that is
best for the nominal and high demand scenario performs quite well also in the
low demand scenario (marked on the scatter plots with a star).

Similar observations are also valid for the case300B IEEE family of instances.
However, the class of instances with low demand scenario (300L) turn out to be
“easy" to solve and all the tested parameters perform reasonably well (see
scatter plot in Figure 3(d)). For the nominal and high demand scenarios, the
effect of parameter $l$ becomes evident. Indeed, we observe again data in
Figures 3(e) and 3(f) are clustered by values of $l$.

As far as the kNN method, see Figure 4, it computes good solutions of case118B
instances. Indeed, for all the demand scenarios, it computes solutions with
(very) small average gap though the computational times are not always very
competitive at least in comparison with the kSP method. On the case300B test
instances, though the performance with specific settings of the parameters can
be good, it appears that several of the tested settings do not provide
competitive solutions demonstrating the greater challenge of identifying the
right settings to use.

(a) data set: 118L

(b) data set: 118N

(c) data set: 118H

(d) data set: 300L

(e) data set: 300N

(f) data set: 300H

Figure 3: parameter tuning for kSP method

(a) data set: 118L

(b) data set: 118N

(c) data set: 118H

(d) data set: 300L

(e) data set: 300N

(f) data set: 300H

Figure 4: parameter tuning for kNN simulation method

###  4.2 Computational results

In this section, we report the computational results of the proposed methods.
For the sake of accuracy, we here consider a cardinality constrained OTS
problem, i.e., the formulation 2.1 with the addition of an upper bound on the
number of edges that can be switched off. In formula $\sum_{(i,j)\in
E}(1-x_{ij})\leq L$. For the experiments herein presented, we set the the
upper bound $L$ to 45.

In Table 3, we report a summary of the computational performances of the
proposed approaches using the “best" parameters’ setting displayed in Table 1.
For comparison purposes, we also display the performances of the LWP method.
All the displayed values are average values computed on the set of twenty
instances in the specific class. For each method, Table3 summarizes the
following statistics:

  * •

"rel_gap", i.e., the average relative gap;

  * •

"time", i.e., the average solution time in seconds;

  * •

"$\\#$ u", i.e., the number of instances not solved to optimality within the
time-limit;

  * •

"opt_gap", i.e., the average optimality gap.

Please note that "rel_gap" is not displayed for the LWP method because it
provides the reference point and is zero by definition. Character "_" denote
the case where no observation is available. The best statistic value over the
three method is highlighted with bold fonts for each class of instances. By
and large, the kSP method is outperforming the other approaches. It solved all
the instances at optimality, with short computational times (being the fastest
on five classes of instances out of six), and without compromising the
solution quality. Indeed, it provides better solutions than LWP for two
classes of instances and for two others it provides solutions of the same
quality. The average speedup of kSP over all 120 instances is almost 13 times
with respect to LWP with a peak of 39.8 on the most difficult to solve
instances (e.g., class 118H). The kNN shows very good computational
performances on the classes of instances derived from case118B, especially in
terms of solutions’ quality. However, on instances of the case300B type, we
didn’t get results of the same quality, likely due to a not fully calibrated
tuning phase. Indeed, if we were equipped with an oracle that returns the best
parameters’ setting for each instance we could achieve much better
performance, as demonstrated by the statistics displayed in Table 3. This is
valid for both the methods. For the sake of clarity, for each instance the
best parameters’ setting has been identified as follows: first, we filter out
all the solutions whose objective value is greater than the smallest objective
value obtained for the specific class increased by 0.1%. From the remaining
solutions, we further filter out all the solutions whose computation time is
longer than the minimum computation time plus 1 second. Finally, if there is
more than one solution remaining, we choose the one with the best objective
function value. This table further highlight the potentials of the proposed
approaches.

| kSP | kNN simulation | lwp
---|---|---|---
CoI | rel_gap | time | # u | rel_gap | time | # u | opt_gap | time | # u | opt_gap
118L | 0.00% | 6.5 | 0 | 0.00% | 73.9 | 1 | 0.41% | 126.1 | 4 | 0.35%
118N | -0.01% | 12.9 | 0 | 0.00% | 123.7 | 3 | 0.24% | 138.9 | 4 | 0.23%
118H | 0.03% | 5.2 | 0 | -0.09% | 90.9 | 1 | 0.11% | 207.1 | 6 | 0.48%
300L | -0.01% | 8.5 | 0 | 0.74% | 5.5 | 0 | _ | 16.4 | 0 | _
300N | 0.00% | 64.2 | 0 | 0.53% | 64.9 | 1 | 0.39% | 204.4 | 5 | 0.16%
300H | -0.01% | 50.5 | 0 | 1.5% | 80.8 | 1 | _ | 80.0 | 0 | _

Table 2: Computational results with “best" parameters’ setting per class of
instances.

| kSP | kNN simulation | lwp
---|---|---|---
CoI | rel_gap | time | # u | rel_gap | time | # u | time | # u | opt_gap
118L | -0.01% | 2.6 | 0 | 0.00% | 13.5 | 0 | 126.1 | 4 | 0.35%
118N | -0.03% | 8.9 | 0 | -0.02% | 17.5 | 0 | 138.9 | 4 | 0.23%
118H | 0.02% | 2.9 | 0 | -0.10% | 23.1 | 0 | 207.1 | 6 | 0.48%
300L | -0.01% | 4.6 | 0 | 0.04% | 7.0 | 0 | 16.4 | 0 | _
300N | 0.02% | 13.5 | 0 | 0.12% | 105.2 | 0 | 204.4 | 5 | 0.16%
300H | -0.01% | 17.1 | 0 | 0.16% | 32.6 | 0 | 80.0 | 0 | _

Table 3: Computational results with instance specific parameters’ setting.

###  4.3 Computational analysis of the cardinality constraints.

In practice, power system operators may want to limit the option of switching
off branches of the grid for system’s reliability reasons. We here further
analyze the impact of bounding the number of branches that can be feasibly
switched-off on the computational performances of the proposed methods.
Because these bounds may have considerable influence on the DC-OTS solution
times, it is desirable that the proposed methods scale well with this limit.
For this analysis, we only consider the nominal scenarios, i.e., the classes
of instances 118N and 300N. Our findings confirm the results of [6, 17] about
the increased computational challenge of solving the DC-OTS problem with a
restricted number of edges that can be feasibly switched off. Indeed, smaller
is the parameter $L$ (i.e., the number of edges that can be switched off) and
more difficult are the instances to be solved. The results are summarized in
Table 5 and 5.

kSP outperforms LWP in both quality of solutions - with a negative relative
gap in seven out of eight cases across the two classes of instances considered
- and computational time. The more challenging it is the instance of the
problem, the larger it is the potential benefit of using kSP. kNN finds more
difficulties in taming the complexity of the problem, although for some of the
cases considered it computes good quality solution with negative relative
gaps.

| kSP | kNN simulation | lwp
---|---|---|---
L | rel_gap | time | # u | opt_gap | rel_gap | time | # u | opt_gap | time | # u | opt_gap
10 | 0.04% | 63.7 | 1 | 0.48% | -0.03% | 412.9 | 9 | 2.78% | 407.9 | 11 | 6.50%
15 | -0.06% | 220.6 | 6 | 2.26% | 0.03% | 350.1 | 10 | 5.22% | 345.1 | 9 | 4.06%
20 | -0.11% | 65.2 | 1 | 0.51% | -0.04% | 275.6 | 8 | 2.99% | 212.5 | 6 | 3.18%
25 | -0.04% | 41.3 | 1 | 0.26% | -0.01% | 296.2 | 8 | 2.22% | 177.5 | 4 | 1.99%

Table 4: Computational comparison across different values of L for the 118N
class of instances.

| kSP | kNN simulation | lwp
---|---|---|---
L | rel_gap | time | # u | opt_gap | rel_gap | time | # u | opt_gap | time | # u | opt_gap
10 | -0.77% | 600.0 | 20 | 23.3% | -0.51% | 573.3 | 19 | 27.00% | 600.0 | 20 | 118%
15 | -0.16% | 600.0 | 20 | 6.7% | 0.44% | 600.0 | 20 | 10.43% | 600.0 | 20 | 37%
20 | -0.09% | 180.4 | 2 | 0.3% | 0.52% | 411.1 | 1 | 5.27% | 473.5 | 12 | 29%
25 | -0.01% | 86.1 | 1 | 0.1% | 0.53% | 245.1 | 0 | 0.63% | 384.9 | 7 | 3%

Table 5: Computational comparison across different values of L for the 300N
class of instances.

##  5 Conclusion and Further Work

In this paper, we present two data-driven heuristic methods to compute big M
values with the purpose of tightening the mathematical programming formulation
of the DC-OTS problem. Having small enough big M values, it may improve the
performance of optimization solvers without scarifying the solution quality.
The methods herein proposed exploit information on power demands and
generation costs. The computational experience on a set of 120 instances -with
different features in terms of power demand and network structure -
demonstrates the viability of the proposed approaches. More in particular, the
kSP method gives competitive solutions in short computational times with
respect to the current approach, which is based on computing big M values with
the longest-weighted path. Although, not at the same level of effectiveness,
the kNN method also shows good performances especially on certain classes of
instances. The relatively superior performance of kSP method can be attributed
to the fact that it uses problem-specific knowledge. We also note that kSP is
quit robust with respect to parameters’ setting, thus making this method easy
to implement and resilient to demand forecasting errors. On the other hand,
kNN seems more sensitive to parameters’ setting. It is important to
investigate strategies that provide instance specific settings to achieve the
full potential of the proposed approaches. We also verified that both kSP and
kNN simulation method provide better solutions - with an improvement of over
$20\%$ \- than those computed by the state of the art heuristic method [9].
Moreover it will be of interest to consider a more accurate model including
unit commitment constraints and multiple time periods adapting the proposed
heuristics.

## Acknowledgement

We would like to thank to Prof. Adam Letchford for his comments and
suggestions on a preliminary draft of the paper.

## References

  * [1] "description of case 14".  https://matpower.org/docs/ref/matpower5.0/case14.html.  Accessed: 2021-07-16 [Online].
  * [2] Clayton Barrows, Seth Blumsack, and Russell Bent.  Computationally efficient optimal transmission switching: Solution space reduction.  In 2012 IEEE Power and Energy Society General Meeting, pages 1–8. IEEE, 2012.
  * [3] Silvio Binato, Mário Veiga F Pereira, and Sérgio Granville.  A new benders decomposition approach to solve power transmission network design problems.  IEEE Transactions on Power Systems, 16(2):235–240, 2001.
  * [4] J Carpentier.  Optimal power flows.  International Journal of Electrical Power & Energy Systems, 1(1):3–15, 1979.
  * [5] Gianni Codato and Matteo Fischetti.  Combinatorial benders’ cuts for mixed-integer linear programming.  Operations Research, 54(4):756–766, 2006.
  * [6] Salar Fattahi, Javad Lavaei, and Alper Atamtürk.  A bound strengthening method for optimal transmission switching in power systems.  IEEE Transactions on Power Systems, 34(1):280–291, 2018.
  * [7] T. A. Feo and M.G.C. Resende.  Greedy randomized adaptive search procedures.  Journal of Global Optimization, (6):109–133, 1995.
  * [8] Emily B Fisher, Richard P O’Neill, and Michael C Ferris.  Optimal transmission switching.  IEEE Transactions on Power Systems, 23(3):1346–1355, 2008.
  * [9] J David Fuller, Raynier Ramasra, and Amanda Cha.  Fast heuristics for transmission-line switching.  IEEE Transactions on Power Systems, 27(3):1377–1386, 2012.
  * [10] John Gleeson and Jennifer Ryan.  Identifying minimally infeasible subsystems of inequalities.  ORSA Journal on Computing, 2(1):61–63, 1990.
  * [11] Kory W Hedman, Michael C Ferris, Richard P O’Neill, Emily Bartholomew Fisher, and Shmuel S Oren.  Co-optimization of generation unit commitment and transmission switching with n-1 reliability.  IEEE Transactions on Power Systems, 25(2):1052–1063, 2010.
  * [12] Kory W Hedman, Richard P O’Neill, Emily Bartholomew Fisher, and Shmuel S Oren.  Optimal transmission switching—sensitivity analysis and extensions.  IEEE Transactions on Power Systems, 23(3):1469–1479, 2008.
  * [13] Kory W Hedman, Shmuel S Oren, and Richard P O’Neill.  A review of transmission switching and network topology optimization.  In 2011 IEEE power and energy society general meeting, pages 1–7. IEEE, 2011.
  * [14] Kory W Hedman, Shmuel S Oren, and Richard P O’Neill.  Optimal transmission switching: economic efficiency and market implications.  Journal of Regulatory Economics, 40(2):111, 2011.
  * [15] Rabih A Jabr, Ravindra Singh, and Bikash C Pal.  Minimum loss network reconfiguration using mixed-integer convex programming.  IEEE Transactions on Power systems, 27(2):1106–1115, 2012.
  * [16] Emma S Johnson, Shabbir Ahmed, Santanu S Dey, and Jean-Paul Watson.  A k-nearest neighbor heuristic for real-time dc optimal transmission switching.  arXiv preprint arXiv:2003.10565, 2020.
  * [17] B. Kocuk, H. Jeon, S. S. Dey, J. Linderoth, J. Luedtke, and X. A. Sun.  A cycle-based formulation and valid inequalities for dc power transmission problems with switching.  Operations Research, (64):922–938, 2016.
  * [18] Karsten Lehmann, Alban Grastien, and Pascal Van Hentenryck.  The complexity of dc-switching problems.  arXiv preprint arXiv:1411.4369, 2014.
  * [19] Cong Liu, Jianhui Wang, and James Ostrowski.  Heuristic prescreening switchable branches in optimal transmission switching.  IEEE Transactions on Power Systems, 27(4):2289–2290, 2012.
  * [20] Miles Lubin and Iain Dunning.  Computing in operations research using julia.  INFORMS Journal on Computing, 27(2):238–248, 2015.
  * [21] Pablo A Ruiz, Justin M Foster, Aleksandr Rudkevich, and Michael C Caramanis.  Tractable transmission topology control using sensitivity analysis.  IEEE Transactions on Power Systems, 27(3):1550–1559, 2012.
  * [22] Jiaying Shi and Shmuel S Oren.  Stochastic unit commitment with topology control recourse for power systems with large-scale renewable integration.  IEEE Transactions on Power Systems, 33(3):3315–3324, 2017.
  * [23] Leonardo Taccari.  Integer programming formulations for the elementary shortest path problem.  European Journal of Operational Research, 252(1):122–130, 2016\.
  * [24] Jun Wu and Kwok W Cheung.  On selection of transmission line candidates for optimal transmission switching in large power networks.  In 2013 IEEE Power & Energy Society General Meeting, pages 1–5. IEEE, 2013.
  * [25] Jin Y Yen.  Finding the k shortest loopless paths in a network.  management Science, 17(11):712–716, 1971.

##  Appendix A Computational Results on Relative Gaps and Big-M Values

Some information on relative gaps are shown in Table 9, 9 and 9. "max_gap",
"min_gap", "avg_gap", "std_gap" denote maximum, minimum, average and standard
deviation of relative gaps, respectively. "# non_neg" denotes the number of
instances that have a relative gap no smaller than $0.001\%$. Note that larger
value indicates inferior performance.

Some information on big-M values is shown in Table 9. We select one instance
in each data set and calculate the ratio between big M parameter computed by
our proposed method and big M parameter computed by the longest path method
for each edge in that instance. "max_ratio", "avg_ratio", "min_ratio" and
"std_ratio" denote the maximum, average, minimum and standard deviation of
those ratios, respectively.

| kSP | kNN simulation
---|---|---
CoI | max_gap | min_gap | avg_gap | std_gap | # non_neg | max_gap | min_gap | avg_gap | std_gap | # non_neg
118L | 0.09% | -0.25% | 0.00% | 0.07% | 12 | 0.09% | -0.26% | 0.00% | 0.07% | 9
118N | 0.07% | -0.21% | -0.01% | 0.06% | 10 | 0.17% | -0.13% | 0.00% | 0.07% | 12
118H | 1.50% | -1.52% | 0.03% | 0.58% | 8 | 0.28% | -1.52% | -0.09% | 0.36% | 9
300L | 0.01% | -0.04% | -0.01% | 0.01% | 5 | 4.93% | -0.02% | 0.74% | 1.06% | 17
300N | 0.05% | -0.07% | 0.00% | 0.03% | 10 | 1.50% | -0.02% | 0.53% | 0.47% | 18
300H | 1.67% | -0.07% | 0.16% | 0.46% | 11 | 6.72% | -0.01% | 1.52% | 2.25% | 17

Table 6: Relative Gaps under Different Power Demands

| kSP | kNN simulation
---|---|---
# off | max_gap | min_gap | avg_gap | std_gap | # non_neg | max_gap | min_gap | avg_gap | std_gap | # non_neg
10 | 0.04% | -0.56% | -0.04% | 0.14% | 7 | 0.08% | -0.32% | -0.03% | 0.10% | 1
15 | 0.32% | -0.50% | -0.06% | 0.21% | 8 | 0.77% | -0.36% | 0.03% | 0.24% | 11
20 | 0.04% | -0.81% | -0.11% | 0.20% | 3 | 0.34% | -0.73% | -0.04% | 0.23% | 7
25 | 0.09% | -0.53% | -0.04% | 0.15% | 6 | 0.20% | -0.32% | -0.01% | 0.13% | 12

Table 7: Relative Gaps on Data Set 118N with Different Limits on Number of
Switched-off Edges

| kSP | kNN simulation
---|---|---
# off | max_gap | min_gap | avg_gap | std_gap | # non_neg | max_gap | min_gap | avg_gap | std_gap | # non_neg
10 | 0.00% | -5.31% | -0.77% | 1.50% | 0 | 4.24% | -5.31% | -0.51% | 1.88% | 2
15 | 0.09% | -1.02% | -0.16% | 0.24% | 2 | 1.56% | -0.71% | 0.44% | 0.54% | 15
20 | 0.03% | -0.74% | -0.09% | 0.17% | 5 | 1.43% | -0.63% | 0.52% | 0.55% | 18
25 | 0.08% | -0.10% | -0.01% | 0.04% | 6 | 1.43% | -0.03% | 0.53% | 0.46% | 19

Table 8: Relative Gaps on Data Set 300N with Different Limits on Number of
Switched-off Edges

| kSP | kNN simulation
---|---|---
CoI | max_ratio | avg_ratio | min_ratio | std_ratio | max_ratio | avg_ratio | min_ratio | std_ratio
118L | 100.00% | 24.38% | 1.91% | 33.65% | 100.00% | 31.92% | 0.26% | 31.67%
118N | 100.00% | 24.30% | 1.91% | 33.70% | 100.00% | 32.09% | 0.28% | 31.85%
118H | 100.00% | 24.20% | 1.50% | 33.76% | 100.00% | 37.89% | 0.18% | 41.58%
300L | 100.00% | 42.80% | 0.20% | 45.09% | 100.00% | 27.29% | 0.03% | 38.99%
300N | 100.00% | 42.15% | 1.11% | 45.21% | 100.00% | 75.80% | 0.01% | 41.35%
300H | 100.00% | 41.98% | 1.11% | 45.22% | 100.00% | 81.57% | 0.003% | 37.27%

Table 9: Big-M Values under Different Power Demands


# Quantum Fisher information as a probe for Unruh thermality

Jun Feng  Corresponding author:
<EMAIL_ADDRESS>School of Physics, Xi’an Jiaotong University, Xi’an,
Shaanxi 710049, China  Institute of Theoretical Physics, Xi’an Jiaotong
University, Xi’an, Shaanxi 710049, China  Jing-Jun Zhang  School of Physics,
Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China

###### Abstract

A long-standing debate on Unruh effect is about its obscure thermal nature. In
this Letter, we use quantum Fisher information (QFI) as an effective probe to
explore the thermal nature of Unruh effect from both local and global
perspectives. By resolving the full dynamics of UDW detector, we find that the
QFI is a time-evolving function of detector’s energy gap, Unruh temperature
$T_{U}$ and particularities of background field, e.g., mass and spacetime
dimensionality. We show that the asymptotic QFI whence detector arrives its
equilibrium is solely determined by $T_{U}$, demonstrating the global side of
Unruh thermality alluded by the KMS condition. We also show that the local
side of Unruh effect, i.e., the different ways for the detector to approach
the same thermal equilibrium, is encoded in the corresponding time-evolution
of the QFI. In particular, we find that with massless scalar background the
QFI has unique monotonicity in $n=3$ dimensional spacetime, and becomes non-
monotonous for $n\neq 3$ models where a local peak value exists at early time
and for finite acceleration, indicating an enhanced precision of estimation on
Unruh temperature at a relative low acceleration can be achieved. Once the
field acquiring mass, the related QFI becomes significantly robust against the
Unruh decoherence in the sense that its local peak sustains for a very long
time. While coupling to a more massive background, the persistence can even be
strengthened and the QFI possesses a larger maximal value. Such robustness of
QFI can surely facilitate any practical quantum estimation task.

##  I Introduction

The celebrated Unruh effect predicts that an observer with uniform
acceleration $a$ may perceive the ordinary Minkowski vacuum as a thermal state
of temperature $T_{U}=a/2\pi$ UNRUH1 . It manifests the observer-dependence of
the particle content of quantum fields, which happens naturally in various
scenarios, such as Hawking radiation of black hole UNRUH2  or particle
creation in inflationary universe UNRUH3 . Thus, aside from its own interest,
Unruh effect can also shed light on yet unknown issue of quantum gravity that
are much less tractable (see UNRUH4  for a recent review).

The Unruh effect can be explicitly revealed by the response spectrum of a
moving Unruh-DeWitt (UDW) detector that couples to the field modes of interest
UNRUH5 . In $n=4$ dimensional flat spacetime, once being in equilibrium with a
massless scalar background, the power spectrum of the accelerating detector
exhibits a Planckian form, agreeing precisely with the conventional thermal
spectrum in a heat bath, therefore strongly suggests the interpretation of
$T_{U}$ as a temperature. In an axiomatic way, it has been proven UNRUH6 ;
UNRUH7  that the Unruh effect is universal in arbitrary dimensions
independently of whether or not particles are interacting.

Nevertheless, the thermal nature of Unruh effect is obscure and remains
puzzles. Firstly, an oft-heard explanation on the thermal character of Unruh
radiation relies on so-called thermalization theorem THERMAL1 ; THERMAL2 ,
which indicates that an accelerating detector is causally disconnected from
the degrees of freedom screened by the _global_ Rindler horizon, causing an
information loss presented by a Kubo-Martin-Schwinger (KMS) state, described
by a canonical ensemble. However, this does not naively lead to an observed
detector’s spectrum in Planck form. For example, within a massive scalar
background, the response rate of detector can be far from the form with
Planckian factor. Indeed, the thermalization theorem only guarantees the KMS
condition TAKAGI1 , or equivalently the principle of detailed balance, and
nothing more. In this light, instead of detector’s response spectrum, various
measures of quantum correlation across the horizon, such as entanglement
entropy ee1 ; ee2 , quantum discord Discord1 ; Discord2 , uncertainty relation
F2 ; F3  and nonlocality BF3 ; F5 , have been employed to encapsulate the
_global_ trait of Unruh thermality.

In generic $n$-dimensional spacetime, the detector’s local response may also
exhibit an anomalous ”statistics inversion”, an apparent interchange between
Bose-Einstein and Fermi-Dirac distributions of the spectrum TAKAGI2 ; TAKAGI3
. In particular, for an UDW detector coupled to a massless scalar field,
Takagi shown TAKAGI2  that the excited power spectrum obeys Bose-Einstein
distribution with even $n$, but turns to be Fermi-Dirac for odd $n$, even
without any fermionic degrees of freedom existing. Later pointed in TAKAGI4 ,
this was ascribed to the absence of Huygens principle in odd dimensional
spacetime, a result of the analytic structure of Wightman function, leading to
anticommutativity for timelike separated scalar correlator. Recent works
further argue TAKAGI5 ; TAKAGI6  that the spurious statistics inversion should
be viewed as a _local_ feature of the detector, encoding the nontrivial
density of states on the Rindler wedge, or the way of the detector approaching
its thermal equilibrium TAKAGI7 . Although such particularities of the
response spectrum would not spoil the thermal character of the Unruh effect,
which is solely determined by the KMS condition, it is interesting to ask
TAKAGI8  if any physical meaningfulness can be extracted from those local
response.

In this Letter, we use quantum Fisher information (QFI), a key concept in both
quantum information theory and metrology, to explore the thermal nature of
Unruh effect from both local and global perspectives. In particular, we
estimate the QFI for an UDW detector, linearly coupled to scalar field
background in arbitrary dimensional flat spacetime. We treat the UDW detector
as an open quantum system OP1 , whose dynamics can then be derived from the
time evolution of combined system (detector+external field) by integrating
over the field degrees. In terms of its response function, the density matrix
of detector can be explicitly resolved, which asymptotically reaches a thermal
mixed state OP2 ; OP3 . This approach enables us to investigate Unruh effect
from a refined perspective, rather than the asymptotic excitation rate of the
detector only.

Our interests in using QFI to diagnose Unruh effect are three-fold. Firstly,
as an operational measure on the distinguishability of quantum states, Fisher
information can discriminate parameterized states $\rho({X})$ and
$\rho(X+\delta X)$ with infinitesimal change of parameter $X$. We will argue
that, given the time-dependent state of UDW detector, the Fisher information
$\mathcal{F}_{X}(\rho)$ may encode the particularities of detector’s response
function, e.g., those inherited from the specific coupling ways FI1 ; FI5  or
the statistics inversion in different spacetime dimensions, therefore reveals
the local trait of the Unruh effect. Secondly, as a consequence of the Cramér-
Rao inequality, the QFI quantifies the highest-precision of parameter
measurement achieved in a metrological scheme. For $\theta$ as the parameter
governing relativistic phenomena, intensive works have been made on utilizing
QFI to probe extremely sensitive quantum gravity effects FI2 ; FI3 ; FI3+  or
spacetime structure FI4 ; FI5 ; FI6 . In our scheme, by choosing
$X=\beta\equiv 1/T_{U}$, the QFI bounds the metrological accuracy of any
quantum estimation on Unruh temperature. By efficiently designing an UDW
detector, we may significantly enhance the QFI and preserve it against
environment noise. Finally, since an entangled state provides a higher
precision than an ensemble of uncorrelated particles, the inverse of QFI can
serve as a _genuine_ entanglement probe, based on the metrological performance
of the quantum states FI7 ; FI8 . This suggests that once obtaining the QFI of
UDW detector, we can find the imprint of the non-local nature of Unruh effect
encoding the quantum correlations across the Rindler horizon 111It is
interesting to note FI9 ; FI10  that such connection between Fisher
information and entanglement has recently been identified in holographic
picture, where the holographic dual of QFI near vacuum in CFT gives bulk
entanglement and canonical energy., as well as its surviving in a noisy
environment. Recall that the entangled side of Unruh effect is directly
related to the universal thermalization, and has nothing to do with having a
strictly Planckian response function. This may seem in conflict to the fact
that QFI behaves differently for the same accelerating detector with different
local responses. Nevertheless, we will show that the universal thermal nature
of Unruh effect can be represented by the asymptotically convergent of QFI,
which is ascribed to the KMS condition FI5 , for times much longer than the
scale of thermalized process.

The Letter is organized as follows. In Section II, we solve the master
equation of the UDW detector as an open quantum system in $n$-dimensional flat
spacetime, and give explicitly its time-dependent state. In Section III, we
evaluate the QFI of UDW detector, whose time evolution encodes the local
response in different scenarios and asymptotically convergent represents the
universal thermal nature of Unruh effect. In Section IV, the summary and
discussion are given. Throughout this Letter, we use units with
$G=c=\hbar=k_{B}=1$.

##  II Dynamical evolution of UDW detector in $n$-dimensional spacetime

###  II.1 Master equation and KMS condition

To proceed, we should first explore the full dynamics of an accelerating UDW
detector in $n$-dimensional Minkowski spacetime. Treating the detector as an
open quantum system coupled to a bath of fluctuating quantum scalar field, its
density matrix is governed by a master equation in Lindblad form and evolves
non-unitarily, due to the environment decoherence and dissipation on the
system.

The total Hamiltonian of the combined system of detector and quantum field is

$H=H_{\text{atom}}+H_{\Phi}+\mu H_{I}$ (1)

where UDW detector is modeled by a two-level atom
$H_{\text{atom}}=\frac{1}{2}\omega\sigma_{3}$, $H_{\Phi}$ is the Hamiltonian
of free scalar field $\Phi(x)$, which satisfies the standard Klein-Gordon
equation in $n$-dimensional Minkowski spacetime. The atom-field interaction is
given by $H_{I}=(\sigma_{+}+\sigma_{-})\Phi(x)$. Here, we denote the atomic
raising and lowering operators as $\sigma_{\pm}$, and $\omega$ is the energy
level spacing of the atom.

Assuming the detector and environment are weakly coupled, such that the
initial state of combined system can be approximated as
$\rho_{tot}(0)=\rho(0)\otimes|0\rangle\langle 0|$, where $\rho(0)$ is the
initial state of the atom and $|0\rangle$ is the field vacuum respecting
spacetime symmetry. Obviously, the dynamics of $\rho_{tot}$ should be governed
by a unitary evolution via von Neumann equation
$\dot{\rho}_{tot}(\tau)=-i[H,\rho_{tot}(\tau)]$, where $\tau$ is the proper
time of the detector. Whenever the typical time scale of the environment is
much smaller than that of the detector, we can further assume it is under a
Markovian evolution. The reduced dynamics of the detector is then obtained by
integrating over the background field degrees from the $\rho_{tot}(\tau)$,
governed by a quantum dynamical semigroup of completely positive map222The
Markov approximation enables the upper limit of the integral can be lifted to
infinite proper time, which is permissible as the integrand disappear for the
large time over which the field correlation function decay sufficiently fast.
This leads the Kossakowski coefficients (6) to be integrated over an infinite
range of detector’s proper time OP1 ., and is represented by a Kossakowski-
Lindblad master equation OP4 ; OP5

$\frac{d\rho}{d\tau}=-i\left[H_{\mbox{\tiny
eff}},\rho(\tau)\right]+\mathcal{L}\left[\rho(\tau)\right]$ (2)

where

$\mathcal{L}\left[\rho\right]=\sum^{3}_{i,j=1}C_{ij}\left[\sigma_{j}\rho\sigma_{i}-\frac{1}{2}\left\\{\sigma_{i}\sigma_{j},\rho\right\\}\right]$
(3)

represents a dissipative evolution attributed to the interaction between the
detector and external fields. The Kossakowski matrix $C_{ij}$ can be
explicitly resolved. After introducing the Wightman function of scalar field
$G^{+}\left(x,x^{\prime}\right)=\left\langle
0\left|\Phi(x)\Phi\left({x}^{\prime}\right)\right|0\right\rangle$, its Fourier
transform

$\mathcal{G}(\lambda)=\int_{-\infty}^{\infty}d\tau~{}e^{i\lambda\tau}G^{+}(\tau)=\int_{-\infty}^{\infty}d\tau~{}e^{i\lambda\tau}\left\langle\Phi(x(\tau))\Phi(x(0))\right\rangle$
(4)

determines the coefficients $C_{ij}$ through the decomposition

$C_{ij}=\frac{\gamma_{+}}{2}\delta_{ij}-i\frac{\gamma_{-}}{2}\epsilon_{ijk}n_{k}+\gamma_{0}\delta_{3,i}\delta_{3,j}$
(5)

where

$\gamma_{\pm}=\mathcal{G}(\omega)\pm\mathcal{G}(-\omega),~{}~{}~{}\gamma_{0}=\mathcal{G}(0)-\gamma_{+}/2$
(6)

up to an overall factor $\mu^{2}$ from detector-field coupling, which usually
can be set to 1 for the sake of simplicity. Moreover, the interaction with
external scalar field would also induce a Lamb shift contribution for the
detector effective Hamiltonian $H_{\mbox{\tiny
eff}}=\frac{1}{2}\tilde{\omega}\sigma_{3}$, in terms of a renormalized
frequency $\tilde{\omega}=\omega+i[\mathcal{K}(-\omega)-\mathcal{K}(\omega)]$,
where
$\mathcal{K}(\lambda)=\frac{1}{i\pi}\mbox{P}\int_{-\infty}^{\infty}d\omega\frac{\mathcal{G}(\omega)}{\omega-\lambda}$
is Hilbert transform of Wightman functions.

For later purposes, we can rewrite (5) in a more instructive form. Following a
trajectory of the accelerating detector, one can find that the field Wightman
function fulfills the KMS condition, i.e., $G^{+}(\tau)=G^{+}(\tau+i\beta)$.
Translating it into frequency space, one has

$\mathcal{G}(\lambda)=e^{\beta\omega}\mathcal{G}(-\lambda)$ (7)

Further using the translation invariance $\langle
0|\Phi(x(0))\Phi(x(\tau))|0\rangle=\langle
0|\Phi(x(-\tau))\Phi(x(0))|0\rangle$, after some algebras, we find that (6)
can be represented as

$\displaystyle\gamma_{+}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}d\tau~{}e^{i\lambda\tau}\langle
0|\left\\{\Phi(\tau),\Phi(0)\right\\}|0\rangle=\left(1+e^{-\beta\omega}\right)\mathcal{G}(\omega),$
$\displaystyle\gamma_{-}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}d\tau~{}e^{i\lambda\tau}\langle
0|\left[\Phi(\tau),\Phi(0)\right]|0\rangle=\left(1-e^{-\beta\omega}\right)\mathcal{G}(\omega),$

holding true for generic interacting fields with or without mass.

It also useful to introduce the ratio

$\gamma\equiv\gamma_{-}/\gamma_{+}=\frac{1-e^{-\beta\omega}}{1+e^{-\beta\omega}}=\tanh(\beta\omega/2)$
(9)

which depends only on the Unruh temperature $T_{U}=1/\beta$ due to the
frequency KMS condition (7), but has nothing to do with local correlator of
the background.

For the UDW detector modeled by a two-level atom, its density matrix can be
expressed in a Bloch form

$\rho(\tau)=\frac{1}{2}\Big{(}1+\sum_{i=1}^{3}\rho_{i}(\tau)\sigma_{i}\Big{)}$
(10)

Substituting it into the master equation (3), the Bloch vector
$\bm{\rho}=(\rho_{1},\rho_{2},\rho_{3})^{T}$ satisfies OP6

$\frac{\partial\bm{\rho}}{\partial\tau}+2\bm{K}\cdot\bm{\rho}+\bm{\eta}=0$
(11)

where $\bm{\eta}=(0,0,-2\gamma_{-})^{T}$ and $\bm{K}$ is

$\bm{K}=\left(\begin{array}[]{ccc}\gamma_{+}+\gamma_{0}&{\tilde{\omega}}/2&0\\\
-{\tilde{\omega}}/2&\gamma_{+}+\gamma_{0}&0\\\
0&0&\gamma_{+}\end{array}\right)$ (12)

Assuming a general initial state
$|\psi\rangle=\sin\frac{\theta}{2}|0\rangle+\cos\frac{\theta}{2}|1\rangle$,
the master equation (2) can be analytically resolved. The time-dependent
components of the state are given as

$\displaystyle\rho_{1}(\tau)$ $\displaystyle=$ $\displaystyle
e^{-\frac{1}{4}\gamma_{+}\tau}\sin\theta\cos\tilde{\omega}\tau$
$\displaystyle\rho_{2}(\tau)$ $\displaystyle=$ $\displaystyle
e^{-\frac{1}{4}\gamma_{+}\tau}\sin\theta\sin\tilde{\omega}\tau$
$\displaystyle\rho_{3}(\tau)$ $\displaystyle=$ $\displaystyle
e^{-\frac{1}{2}\gamma_{+}\tau}\left(\cos\theta+\gamma\right)-\gamma$ (13)

It is easy to see that for large time $\tau\gg 1$, the detector’s state
reaches at a thermal equilibrium with Unruh temperature like

$\rho_{\text{equil}}(\beta)=\frac{1}{2}\left(\begin{array}[]{cc}1-\gamma&0\\\
0&1+\gamma\end{array}\right)=\frac{e^{-\beta
H_{\text{atom}}}}{\operatorname{Tr}\left[e^{-\beta H_{\text{atom}}}\right]}$
(14)

An important insight TAKAGI7  is that even for the same end of thermalization,
the detector can follow different paths within its space of states (13) to
reach thermal equilibrium, parameterized by its local features such as the
initial state preparation and the spacetime dimensionality.

With the reduced density matrix (13), the full information of the nonunitary
evolution of the atom can be extracted. A best known example is the transition
rate of the detector between its initial and final states, i.e.,
$\mathcal{P}_{i\rightarrow
f}(\tau)=\operatorname{Tr}\left[\rho_{f}\rho(\tau)\right]$. When initially
$\rho(0)\equiv\rho_{i}$ and $\rho_{f}$ are the ground and excited states of
the detector, i.e., the eigenstates of the Hamiltonian
$H_{\text{atom}}=\frac{1}{2}\omega\sigma_{3}$ in (1), then the probability for
a spontaneous transition of the detector is OP2

$\mathcal{P}_{\uparrow}=\frac{\gamma_{+}-\gamma_{-}}{2\gamma_{+}}\left(1-e^{-2\gamma_{+}\tau}\right).$
(15)

which gives the transition probability per unit time at $\tau=0$ as

$\Gamma_{\uparrow}=\left.\frac{\partial\mathcal{P}_{\uparrow}}{\partial\tau}\right|_{\tau=0}=\mathcal{G}\left(-\omega\right)$
(16)

This spontaneous excitation rate is what usually referred as the Unruh effect.
We note that $\Gamma_{\uparrow}$ does not contain the renormalized frequency
$\tilde{\omega}$ since it indeed manifests the detection rate of the Rindler
quanta from thermalized Minkowski vacuum, while the atom’s Lamb shift is
responsible for the AC Stark shift associated with these quanta OP7 .

###  II.2 State dynamics with various responses

To derive the complete dynamics of detector’s state (13), we need to specify
$\mathcal{G}(\omega)$ following the trajectory of atom with uniformly
acceleration, like

$x^{0}(\tau)=a^{-1}\sinh a\tau,~{}x^{1}(\tau)=a^{-1}\cosh
a\tau,~{}\bm{x}(\tau)=\text{const}$ (17)

in $n$-dimensional Minkowski spacetime. Expending the free scalar field as

$\Phi(x)=\int d^{n-1}k\left\\{a_{k}U_{k}(x)+{a}_{k}^{+}U_{k}^{*}(x)\right\\}$
(18)

the positive field mode can be read as

$U_{k}=\left[2\Omega_{k}(2\pi)^{n-1}\right]^{-1/2}\exp
i\left(k_{1}x^{1}+\bm{k}\cdot\bm{x}-\Omega_{k}x^{0}\right)$ (19)

where $\Omega_{k}^{2}=m^{2}+k_{1}^{2}+|\bm{k}|^{2}$ and the
creation/annihilation operators obey canonical communication relation
$[a_{k},a_{k^{\prime}}^{{\dagger}}]=\delta^{n-1}\left(k-k^{\prime}\right)$.

Substituting (18) and (17) into (4), for an accelerating UDW detector, we can
formally obtain

$\mathcal{G}_{n}(\omega)=-\frac{\pi}{\omega}\frac{D_{n}(\omega)}{e^{-\beta\omega}-1}$
(20)

where the profile function $D_{n}(\omega)$ is given by an integral of modified
Bessel function TAKAGI1

$D_{n}(\omega)=\frac{2}{\pi|\Gamma(i\omega/a)|^{2}}\int\frac{d^{n-2}k}{(2\pi)^{n-2}}\left|K_{i\omega/a}\left(\sqrt{m^{2}+|\bm{k}|^{2}}/a\right)\right|^{2}$
(21)

For a free massless scalar field in $n$-dimension, (20) has an analytical form

$\mathcal{G}^{m=0}_{n}(\omega)=\frac{\pi^{\frac{n-5}{2}}\beta^{3-n}}{2\Gamma\left(\frac{n-1}{2}\right)}\left|\Gamma\left(\frac{n}{2}-1+\frac{\beta\omega}{2\pi}i\right)\right|^{2}\frac{f_{n}(\omega)}{e^{-\beta\omega}-(-1)^{n}}$
(22)

with

$f_{n}(\omega)=\left\\{\begin{array}[]{ll}-\operatorname{sinh}(\beta\omega/2)&\text{
if }n\text{ is even }\\\ \operatorname{cosh}(\beta\omega/2)&\text{ if }n\text{
is odd }\end{array}\right.$ (23)

We observe a Plankian factor with Bose-Einstein distribution for even $n$, but
Fermi-Dirac distribution for odd $n$, that clearly manifests Takagi’s
statistics inversion. This phenomenon also convinces that Unruh effect is
essentially distinct from a simple thermal noise. To see it, one only need to
evaluate the Wightman function in a thermal state with Unruh temperature
$T_{U}$ TAKAGI2 , i.e.
$\left\langle\Phi(x)\Phi\left({x}^{\prime}\right)\right\rangle_{\beta}$, from
which we evaluate (4) as

$\mathcal{G}_{thermal}\sim-\frac{2^{2-n}\pi^{\frac{1-n}{2}}\omega^{n-3}}{\Gamma\left(\frac{n-1}{2}\right)}\frac{1}{e^{-\beta\omega}-1}$
(24)

Apparently, the statistics perceived by the detector remains in arbitrary
spacetime dimensions.

For a free massive scalar field in $n$-dimension, unfortunately, the integral
(21) admits no simple analytic expression. Nevertheless, we can alternatively
evaluate this asymptotically TAKAGI1 , which will be sufficient for our
concern on the possible mass effect encoded in QFI. In particular, for large
field mass with $m\gg T_{U}$, using the asymptotic form of the modified Bessel
function for large argument, the profile function $D_{n}(\omega)$ can be
evaluated, which gives (20) like

$\mathcal{G}^{\text{massive}}_{n}(\omega)\approx\frac{m^{n/2-2}e^{-m\beta/\pi}}{2^{n/2}\beta^{n/2-1}}e^{\beta\omega/2}$
(25)

Although does not contain Planck factor in an explicit way, (25) is still
thermal in the sense that it fulfills the frequency KMS condition, i.e.,
$\mathcal{G}^{\text{massive}}_{n}(-\omega)=e^{-\beta\omega}\mathcal{G}^{\text{massive}}_{n}(\omega)$.

Substituting (22) and (25) into (LABEL:eq8), one can obtain the Kossakowski
coefficients for scalar background with large mass

$\gamma_{+,~{}n}^{\text{massive}}=\frac{m^{n/2-2}e^{-m\beta/\pi}}{2^{n/2-1}\beta^{n/2-1}}\operatorname{cosh}(\beta\omega/2)$
(26)

and for massless scalar field

$\gamma_{+,~{}n}^{m=0}=\frac{\pi^{\frac{n-5}{2}}\beta^{3-n}}{2\Gamma\left(\frac{n-1}{2}\right)}\left|\Gamma\left(\frac{n}{2}-1+\frac{\beta\omega}{2\pi}i\right)\right|^{2}\operatorname{cosh}(\beta\omega/2)$
(27)

Using Euler’s reflection formula $|\Gamma(1/2+ix)|^{2}=\pi/\cosh(\pi x)$ and
recurrence relation $\Gamma(z+1)=z\Gamma(z)$, we give the first several cases
of (27) as

$\displaystyle\gamma_{+,~{}3}^{m=0}=\frac{1}{2}$ ,
$\displaystyle\gamma_{+,~{}4}^{m=0}=\frac{\omega}{2\pi}\gamma^{-1}$
$\displaystyle\gamma_{+,~{}5}^{m=0}=\frac{\pi}{8\beta^{2}}+\frac{\omega^{2}}{8\pi}$
,
$\displaystyle\gamma_{+,~{}6}^{m=0}=\left(\frac{\omega}{3\beta^{2}}+\frac{\omega^{3}}{12\pi^{2}}\right)\gamma^{-1}$
(28)

for later use. The hyperbolic factor $\gamma^{-1}$ appears only for even $n$,
inherited from the statistics inversion.

We note that dimension $n=3$ is quite unique in the sense that its Kossakowski
coefficient is a constant, while those for other dimension are modulated by a
polynomial on $\beta$. Later, we will see that it is such speciality of $n=3$
making the derived QFI evolves drastically different from the models with
other dimension.

##  III QFI for UDW detector and Unruh thermality

Our aim is to explore the QFI of detector’s state (13), with Unruh temperature
$\beta=1/T_{U}$ as a parameter chosen to be optimally estimated. We anticipate
that the different time evolution of QFI can manifest the local feature of
Unruh effect, ascribed to various detector’s response under the interaction
with scalar background. Moreover, within a metrological task on estimating
Unruh temperature, QFI bounds the maximal precision of estimator. Therefore,
how to effectively enhance the QFI by properly design an UDW detector, e.g.,
carefully selecting its initial state, level spacing as well as its
interacting way to the background, may also be revealed for the particular
cases fixed by (26) and (27). Finally, as the equilibrium state (14) depends
solely on Unruh temperature, we would naturally expect that asymptotic
behavior of QFI should represent the thermal nature of Unruh effect which is
guaranteed by the KMS condition alone.

###  III.1 Fisher information for UDW detector

In any quantum metrological tasks, QFI gives a lower bound to the mean-square
error in the parameter estimation via the Cramér-Rao inequality
$\operatorname{Var}(X)\geqslant\left[N\mathcal{F}_{X}\right]^{-1}$, where $N$
is the number of repeated measurements. In terms of the symmetric logarithmic
derivative (SLD) operator $L_{X}$, which satisfies
$\partial_{X}\rho=\frac{1}{2}\\{\rho,L_{X}\\}$, the QFI is defined as
$\mathcal{F}_{X}=\mbox{Tr}[\rho(X)L^{2}_{X}]$. For a diagonalized density
matrix like $\rho=\sum_{i=\pm}\lambda_{i}|\psi_{i}\rangle\langle\psi_{i}|$,
the related QFI can be further written as FI11 ; FI12

$\mathcal{F}_{X}=\sum_{i=\pm}\frac{(\partial_{X}\lambda_{i})^{2}}{\lambda_{i}}+\sum_{i\neq
j=\pm}\frac{2(\lambda_{i}-\lambda_{j})^{2}}{\lambda_{i}+\lambda_{j}}|\langle\psi_{i}|\partial_{X}\psi_{j}\rangle|^{2}$
(29)

where the summations run over all eigenvalues satisfying $\lambda_{i}\neq 0$
and $\lambda_{i}+\lambda_{j}\neq 0$.

In our context, by diagonalizing the detector’s density matrix (10), the
related QFI on parameter $\beta=1/T_{U}$ for arbitrary initial state can be
straightforwardly calculated as FI5

$\mathcal{F}_{\beta}=\frac{\left(\partial_{\beta}\ell\right)^{2}}{1-\ell^{2}}+\frac{\left(\partial_{\beta}(\rho_{3}/\ell)\right)^{2}}{1-(\rho_{3}/\ell)^{2}}+\tau^{2}\left(\rho_{1}^{2}+\rho_{2}^{2}\right)\left(\partial_{\beta}\tilde{\omega}\right)^{2}$
(30)

where $\ell\equiv\sqrt{\sum_{i}\rho_{i}^{2}}$ is the length of Bloch vector.
For given mixed state (13), the length of Bloch vector should satisfy
$\ell<1$, which means each term of (30) is positive and be periodic in
$\theta$. By varying $\theta$ only, it is not hard to find FI5  that the QFI
can achieve its maximum once the detector be prepared in the initial state
with $\theta=(2k+1)\pi$ $(k\in\mathbb{Z})$.

The full dynamics of an UDW detector in $n$-dimensional flat spacetime
indicates that the QFI is a time-evolving function as
$\mathcal{F}_{\beta}\left(\tau,\omega,\theta;n\right)$, which relays on the
initial state preparation fixed by $\theta$, the energy level spacing of atom
$\omega$, as well as the particularities of the background field, such as its
mass and spacetime dimensionality $n$, encoded in detector’s local response.
To access the highest precision of estimation, we have to maximize the value
of QFI over all relevant parameters it depends on.

For simplicity, we can choose the initial state with $\theta=\pi$, i.e., the
detector is initially prepared in ground state $|0\rangle$, which gives
$\ell=|\rho_{3}|$. By straightforwardly calculation on (30), we obtain

$\mathcal{F}_{\beta}\left(\tau,\omega;n\right)=\frac{\left[\left(1-e_{n}\right)\partial_{\beta}\gamma-\frac{\tau}{2}e_{n}(1-\gamma)\partial_{\beta}\gamma_{+,n}\right]^{2}}{1-\left[e_{n}+\gamma(1-e_{n})\right]^{2}}$
(31)

with abbreviation ${e}_{n}\equiv\exp\left(-\gamma_{+,n}\tau/2\right)$. This
expression would be used later to analyze the QFI in various scenarios.

In later discussion, we will work with dimensionless parameters by rescaling
the Unruh temperature and field mass apparently as

$\beta\longmapsto\tilde{\beta}\equiv\beta\omega,\quad
m\longmapsto\tilde{m}\equiv m/\omega$ (32)

Particularly cautious is needed when we rescale the proper time to a
dimensionless parameter. Note that the stationary state (14) is reached
asymptotically via the exponential decay of the $\tau$-dependent term in (13).
This means that $\gamma_{+}^{-1}$ has a characteristic time scale for
arbitrary dimension, and its combination with $\tau$ gives a dimensionless
parameter. Therefore, we define a dimensionless time parameter as333This may
seems strange since only $\mu_{4}$ has the dimension $[\text{Time}]^{-1}$.
However, the rescaling should actually be $\mu_{n}=\mu^{2}\omega^{n-3}$ where
the detector-field coupling $\mu^{2}$ has been absorbed before. For $n\neq 4$,
the coupling constant is dimensionful to maintain the combination
$\gamma_{+}\tau$ dimensionless TAKAGI7 .

$\tau\longmapsto\tilde{\tau}\equiv\mu_{n}\tau,\quad$ (33)

where $\mu_{n}=\omega^{n-3}$. However, for convenience, we continue to term
$\tilde{\beta}$, $\tilde{\tau}$ and $\tilde{m}$, respectively, as $\beta$,
$\tau$ and $m$.

###  III.2 Asymptotic QFI and KMS condition

After evolving for enough long time, the QFI (30) of the reached equilibrium
state (14) has a simple form

$\mathcal{F}_{\beta}\Big{|}_{\tau\rightarrow\infty}=\frac{1-\gamma^{2}}{4}$
(34)

as depicted in Fig.1. One can observe that the asymptotic QFI has nothing to
do with the background field, i.e., be irrelevant to the local response of the
detector, therefore representing the global thermal nature of Unruh effect.

Figure 1: The QFI for the asymptotic equilibrium state of an accelerating UDW
detector, which is independent on the its local response to the background and
represents the global thermal nature of Unruh effect.

As Unruh temperature growing, the QFI value increases which means higher
precision of parameter-estimation can be achieved for the detector with larger
acceleration. This is reasonable in the sense that it is always easy to figure
a distinctive signal, e.g., high $T_{U}$, out from the average background. In
a combined limit of infinite acceleration and long evolving time, we have the
asymptotic QFI as $\mathcal{F}_{\text{asym}}\left(\beta\rightarrow
0\right)=0.25$.

###  III.3 QFI and spacetime dimensionality

We now estimate the QFI of an accelerating detector within massless field
background. For arbitrary dimensional flat spacetime, the detector’s response
function exhibits a statistics inversion manifested in (22). Our aim now is to
exploit that how spacetime dimensionality may influence the QFI, through its
evolution in various dimension.

Without loss of generality, we work in particular cases with $n=3,4,5,6$ and
substitutes the Kossakowski coefficients (28) into (31). We depict the time-
evolving QFI as a function of Unruh temperature in Fig.2.

The particularity of $n=3$ model immediately steps out via an unique behavior
of its QFI. This can be ascribed to the $\beta$-independence of the
Kossakowski coefficient, i.e., $\partial_{\beta}\gamma_{+,3}=0$, which
simplify (31) into a monotonic function of $\beta$ and proper time $\tau$.
With growing time and increasing Unruh temperature, the QFI increases and
monotonically approaches an asymptotic value, e.g.,
$\mathcal{F}_{\text{asym}}=0.25$. In other word, for $n=3$ case, the optimal
precision of estimation can only be achieved for both long enough evolution
time and extreme large acceleration.

For the detector moving in $n=4,5,6$ dimension, however, we find that the QFI
has a peak value located at early time and finite acceleration, and be
significantly larger than their asymptotic value. This suggests that one may
achieve an enhanced precision of estimation on Unruh temperature at a relative
low acceleration. As we mentioned, such non-monotonicity of the QFI is
inherited directly from the polynomial factor on $\beta$ appearing in
Kossakowski coefficient (29), therefore hold for all models with massless
scalar background in $n\neq 3$ dimension.

(a) $\mathcal{F}_{\beta}(n=3)$ (b) $\mathcal{F}_{\beta}(n=4)$

(c) $\mathcal{F}_{\beta}(n=5)$ (d) $\mathcal{F}_{\beta}(n=6)$

Figure 2: The QFI for the accelerating UDW detector coupled with a massless
scalar background in $n=3,4,5,6$-dimensional Minkowski spacetime, as a
function of dimensionless time parameter $\tau$ and temperature parameter
$T_{U}/\omega$. For large acceleration and long time evolution, the QFI
approaches an asymptotic value $\mathcal{F}_{\text{asym}}=1/4$.

It is interesting to note that the location and magnitude of the QFI peak
value can be used to distinguish the different spacetime dimension. We
illustrate this in Fig.3 by comparing the refined evolution of QFI in various
spacetime dimension. Apart from the striking monotonic behavior in $n=3$, we
see that the QFI in higher dimensional spacetime is more sensitive to the
Unruh temperature and arrives at its maximum at a relative earlier time, which
indicates a more sharp peak during its time evolution. Numerically analysis
also shows that the peak value of QFI would increase as the spacetime
dimension growing.

Figure 3: The peaks of the QFI for the accelerating detector in $n=3,4,5,6$
dimensional Minkowski spacetime. The time-evolution of QFI for specific
dimensionless parameter $T_{U}/\omega=0.15,0.5,1$ are given respectively.
Within higher dimensional spacetime, the QFI has a bigger maximal value which
is approached in an earlier time.

Finally, comparing the QFI evolution and detector’s response (22), a crucial
difference between them need to be highlighted, where no inversion behavior of
QFI as a function of dimensionality was found. This may not be surprising
since the QFI evolution essentially distinguishes the way of the detector
thermalizing, which is irrelevant with spurious interchange of statistics
encoded in the functional form of response function TAKAGI3 .

###  III.4 QFI and mass effect

For a massive scalar background, the Kossakowski coefficients (26) is mass-
dependent, indicating that the related QFI of an accelerating detector may
encode somewhat mass effect of the field.

The first mass effect we note is that, with a factor $e^{-m\beta/\pi}$
appearing, the coefficients $\gamma_{+,~{}n}^{\text{massive}}$ are always
modulated by a polynomial on $\beta$. Therefore, we conclude that the QFI
should now be a non-monotonic function of $\beta$ for arbitrary spacetime
dimension, i.e., no such monotonic QFI as shown in Fig.II.(a) within massless
background can be found once the field acquires mass.

To further explore the mass effect encoded in (31), we depict the time-
evolving QFI for $n=3,4,5,6$ massive background with fixed Unruh temperature
in Fig.4. Since analysis should be enrolled in the large field mass limit
$m\gg T_{U}$, we restrict the estimation with low Unruh temperature, e.g.,
$T_{U}/\omega=0.5$, $m/\omega=50$. We observe that the non-monotonic QFI can
reach a maximal value much larger than the asymptotic value
$\mathcal{F}_{\text{asym}}=0.25$, indicating a significant enhancement of the
precision of quantum estimation. Similar as in massless background, as the
spacetime dimension growing, the QFI has a more sharp peak which has a larger
maximum located at an earlier time.

Figure 4: The robust QFI for the accelerating detector coupled with a massive
background in $n=3,4,5,6$ dimensional spacetime. The QFI has peak value
located at early time, which increases for the growing spacetime dimension.
For the scalar background with larger mass, the QFI is significantly enhanced.

Nevertheless, it is interesting to note that the mass-dependent QFI can be
very robust against the environment decoherence, as its peak sustain for a
significantly longer time than it in massless case, whenever the quantum
estimation process can be implemented. This is ascribed to the fact that the
response of the UDW detector is now suppressed by an exponential factor of the
mass, which makes the detector attain to its equilibrium in a very long time
TAKAGI1 . To demonstrate this, we have compared the QFI for the backgrounds
with different mass in Fig.4, by plotting the case of dimensionless parameter
$m/\omega=50,51$ with solid and dot-dashed line respectively. We find that,
for a more massive field, the persistence of the related QFI is strengthened
while its maximal value increase.

##  IV Conclusions

We have reinvestigated the thermal nature of Unruh effect by using quantum
Fisher information as an effective probe. We treat an accelerating UDW
detector as an open quantum system coupled to $n$-dimensional scalar
background, whose complete dynamics is resolved from a Lindblad master
equation of density matrix. The related time-evolving QFI (31) depends on
detector’s energy gap, the Unruh temperature $T_{U}$, as well as the features
of the scalar background, such as its mass and dimensionality. We find that
the asymptotic QFI whence detector arrives equilibrium is irrelevant to the
local response of the detector but only determined by the Unruh temperature
$T_{U}$, thus encapsulates the global side of Unruh thermality guaranteed by
the KMS condition. On the other hand, we show that the local side of Unruh
effect, i.e., the different ways to approach the same thermal equilibrium, is
encoded in the time-evolution of the related QFI. Within massless scalar
background, the QFI has unique monotonicity in $n=3$ dimensional spacetime,
and becomes non-monotonous for $n\neq 3$ models, where a local peak value at
finite acceleration can exist. In a metrological view, this means that an
enhanced precision of estimation on Unruh temperature can be achieve at a
relative low acceleration. Finally, we emphasized that once the field
acquiring mass, the related QFI is robust against the Unruh decoherence, in
the sense that its local peak can sustain for a very long time. For the
background field with larger mass, the persistence of the QFI can even be
significantly strengthened with larger maximal value. Such robustness of QFI
can surely facilitate any practical quantum estimation task.

Our results agree with other previous studies TAKAGI5 ; TAKAGI6 ; TAKAGI7 ,
claiming that the thermal nature of Unruh effect is correct and has nothing to
do with the particularities of detector’s local response. With utilizing of
QFI, we can further refine the physical information extracted from the local
response as the maximal accuracy of possible quantum estimation process FI5 .
It would be interesting to extend our study to other background field, like
the fermionic fields considered in TAKAGI3  or those related to intermediated
statistics TAKAGI9 , where the analysis presented in this Letter can apply.

## Acknowledgement

This work is supported by the National Natural Science Foundation of China
(No.12075178). J. F. would like to thank Yao-Zhong Zhang for enlightening
discussions and kind hospitality at the University of Queensland where part of
this work started.

## References

## References

  * (1) W. G. Unruh, Phys. Rev. D 14, 870 (1976).
  * (2) S. W. Hawking, Commun. Math. Phys. 43, 199 (1975).
  * (3) G. W. Gibbons and S. W. Hawking, Phys. Rev. D 15, 2738 (1977).
  * (4) L. C. B. Crispino, A. Higuchi and G. E. A. Matsas, Rev. Mod. Phys. 80, 787 (2008).
  * (5) B. S. DeWitt, Quantum gravity: the new synthesis, in _General Relativity: An Einstein Centenary Survey_ , edited by S.W. Hawking and W. Israel (Cambridge University Press, Cambridge, 1979), pp. 680.
  * (6) G. L. Sewell, Ann. Phys. 141, 201 (1982).
  * (7) K. Fredenhagen and R. Haag, Commun. Math. Phys. 127, 273 (1990).
  * (8) S. A. Fulling, Phys. Rev. D 7, 2850 (1973).
  * (9) P. C. W. Davies, J. Phys. A 8, 609 (1975).
  * (10) S. Takagi, Prog. Theor. Phys. Suppl. 88, 1 (1986).
  * (11) L. Bombelli, R. Koul, J. Lee and R. Sorkin, Phys. Rev. D 34, 373 (1986).
  * (12) A. Pozas-Kerstjens and E. Martín-Martínez, Phys. Rev. D 92, 064042 (2015).
  * (13) E.G. Brown, K. Cormier, E. Martín-Martínez and R.B. Mann, Phys. Rev. A 86, 032108 (2012).
  * (14) J. Doukas, E. G. Brown, A. Dragan and R. Mann, Phys. Rev. A 87, 012306 (2013).
  * (15) J. Feng, Y. Z. Zhang, M. D. Gould and H. Fan, Phys. Lett. B 726, 527 (2013).
  * (16) J. Feng, Y. Z. Zhang, M. D. Gould and H, Fan, Europhys. Lett. 122, 60001 (2018).
  * (17) N. Friis, P. Köhler, E. Martín-Martínez and R. A. Bertlmann, Phys. Rev. A 84, 062111 (2011).
  * (18) J. Feng, X. Huang, Y. Z. Zhang and H. Fan, Phys. Lett. B 786, 403 (2018).
  * (19) S. Takagi, Prog. Theor. Phys. 74,142 (1985).
  * (20) J. Louko and V. Toussaint, Phys. Rev. D 94, 064027 (2016).
  * (21) H. Ooguri, Phys. Rev. D 33, 3573 (1986).
  * (22) L. Sriramkumar, What do detectors detect, in _Gravity and the Quantum_ , edited by J. S. Bagla and S. Engineer (Springer Publishing 2017)
  * (23) L. Sriramkumar, Gen. Rel. Grav. 35, 1699 (2003).
  * (24) J. Arrechea, C. Barceló, L. J. Garay and G. García-Moreno, Phys. Rev. D 104, 065004 (2021).
  * (25) S. A. Fulling, Particles do not exist, in _Quantum Theory of Gravity_ , edited by S. M. Christensen (CRC Press 1984).
  * (26) H. P. Breuer and F. Petruccione, _The Theory of Open Quantum Systems_ (Oxford University Press 2002).
  * (27) F. Benatti and R. Floreanini, Phys. Rev. A 70, 012112 (2004).
  * (28) D. Moustos and C. Anastopoulos, Phys. Rev. D 95, 025020 (2017).
  * (29) B. Garbrecht and T. Prokopec, Class. Quantum Grav. _21_ , 4993 (2004).
  * (30) M. Aspachs, G. Adesso and I. Fuentes, Phys. Rev. Lett. 105, 151301 (2010).
  * (31) M. Ahmadi, D.E. Bruschi, C. Sabín, G. Adesso and I. Fuentes, Sci. Rep. 4 4996 (2014).
  * (32) Z. Zhao, Q. Pan and J. Jing, Phys. Rev. D 101 056014 (2020).
  * (33) J. Wang, Z. Tian, J. Jing and H. Fan, Nucl. Phys. B 892, 390 (2015).
  * (34) X. Huang, J. Feng, Y. Z. Zhang and H. Fan, Ann. Phys. 397, 336 (2018).
  * (35) H. Du and R. B. Mann, J. High Energy Phys. 05, 112 (2021).
  * (36) L. Pezzé and A. Smerzi, Phys. Rev. Lett. 102, 100401 (2009).
  * (37) P. Hyllus, _et al_., Phys. Rev. A 85, 022321 (2012).
  * (38) N. Lashkari and M. Van Raamsdonk, J. High Energy Phys. 04, 153 (2016).
  * (39) S. Banerjee, J. Erdmenger and D. Sarkar, J. High Energy Phys. 08, 001 (2018).
  * (40) V. Gorini, A. Kossakowski and E.C.G. Surdarshan, J. Math. Phys. 17, 821 (1976).
  * (41) G. Lindblad, Commun. Math. Phys. 48, 119 (1976).
  * (42) R. Alicki and K. Lendi, _Quantum Dynamical Semigroups and Applications_ (Springer 1987).
  * (43) J. Audretsch, R. Müller and M. Holzmann, Class. Quantum Grav. 12, 2927 (1995)
  * (44) S. L. Braunstein and C. M. Caves, Phys. Rev. Lett. 72, 3439 (1994).
  * (45) M. G. A. Paris, Int. J. Quant. Inf. 7, 125 (2009).
  * (46) S. Ohya, Phys. Rev. D 96, 045017 (2017).