\section{Introduction}

Black hole pair creation has been an interesting subject, in the past few decades, partly because it leads to a deeper understanding of the black hole entropy~\cite{garfinkle1991,garfinkle1994,hawking1995,booth1998}. In fact, it has been shown in~\cite{garfinkle1994} that the pair production amplitude contains a factor of $e^{S_{bh}}$, where $S_{bh}$ denotes the black hole entropy. This is consistent with the view that $e^{S_{bh}}$ measures the number of black hole (internal or surface) states.

Pair creation through quantum mechanical tunneling effect can also explain the production of black holes less massive than a few Solar masses. Due to fermionic degeneracy pressure, gravitational collapse cannot be responsible for the production of these small black holes. A recent upper bound on the mass of cold spherical neutron stars is obtained to be $2.3\, M_{\astrosun}$, obtained using the results of GW170817, where it has been shown that neutron stars more massive than this limit would gravitationally collapse into a black hole ~\cite{shibata2019}.

A possible mechanism to produce black hole pairs is through breaking of a cosmic string~\cite{hr} (see also~\cite{achucarro1995}). The string could also ``fray'' to create black hole pair~\cite{eardley1995}. When there is no background magnetic field, it has been shown in~\cite{ashoorioon2014} that the tension of the strut between the two black holes, $\mu_4$, must be smaller than the tension of the string, $\mu_3$, from each black hole to infinity. Given the observational constraint on the tension, $\mu\leq 2\times 10^{-7}$~\cite{seljak2006}, either case of breaking or fraying of the cosmic string to produce black holes of mass larger than the Planck mass is rare~\cite{hr,eardley1995}.

On the other hand, pair creation of black holes in a background magnetic field, with no string, has been widely studied in the literature~\cite{garfinkle1991,dowker1994,hawking1995}. Using the instanton action, the rate of this process is found in the semiclassical approximation. As we show explicitly in this paper, for observationally viable values of the intergalactic magnetic field, the production rate of black holes of mass larger than the Planck mass is negligibly small.

In this paper we investigate the pair production of magnetically charged black holes on a cosmic string in the presence of a background magnetic field. To find the rate of this process, we use the action, $I$, of instanton interpolating between two following states. The first state is a spacetime with magnetic field and a cosmic string of tension $\mu_3$. We refer to this state as background spacetime. In the second state, which we refer to as physical spacetime, we have two black holes, each connected to a string of tension $\mu_3$ that runs to infinity. There may also be a strut between these two black holes with tension $\mu_4$. This configuration is in the presence of a magnetic field.

For a general value of the background magnetic field, we use numerical analysis and we find that the instanton action decreases as the background magnetic field increases. Since the pair production rate is given by $e^{-I}$, we find that for larger values of the background magnetic field the pair production is more probable.

We find an explicit relation for the action $I$ in the weak field limit. We compare this action, with that of the pair creation of black holes in a background magnetic field and no string. The setup we focus on is theoretical, however, we investigate the cosmological implications with real parameters. Using the cosmologically viable values for the tension of the cosmic string and the background magnetic field, we show that the presence of the string would significantly decrease the value of the action and, hence, increase the production rate.

For the sake of comparison with the black hole pair production, we study magnetically charged monopole-antimonopole pair production on a cosmic string in a background magnetic field. Here the configuration of physical spacetime is like that of the black hole case with black holes replaced by magnetic monopoles. To obtain the action of monopole pair production we use the WKB approximation.

We also heuristically study how the presence of a positive cosmological constant affects the black hole pair creation rate in the presence of a background magnetic field and cosmic string. We find that the pair creation rate is higher in a de Sitter (dS) background. During the inflation, where the cosmological constant has its highest value, black hole pair production is the most probable.

The outline of the paper is as follows. In the next section we review the Ernst metric which describes a pair of black holes connected to semi-infinite cosmic strings, in a background magnetic field. We study the production rate of such pairs in Sec.~\ref{sec:rate} by using the instanton action. The weak field limit of this action is studied in Sec.~\ref{sec:wf} where we provide some examples to show how the presence of a cosmic string would decrease the action. In Sec.~\ref{sec:mp} we study the pair production of monopole-antimonopole on a cosmic string in a background magnetic field.  A heuristic derivation of the action in a dS background is presented in Sec.~\ref{sec:ds}. We conclude our paper in Sec.~\ref{sec:con}. We use the units in which $G=c=\hbar=\epsilon_0=1$.

\section{Black hole pair on a cosmic string}\label{sec:ernst}

Before the pair creation, the spacetime consists of a cosmic string with tension $\mu_3$ in a magnetic field. We call this state, the background spacetime and it will be described by Melvin solution~\cite{melvin1964}
\begin{align}
	ds^2&=\bar{\Gamma}^2\left[-dt^2+dz^2+d\rho^2\right]+\bar{\Gamma}^{-2}\rho^2d\bar{\phi}^2, \label{melvin}\\
	\bar{\Gamma}&=1+\frac{1}{4}\hat{B}_{M}^2\rho^2, \nonumber
\end{align}
where $\hat{B}_{M}$ denotes the magnetic field, and the gauge potential is given by $2\bar{\Gamma}A_{\bar{\phi}}=\hat{B}_{M}\rho^2$. The cosmic string results in a deficit angle $\delta_3=8\pi\mu_3$ in the azimuthal coordinate, so $\Delta\bar{\phi}=2\pi-\delta_3$.

After the pair creation, the spacetime, which now we refer to as physical spacetime, will be given by the Ernst metric~\cite{ernst1976}
\begin{align}
	ds^2&=r^2\Gamma^2\left[G(y)dt^2-G^{-1}(y)dy^2+G^{-1}(x)dx^2\right] \nonumber\\ &+\frac{r^2G(x)}{\Gamma^2}d\phi^2, \label{ernst}\\
	A_{\phi}&=-\frac{2}{B\Gamma}\left(1+\frac{1}{2}Bqx\right)+k, \nonumber
\end{align}
where $r=A^{-1}(x-y)^{-1}$ represents the radial coordinate~\cite{kinnersley1970}, and
\begin{align}
	\Gamma&=\left(1+\frac{1}{2}B q x\right)^2+\frac{1}{4}r^2B^2G(x),\\
	G(\xi)&=(1+r_{-}A\xi)(1-\xi^2-r_{+}A\xi^3). \nonumber
\end{align}
The parameter $q^2=r_{-}r_{+}$ is related to the charge of the black hole and the mass of the black hole is given by $m=(r_{-}+r_{+})/2$. This solution describes two charged black holes accelerating from each other. The acceleration is given by the parameter $A$ with its inverse representing the typical distance between the two black holes at the time of their creation~\cite{dias2004}.

The physically interesting case is the one in which the function $G(\xi)$ has four real roots. One root is given by $\xi_1=-1/(r_{-}A)$. To have three other real roots we need to restrict the parameters so that $r_{+}A\leq\sqrt{4/27}$, which means that the black holes should be small and/or they should be far apart from each other. For $r_{+}A\ll 1$ one finds the remaining roots as follows
\begin{align}
	\xi_2&=-\frac{1}{r_{+}A}+r_{+}A+\cdots,\label{roots}\\
	\xi_3&=-1-\frac{r_{+}A}{2}+\cdots, \nonumber\\
	\xi_4&=1-\frac{r_{+}A}{2}+\cdots. \nonumber
\end{align}

To obtain the correct signature of the metric (\ref{ernst}), we require $\xi_3\leq x \leq \xi_4$ and $-\infty<y\leq x$. The surfaces $y = \xi_1$, $\xi_2$ , $\xi_3$ correspond to black hole inner horizon, event horizon, and acceleration horizon, respectively. Also $x=\xi_3$ and $x=\xi_4$ are axes that point to spatial infinity and to the other black hole, respectively. The parameter $k$, in the gauge potential of Eq. (\ref{ernst}), will be chosen so as to confine the Dirac string singularities to the $x=\xi_4$ axis. It should also be mentioned that the parameter $B$ in the Ernst metric is related to the physical magnetic field $\hat{B}$ by~\cite{dowker1994}
\begin{equation}
\hat{B}=\frac{B\,G'(\xi_3)}{2L^{3/2}}, \label{physb}
\end{equation}
where $L=\Gamma(\xi_3)$. It could be easily shown that in the weak field limit $\hat{B}=B$ to the leading order of background magnetic field. Also, the physical magnetic charge of the black hole is defined by $\hat{q} = 1/(4\pi)\int F$, where $F$ denotes the Maxwell field~\cite{dowker1994}. Integration over a two sphere surrounding the black hole yields
\begin{equation}
\hat{q} = q \frac{\Delta\phi\left(\xi_4-\xi_3\right)}{4\pi L^{1/2}\left(1+\frac{1}{2}B q \xi_4\right)}, \label{eqn:charge}
\end{equation}
where $\Delta\phi$ is the period of the azimuthal coordinate.

Consider the angular part of the Ernst metric (\ref{ernst})
\begin{equation}
d\Phi^2=r^2\left[\frac{\Gamma^2}{G(x)}dx^2+\frac{G(x)}{\Gamma^2}d\phi^2\right].
\end{equation}
Using the change of variable $\theta=\int_{\xi_i}^{x}\Gamma G^{-1/2}(x)dx$, near the poles $x=\xi_i$ ($i=3,4$), we obtain
\begin{equation}
d\Phi^2=r^2\left[d\theta^2+\frac{G'^2(\xi_i)}{4\Gamma^4(\xi_i)}\theta^2d\phi^2\right].
\end{equation}
We are interested in a configuration in which each black hole is connected to a string with deficit angle $\delta_3$ that runs to infinity and a strut with deficit angle $\delta_4$ that goes to the other black hole. In such case, one can show that the following relation holds near the poles
\begin{equation}
\frac{|G'(\xi_i)|}{2\Gamma^2(\xi_i)}\Delta\phi=2\pi-\delta_i, \qquad i=3,4. \label{phiatpoles}
\end{equation}
Now we write the period of the azimuthal coordinate in terms of $\xi_3$
\begin{equation}
\Delta\phi=\frac{2L^2}{G'(\xi_3)}(2\pi-\delta_3)\,. \label{deltaphi}
\end{equation}
It is easy to show that $G'(\xi_3)>0$. By using Eqs. (\ref{phiatpoles}) and (\ref{deltaphi}), we find
\begin{equation}
\left|\frac{G'(\xi_4)}{G'(\xi_3)}\right|\frac{L^2}{\Gamma^2(\xi_4)}=\frac{2\pi-\delta_4}{2\pi-\delta_3}. \label{wbnewton}
\end{equation}
This relation would lead to Newton's law in the weak field limit $Bq\ll1$. In observationally important case we have $\delta\ll1$. We are also working in the limit which the size of the black hole is small compared to the typical distance between them (the so called small black hole limit). Then Eqs. (\ref{roots}) and (\ref{wbnewton}) would result in
\begin{equation}
mA=\left(\mu_3-\mu_4\right)+Bq. \label{newton}
\end{equation}

Now consider the first two terms of the Ernst metric (\ref{ernst}). By Wick rotating the time coordinate $t=i\tau$ and following a similar calculations as above, we find that the period of the imaginary time at the surfaces $y=\xi_i$ ($i=2,3$) is $\Delta\tau=4\pi/|G'(\xi_i)|$. For the case of extreme black hole, $\xi_1=\xi_2$, the surface $y=\xi_2$ is infinitely far away from any point in the Euclidean sector of the metric~\cite{hawking1995}. As a result, one does not worry about the periodicity of $\tau$ at this surface (In this case we take $\Delta\tau=4\pi/G'(\xi_3)$). For the non-extreme case we request that $y=\xi_2$ and $y=\xi_3$ have the same periodicity of $\tau$, or equivalently, have the same temperature $1/\beta$. So
\begin{equation}
\Delta\tau=\beta=\frac{4\pi}{G'(\xi_3)}=-\frac{4\pi}{G'(\xi_2)},
\end{equation}
which leads to the following relation for the non-extremal case
\begin{equation}
\xi_4-\xi_3=\xi_2-\xi_1. \label{nerelation}
\end{equation}

\section{Production rate of black hole pair on a cosmic string in the presence of a background magnetic field}\label{sec:rate}

In this section we are going to obtain the rate of black hole pair production on cosmic string in the presence of a background magnetic field. To do so, we use the action of the instanton that interpolates between background and physical spacetimes. Based on the reasoning of~\cite{hr,hawking1996}, this action is equivalent to the Euclidean action
\footnote{We note that in Eqs. (\ref{neaction}) and (\ref{eaction}) all the surface terms have been taken into account~\cite{hawking1996}. The discussion of~\cite{hawking1996} applies to all spacetimes that can be foliated by complete, nonintersecting spacelike surfaces. In derivation of this formulae one needs to keep in mind that, since the Lorentzian action diverges for spatially non-compact geometries, one needs a reference background. One then fixes a boundary near infinity and requires that the metric and matter fields of the physical spacetime induce the same fields on the boundary as those of the background spacetime induce. Also, it is assumed that the periodicity of time coordinate $\beta$ be the same in physical and background spacetimes.}
\begin{equation}
I=\beta H-\frac{1}{4}\Delta\mathcal{A}, \label{eaction}
\end{equation}
in the extremal black hole case, and
\begin{equation}
I=\beta H-\frac{1}{4}\left(\Delta\mathcal{A}+A_{bh}\right), \label{neaction}
\end{equation}
in the non-extremal case.
\footnote{The extra term in Eq. (\ref{neaction}) relative to Eq. (\ref{eaction}) shows that the production rate in non-extreme case is enhanced by $e^{S_{bh}}$, where $S_{bh}=A_{bh}/4$ is the entropy of the non-extremal black hole, as the entropy of an extremal black hole is zero~\cite{hawking1995}. See also~\cite{emparan1995:2} to find out why only the entropy of one of the black holes enter the equations.}
Here $\Delta\mathcal{A}$ is the difference between the area of the acceleration horizon in physical and background spacetimes. $A_{bh}$ is the area of the black hole, and the Hamiltonian $H$ is given by~\cite{hawking1996}
\begin{equation}
H=\int_{\Sigma}N\mathcal{H}-\frac{1}{8\pi}\int_{S_{\infty}^2}N\left(^2K-{}^2K_0\right), \label{hamilton}
\end{equation}
where $N$ is the lapse function and $\Sigma$ is the constant time hypersurface. $\mathcal{H}$ is the Hamiltonian constraint which vanishes on the solution. $S_{\infty}^2$ is the boundary of $\Sigma$ and $^2K$ and $^2K_0$ are the extrinsic curvature of this boundary in physical and background spacetimes, respectively. Note that $S_{\infty}^2$ has been chosen so as to match in physical and background spacetimes. This way, a calculation similar to that of~\cite{hawking1995} would result in the conclusion that the second integral in Eq. (\ref{hamilton}) also vanishes.

Now, we are going to find the area of the black hole horizon, $A_{bh}$, which appears in the action (\ref{neaction}) for non-extremal black holes. The event horizon is at $y=\xi_2$ in the Ernst metric (\ref{ernst}). We find, by using (\ref{deltaphi}),
\begin{align}
	A_{bh}&=\int_{y=\xi_2}\sqrt{g_{xx}g_{\phi\phi}}dxd\phi \nonumber\\ &=\frac{2L^2}{A^2G'(\xi_3)}\left(2\pi-\delta_3\right)\frac{\xi_4-\xi_3}{\left(\xi_3-\xi_2\right)\left(\xi_4-\xi_2\right)}. \label{bharea}
\end{align}

Now, we proceed to find the area of the acceleration horizons. Recalling the definition of the radial coordinate, $r=A^{-1}(x-y)^{-1}$, and the restrictions on $x$ and $y$, we note that equalities $x = y = \xi_3$ show the spatial infinity. To find the area of the acceleration horizon in the Ernst metric, which is located at $y=\xi_3$, we integrate from $x=\xi_4$ to $x=\xi_3+\epsilon$, and will take the limit $\epsilon\rightarrow0$ at the end. We would find
\begin{align}
	\mathcal{A}&=\int_{y=\xi_3}\sqrt{g_{xx}g_{\phi\phi}}dxd\phi\nonumber\\ &=\frac{2L^2}{A^2G'(\xi_3)}\left(2\pi-\delta_3\right)\left(\frac{1}{\epsilon}-\frac{1}{\xi_4-\xi_3}\right). \label{ernstacch}
\end{align}
The area of the acceleration horizon in the Melvin metric (\ref{melvin}) is
\begin{equation}
\bar{\mathcal{A}}=\int\rho d\rho d\bar{\phi}=\frac{2\pi-\delta_3}{2}\bar{\rho}^2. \label{melvinacch}
\end{equation}
By requiring that the proper length of the boundary and the integral of the gauge potential around the boundary be the same in Ernst and Melvin metrics, it has been found in~\cite{hawking1995} that
\begin{equation}
\bar{\rho}^2=\frac{4L^2}{A^2G'(\xi_3)\epsilon}\left(1+\frac{G''(\xi_3)}{4G'(\xi_3)}\epsilon\right)^2.
\end{equation}
Substituting $\bar{\rho}$ from the above equation into Eq. (\ref{melvinacch}), we obtain
\begin{equation}
\bar{\mathcal{A}}\simeq\frac{2L^2}{A^2G'(\xi_3)}\left(2\pi-\delta_3\right)\left(\frac{1}{\epsilon}+\frac{G''(\xi_3)}{2G'(\xi_3)}\right). \label{melvinacch2}
\end{equation}
To find the second term in the above equation we note that one could write the function $G(\xi)$ as
\begin{equation}
G(\xi)=-(r_-A)(r_+A)(\xi-\xi_1)(\xi-\xi_2)(\xi-\xi_3)(\xi-\xi_4).
\end{equation}
By using Eqs. (\ref{ernstacch}) and (\ref{melvinacch2}) we find that
\begin{align}
	\Delta\mathcal{A}=\mathcal{A}-\bar{\mathcal{A}}&=-\frac{2L^2}{A^2G'(\xi_3)}\left(2\pi-\delta_3\right) \nonumber\\
	&\times\left(\frac{1}{\xi_3-\xi_1}-\frac{1}{\xi_3-\xi_2}\right). \label{acch}
\end{align}
By substituting this equation along with Eq. (\ref{bharea}) into Eq. (\ref{neaction}), and using the non-extremality relation (\ref{nerelation}), we obtain
\begin{equation}
I=\frac{L^2}{A^2G'(\xi_3)}\left(\frac{2\pi-\delta_3}{\xi_3-\xi_1}\right). \label{action}
\end{equation}
The same equation would be found for the action of extremal black hole, by using Eqs. (\ref{eaction}) and (\ref{acch}) and the relation $\xi_2=\xi_1$.

Here we would like to study the action (\ref{action}) for the case in which $q = m$ (or equivalently $r_+ = r_-$). We have restricted the black holes to be small on the scale set by the acceleration. That is, we are considering the case in which $r_+A\ll1$. By using this condition, we could show from (\ref{physb}) that
\begin{equation}
\hat{B}\simeq\frac{8B}{(2-Bm)^{3}}. \label{physb:genb}
\end{equation}
We see that $\hat{B}$ goes from zero to infinity as $B$ changes from zero to $2/m$.

On the other hand, for the case $q = m$, Eq. (\ref{action}) could be written as
\begin{equation}
I\simeq\frac{m(B m-2)^4}{32A}\left(2\pi-\delta_3\right), \label{action:genb}
\end{equation}
where the parameter $A$ has to satisfy Eq. (\ref{wbnewton}). We have not been able to solve Eq. (\ref{wbnewton}) analytically, to find $A$ for a general value of $B$. However by using numerical methods, we have used Eqs. (\ref{physb:genb}), (\ref{action:genb}), and (\ref{wbnewton}) to plot the action as function of physical magnetic field $\hat{B}$ for the $q = m$ case in Fig.~\ref{fig:action}. It is obvious from this figure that the instanton action decreases as the physical magnetic field increases. Therefore, $e^{-I}\rightarrow1$ as $\hat{B}\rightarrow\infty$, which means that the black hole pair production rate increases by increasing the magnetic field.

\begin{figure}[htp]
	\centering
	\includegraphics[width=0.45\textwidth]{action.pdf}
	\caption{The action as a function of the physical magnetic field for the $q = m$ case. We see that the action decreases by increasing the physical magnetic field $\hat{B}$. We have set $m\,(=q)=10$, $\mu_3 = 2\times 10^{-7}$, and $\mu_4 = 10^{-7}$.}
	\label{fig:action}
\end{figure}

\section{Weak field limit of the action}\label{sec:wf}

We see from Fig.~\ref{fig:action} that the action for the pair production of black holes of masses larger than the Planck mass is of order unity, only when the physical magnetic field is about 10 or larger. This is much larger than the value of intragalactic magnetic field, which, nowadays, is less than about $10^{-63}$ in natural units~\cite{ashoorioon2005}. It motivated us to study the weak field limit of the action (\ref{action}).

In the weak field limit, the physical magnetic field (\ref{physb}) and charge (\ref{eqn:charge}) can be approximated by
\begin{align}
	\hat{B}&\simeq B\left[1-\frac{1}{2}Bq-2\left(\mu_3-\mu_4\right)\right], \label{eqn:physb:emp}\\
	\hat{q}&\simeq q\left(1-4\mu_3\right)\left[1+2\left(\mu_3-\mu_4\right)\right],
\end{align}
where we have used Eq. (\ref{deltaphi}). These expressions reduce to the ones derived in ~\cite{emparan1995} when $\mu_3=0$. Using these equations, along with the Newton's law (\ref{newton}) and the non-extremality relation (\ref{nerelation}), one can write all the parameters in terms of $\hat{B}$, $\hat{q}$, $\mu_3$, and $\mu_4$
\begin{align}
	r_{\pm}&\simeq \hat{Q}\left[1-2\left(\mu_3-\mu_4\right)\pm\left(\hat{B}\hat{Q}+\mu_3-\mu_4\right)\right], \label{eqn:para:r}\\
	q&\simeq\hat{Q}\left[1-2\left(\mu_3-\mu_4\right)\right], \label{eqn:para:q}\\
	A&\simeq\hat{B}\left[1+\frac{1}{2}\hat{B}\hat{Q}+\frac{\mu_3-\mu_4}{\hat{B}\hat{Q}}\left(1+2\hat{B}\hat{Q}\right)\right], \label{eqn:para:A}
\end{align}
where $\hat{Q}=\hat{q}\left(1+4\mu_3\right)$. Again see similar expressions in ~\cite{emparan1995}, when $\mu_3=0$.

By substituting the parameters (\ref{eqn:para:r})-(\ref{eqn:para:A}) into Eq. (\ref{action}) we find the weak field approximation of the action as
\begin{equation}
I\simeq\pi\hat{Q}^2\left(1-4\mu_3\right)\left[\frac{1}{\hat{B}\hat{Q}}-\frac{1}{2}-\frac{\mu_3-\mu_4}{\hat{B}^2\hat{Q}^2}\right]. \label{eqn:act:emp}
\end{equation}
If we have no string, $\mu_3 = \mu_4 = 0$, one obtains the result of~\cite{dowker1994}. Hence, Eq. (\ref{eqn:act:emp}) is also in agreement with Schwinger result~\cite{schwinger1951}.
\footnote{Eq. \eqref{eqn:act:emp} agrees to the first three terms of the action obtained in~\cite{emparan1995}. However, if we continue our expansion \eqref{eqn:act:emp}, the fourth term would be $-\pi\hat{Q}\left(1-4\mu_3\right)\left(\mu_3-\mu_4\right)/\hat{B}$ which differs from the fourth term of~\cite{emparan1995} by a factor of 2. We only consider the first three terms of the action \eqref{eqn:act:emp} when we are comparing black hole pair production action to that of monopoles in the next section.}

Now we study the action in the case where $q=m$ and we are also in the weak field limit. In this case, the action (\ref{action}) can be written as
\begin{equation}
I\simeq \frac{\pi m (1-2Bm)}{A}\left(1-4\mu_3\right). \label{action:smallb}
\end{equation}
By using Eqs. (\ref{newton}) and (\ref{eqn:physb:emp}) one can rewrite the above equation as
\begin{equation}
I\simeq \frac{\pi m^2 \left[1-2\left(J+\mu_3-\mu_4\right)\right]\left(1-4\mu_3\right)}{J-1+\left[1+2\left(\mu_3-\mu_4\right)\right]\left(\mu_3-\mu_4\right)}, \label{ours}
\end{equation}
where the auxiliary function $J$ is
\begin{equation}
J=\sqrt{\left[1-2\left(\mu_3-\mu_4\right)\right]\left[1-2\left(\hat{B}m+\mu_3-\mu_4\right)\right]}. \nonumber
\end{equation}

Let us see how the presence of the cosmic string would change the rate of the production. Suppose that a pair of black holes, each of mass $m=10$, pop up in the background magnetic field $\hat{B}=10^{-63}$ which is the today's value of the intragalactic magnetic field. Substituting these values into Eq. (\ref{ours}) with $\mu_3 = \mu_4 = 0$, we find $I\simeq 3.141592\times 10^{64}$ for the action of this pair creation. Let us compare this with the case where $\hat{B}=10^{-63}$, $m=10$, and $\mu_3 = \mu_4 =2\times 10^{-7}$ i.e.~when the pair creation of black holes occurs on a cosmic string with uniform tension. In this case we find $I\simeq 3.141590 \times 10^{64}$, which is slightly smaller than when black holes pair is created from vacuum. This is due to the factor $1-4\mu_3$ in Eq. \eqref{ours}. This implies that the presence of the string, even if it does not fray or break, reduces the action, although for the cosmic strings with observationally relevant values of tension, the change in the action will be quite small.

Now let us assume that we have a string of tension $\mu_3=2\times 10^{-7}$ in the Melvin universe (recall that this is the upper limit on the tension of the cosmic string). Assume that the string frays and produces a pair of black holes with the same mass $m=10$. Let us take $\mu_4=10^{-7}$, for the tension of the strut between the two black holes. Also take the same value of the background magnetic field as above, $\hat{B}=10^{-63}$. Substituting these values into the action (\ref{ours}), we find $I\simeq 3.141590 \times 10^{9}$.
\footnote{We note that if we omit the factor $\left(1-4\mu_3\right)$ in Eq. (\ref{ours}), we would find $I\simeq 3.141592 \times 10^{9}$. Therefore, the contribution of this factor is very small compared to the changes of the action upon breaking/fraying of the cosmic string.}
This is smaller than the value of the action in the absence of the string, by 55 order of magnitude!

In the early universe, the background magnetic field can be much stronger. It is supposed to be of the order $\hat{B}=10^{-13}$ in natural units at the end of inflationary era \footnote{We have assumed that the observed intragalactic field, has been produced during inflation but its strength has reduced super-adiabatically, $B\propto a^{-2}$, since then.}. Using this value of the background magnetic field along with $\mu_3 = \mu_4 = 0$ and $m=10$, we find $I\simeq 3.141592 \times 10^{14}$,  which is smaller than the corresponding value in today magnetic ﬁeld strength by 50 orders of magnitude.

For the case of pair creation on the cosmic string in the early universe, we use $\mu_3=2\times 10^{-7}$, $\mu_4=10^{-7}$, $m=10$, and $\hat{B}=10^{-13}$. By substituting these values into Eq. (\ref{ours}), we find $I\simeq 3.141559 \times 10^{9}$. One notices that the value of the action for pair creation on a cosmic string at the end of the inflationary phase does not differ significantly from its current value, although without cosmic string, the action reduces significantly.

We find from Eq. (\ref{ours}) that for fixed values of the parameters, the action is minimized if the string breaks, i.e.~$\mu_4=0$. For $\hat{B}=10^{-63}$, $m=10$, $\mu_3=2\times 10^{-7}$, and $\mu_4=0$ we find $I\simeq 1.570795 \times 10^{9}$. On the other hand if the cosmic string does not fray/break, one finds the maximum value of the action for a given set of the parameters, even though we should stress that the obtained value in this case is slightly smaller than the value of the action when the pair creation of black holes takes place only in presence of background magnetic field without a cosmic string.

We also see from Fig.~\ref{fig:actionwf} that for fixed $\mu_3$, $\mu_4$, and $\hat{B}$, the action (\ref{ours}) would increase by increasing the mass of the black holes. One concludes that the probability of pair producing more massive black holes is smaller.

\begin{figure}[htp]
	\centering
	\includegraphics[width=0.45\textwidth]{actionwf.pdf}
	\caption{The action for the $q=m$ case in the weak field limit as a function of black hole mass. We see that the action increases by increasing the mass of the black holes. We have taken $\mu_3=2\times 10^{-7}$, $\mu_4=10^{-7}$, and $\hat{B}=10^{-63}$.}
	\label{fig:actionwf}
\end{figure}

\section{Monopole pair production on cosmic strings in a background magnetic field}\label{sec:mp}

Monopoles are hypothetical particles with quantized magnetic charge~\cite{dirac1931}. It has been shown that they can be produced in theories with certain non-Abelian gauge fields~\cite{t1974}. In this section we study the production rate of magnetically charged monopole-antimonopole pair on cosmic string in the presence of a background magnetic field.

In the absence of a background magnetic field, the semiclassical probability for a cosmic string to break producing a monopole-antimonopole pair has been obtained and their cosmological evolution has been studied~\cite{vilenkin1982}. The production rate of monopole pairs in a background magnetic field has also been studied in~\cite{affleck1982}. Recently, the rate of monopole pair production has been investigated in thermal baths~\cite{gould2017} and strong background magnetic fields~\cite{rajantie2019,ho2020}. Investigation of monopole pair production on cosmic string in a background magnetic field is missing in the literature.

Here we use WKB approximation, similar to that of~\cite{vilenkin1982}, to find the action of monopole pair production on cosmic string in an external magnetic field. We compare this result with the action of black hole pair production on a cosmic string in a background magnetic field.

Assume a cosmic string with tension $\mu_3$ on an axis $z$. Assume also a background magnetic field $\hat{B}$ parallel to this axis. Suppose the string breaks and produces a monopole and an antimonopole each connected to a semi-infinite string. Because of the symmetry, we only consider the monopole in the following analyze and include the antimonopole rule by multiplying the monopole action by 2, in the end.

The monopole has mass $m$ and magnetic charge $\hat{q}$. The energy of the monopole as a function of $z$ is $m+\left(\mu_3+\hat{B}\hat{q}\right)z$. We note that if the string frays and we have strut of tension $\mu_4$ between the pair, then the force on monopole due to the string-strut system is $\mu_3-\mu_4$ and the energy would be
\begin{equation}
E=m+\left(\mu_3-\mu_4+\hat{B}\hat{q}\right)z.
\end{equation}
The additive constant has been chosen to make the monopole bounce at $z=0$, see \cite{vilenkin1982}. The action for the pair production is then
\begin{equation}
I_{mon}\simeq 4 \int_{\frac{-m}{\mu_3-\mu_4+\hat{B}\hat{q}}}^{0}\sqrt{m^2-E^2}dz.
\end{equation}
The lower bound of the integral has been chosen so as to make $E\geq 0$ throughout the integration domain. Integrating we find
\begin{equation}
I_{mon}\simeq \pi m^2 \left(\frac{1}{\hat{B}\hat{q}}-\frac{\mu_3-\mu_4}{\hat{B}^2\hat{q}^2}\right). \label{eqn:mnac1}
\end{equation}

It has been suggested in~\cite{dowker1994} that, regarding the production rate, monopoles behave more like extremal black holes than non-extremal ones. This suggestion is based on the calculations of~\cite{garfinkle1994} which shows that the action of the instanton describing monopole pair production is greater than the action of non-extreme black holes by the black hole entropy. On the other hand, it was found in~\cite{dowker1994} that the instanton action for the pair production of extremal black holes is larger than the action for the case of non-extreme black holes again by the black hole entropy. For the extreme black holes we have $m=\hat{q}$. Thus if we regard monopoles as extreme black holes we can rewrite Eq. \eqref{eqn:mnac1} as
\begin{equation}
I_{mon}\simeq \pi \hat{q}^2 \left(\frac{1}{\hat{B}\hat{q}}-\frac{\mu_3-\mu_4}{\hat{B}^2\hat{q}^2}\right). \label{eqn:mnac2}
\end{equation}

Eq. \eqref{eqn:act:emp} is the action for pair production of (slightly) non-extremal black holes, and, in $\mu_3\ll 1$ limit, can be written as the following
\begin{equation}
I\simeq \pi \hat{q}^2 \left(\frac{1}{\hat{B}\hat{q}}-\frac{\mu_3-\mu_4}{\hat{B}^2\hat{q}^2}-\frac{1}{2}\right).
\end{equation}
Comparing this equation with Eq. \eqref{eqn:mnac2}, we find that the action obtained by WKB approximation is greater than the action of non-extremal black holes by half the Bekenstein-Hawking entropy $\pi\hat{q}^2/2=S_{bh}/2$. Similar discrepancies between the actions for non-extremal black hole pair creation, \cite{dowker1994}, the exact action for monopole pair creation, \cite{affleck1982}, and the WKB result for the monopole pair generation, all in a uniform magnetic background, is noticed. The difference is probably mainly due to the inadequacy of the WKB approach in obtaining the exact result. More on this and related issues remain to be studied.

\section{Heuristic derivation of the production rate in a de Sitter background}\label{sec:ds}

In this section we study black hole pair production on a cosmic string in the presence of a background magnetic field along with a positive cosmological constant $\Lambda$. The rate of black hole pair creation in dS background has been first studied in~\cite{mellor1989}, and further investigated in~\cite{mann1995}. It was shown that the action of the instanton that mediate between dS spacetime and dS spacetime with a pair of black holes, each of mass $m$, is $I = \pi m \sqrt{3/\Lambda}$. Later, in~\cite{dias2004:2}, a cosmic string was invoked in a dS background and the rate of black hole pair production has been studied.

To study the pair production of dS black holes in the presence of a background magnetic field, one needs a dS version of the Ernst metric. However such a solution does not exist~\cite{dias2003}. In fact, Ernst has used an Ehlers-Harrison type transformation~\cite{ehlers1959,harrison1968} to add the magnetic field to the C metric. But, applying this transformation to the cosmological C metric would not yield a new solution to the Einstein-Maxwell theory~\cite{dias2003}.

Nonetheless we can deal with this problem heuristically. One can think of a positive cosmological constant doing the same role as the cosmic string; in the sense that it provides the energy to materialize the pair of black holes and, also, provides the force to accelerate the black holes apart.

To find the acceleration caused by a positive cosmological constant, we recall that the Newtonian potential of dS spacetime is $\Lambda r^2/6$. The force per unit mass or the acceleration would then be $\Lambda r/3$, where $r$ is the dS radius given by $r=\sqrt{3/\Lambda}$. Therefore, in the weak field limit, we can modify the Newton's law (\ref{newton}) to be
\begin{equation}
mA=\tilde{\mu}+Bq,
\end{equation}
where $\tilde{\mu}=\mu_3-\mu_4+m\sqrt{\Lambda/3}$. Likewise, we can modify the actions (\ref{eqn:act:emp}) and (\ref{ours}) by replacing $\mu_3-\mu_4$ by $\tilde{\mu}$.

Consider, for instance the $q=m$ case in the weak field limit. The dS counterpart of the action (\ref{ours}) is
\begin{equation}
I\simeq \frac{\pi m^2 \left[1-2\left(\tilde{J}+\tilde{\mu}\right)\right]\left(1-4\mu_3\right)}{\tilde{J}-1+\left(1+2\tilde{\mu}\right)\tilde{\mu}}, \label{ours:ds}
\end{equation}
where
\begin{equation}
\tilde{J}=\sqrt{\left(1-2\tilde{\mu}\right) \left[1-2\left(\hat{B}m+\tilde{\mu}\right)\right]}. \nonumber
\end{equation}

In the limit $\Lambda\ll m^{-2}$, and $\hat{B},\, \mu_3,\, \mu_4 \rightarrow 0$, Eq. (\ref{ours:ds}) reduces to $I\simeq \pi m \sqrt{3/\Lambda}$, to leading order of $\Lambda$. This is consistent with the result of~\cite{mann1995} for the case of lukewarm black holes in which one has non-extreme black holes satisfying $q^2=m^2$~\cite{romans1992}.

We note here that if we modify the parameters (\ref{eqn:para:r})-(\ref{eqn:para:A}) to include the cosmological term and substitute them into Eq. (\ref{action}), then, the limit $\Lambda\ll m^{-2}$, and $\hat{B},\, \mu_3,\, \mu_4 \rightarrow 0$ would lead to $I\simeq \pi m \sqrt{3/\Lambda}$ as well. This is because in this limit we have $r_+\rightarrow r_-$, or equivalently $q\rightarrow m$. Therefore, again we are left with lukewarm black holes. This result shows that the lukewarm instanton is the only one available for non-extreme black holes, for which the conical singularities are eliminated at both $y=\xi_2$ and $y=\xi_3$ horizons in a dS background (see~\cite{mann1995} and~\cite{dias2004:2}).

Now let us study the action (\ref{ours:ds}) through some examples, considering the current bounds on the physical quantities. For $\hat{B}=10^{-63}$, $\mu_3=2\times 10^{-7}$, $\mu_4=10^{-7}$, $m=10$, and $\Lambda=10^{-122}$ we find $I\simeq 3.141590 \times 10^{9}$. On the other hand, if we have no string (i.e.~$\mu_3=\mu_4=0$) and keep the other parameters the same, we find $I\simeq 5.348755 \times 10^{62}$.

For the inflationary era, by taking $\hat{B}=0$, $\Lambda=10^{-12}$, $m=10$, and $\mu_3=2\times 10^{-7}$, and $\mu_4=10^{-7}$, we find $I\simeq 5.348751 \times 10^{7}$. For the no string case (i.e.~$\mu_3=\mu_4=0$) it is $I\simeq 5.441398 \times 10^{7}$, which is a bit larger. Although the value of the action for pair creation of black holes would not change significantly during inflation regardless of the presence of the cosmic string, the action is still considerably smaller in this era than the time it frays in presence of today's background magnetic field.

\section{Concluding remarks}\label{sec:con}

The rate of the black hole pair production (per unit time per unit length of the cosmic string) is given by $e^{-I}$. We see from Fig.~\ref{fig:action} that as the background magnetic field increases the black hole pair creation would be more probable. These black holes are accelerated away from each other after their production due to the force of the string and the background magnetic field.

We worked on a theoretical basis, however, we have presented our numerical examples with observationally viable values. Since the background magnetic field in the Universe is very small, we could safely use the weak field limit of the instanton action in cosmological applications. We have compared the pair creation of black holes in a background magnetic field (and no string) with that on a cosmic string in the presence of the background magnetic field. We have shown that the fraying/breaking of the cosmic string, even with the observational bound on its tension, would decrease the action significantly.

In~\cite{emparan1995}, the production rate of black hole pair, connected by a strut but with no string, was obtained in the Melvin universe. It was shown that this rate is suppressed relative to the case where no strut is present. Here we showed that if a cosmic string is present in the Melvin universe from the beginning, the pair production rate is enhanced.

In the early universe at the end of the inflationary era, the background magnetic field is expected to be much larger if they are produced during inflation. We find that the action of the pair production in no string case could be much smaller in that period. However if the pair production took place on a cosmic string, the value of the action at the end of the inflationary era and present day universe do not differ significantly.

We have studied the production rate of monopole-antimonopole pair as well. Here we have used WKB approximation to find the exponent of the rate. We have found that the rate of monopole pair production is suppressed relative to non-extremal black hole pair production by the factor $e^{S_{bh}/2}$. This is not consistent with the suggestion that monopoles behave like extremal black holes in the sense that their production is suppressed relative to non-extremal black hole production by the factor $e^{S_{bh}}$. There might be a problem with the WKB approach. There are various loop corrections with the same order as black hole entropy that have not been taken into account in WKB approximation~\cite{affleck1982,garfinkle1994}. This issue needs further investigations in future studies.

We have also presented a heuristic derivation of the production rate on cosmic string in the presence of a background magnetic field in dS spacetime. Like the cosmic string the presence of a positive cosmological constant decreases the action, hence increases the production rate. By using the action obtained in this way, we find that the pair creation is the most probable in the inflationary era. Contrary to the epoch after inflation, the black hole pair creation rate during inflation is nearly insensitive to whether it takes place on a cosmic string.

In our numerical example we have taken the mass of the black holes to be $m=10$ in Planck units. Such small black holes may tolerate the evaporation~\cite{macgibbon1987}. Besides, for the evaporation of the magnetically charged black holes, one needs copious magnetically charged monopoles to take away the charge of the black holes. This is not observationally acceptable as we do not see many (if any) monopoles around us. In GUT theories with magnetically charged monopoles with with mass $10^{16}~{\rm GeV}\lesssim \mu \lesssim  10^{19}~{\rm GeV}$, it has been shown that the black holes with mass $m\geq 150 \left(\frac{\mu}{10^{17}~{\rm GeV}}\right)^2$ can become stable roughly over the lifetime of the universe \cite{Hiscock:1983ag}. The temperature of the accelerating black holes is also somewhat subtle to define (see~\cite{appels2016,appels2017,appels2020}). More study on the thermodynamics of such black holes is needed.

\bibliographystyle{model1-num-names}


\section*{Results}
\subsection*{Fitting a Logistic Regression Model}
We used the training dataset with 351 patients without missing any values of the three important features (i.e., LDH, lymphocyte, and hs-CRP) to fit a step-wise logistic regression model with all the second-order interaction items. We used the coefficient of determination adjusted $R^2$  \cite{hastie2009elements}, indicating the portion of variance in the dependent variable (i.e., patient outcome) explained by the independent variables, to determine which interaction items to include. We obtained a model with the largest adjusted $R^2 = 0.797$ with all the three variables, i.e., LDH ($P = 0.000$), lymphocyte ($P = 0.000$), and  hs-CRP ($P = 0.000$), and two interaction items, i.e., LDH:lymphocyte ($P = 0.081$) and LDH:hs-CRP ($P = 0.000$).
\subsection*{Prediction Performance}
Table 1 shows the fatality prediction performance in terms of area under the curve (AUC) \cite{walter2005partial} scores of logistic regression using 100-round five-fold cross-validation. It can be seen that the logistic regression models are able to accurately predict the fatality outcomes of COVID-19 patients with three important biomarkers. As the number of the features increases, the performance also improves before the last interaction item (i.e.,lymphocyte:hs-CRP) was added. This was consistent with the results of fitting a logistic regression model found above. The model with optimal performance had two (product) interaction items, i.e., LDH:lymphocyte and LDH:hs-CRP.

\begin{figure*}[hbt!]
\centering
\includegraphics[width=17.8cm]{Figures/PredictTimeAhead.png}
\caption{The performance of multi-day ahead forecasting of the proposed logistic regression model: Histogram of the maximum days ahead of the correctly predicted patients of (a) all the 485 patients and (c) 110 patients in the external test set; The sample size, f1-score, accuracy in each day and cumulative f1-score and cumulative accuracy over the days of (b) all the the samples and (d) the samples in the external test set. Note there were three negative value records excluded in the figure as they were obtained after the patient was discharged or died.}
\label{fig:frog}
\end{figure*}

\subsection*{An Explainable Logistic Regression Model}
Based on the results in Table 1, we developed an explainable logistic regression model with two interaction items to predict the fatality rate of COVID-19 patients. The modeling of prediction using logistic regression models is transparent and the model can produce fatality probabilities between 0 and 1 rather than a binary value. Using $Y=1$ and $Y=0$ to indicate death and survival, respectively, we can formulate the logistic regression model in the following way \cite{hastie2009elements}:
\begin{equation} \label{eu_eqn}
l\ =\ ln\frac{P( Y=1)}{P( Y=0)} =ln\frac{p}{1-p} =\mathbf{B^{T}X} \ =\sum\limits _{i} \beta _{i} x_{i},
\end{equation}
where $\beta_{0}$ is a constant and $\beta_{i}, i=1,..., 5$, are the coefficients of $x_{1}$ (LDH), $x_{2}$ (lymphocyte), $x_{3}$ (hs-CRP), $x_{4}$ (LDH:lymphocyte), and $x_{5}$ (LDH:hs-CRP). We ran 100 rounds of five-fold cross-validation with random search to identify the optimal values of the regularization term $C$ between 0.0001 and 1000 with two types of penalty, i.e., $l1$ and $l2$. This process resulted in 500 sets of the coefficients and we used the median values as the coefficients of the final model in a vector form: $\mathbf{B^{T}} = [-4.976, 1.440e-2, -3.053e-1, 4.378e-2, 4.766e-4, -6.748e-5]$. We also optimized the prediction performance of the model by adjusting the threshold of the death probability both for the external test data and multi-day ahead forecasting below. We found the optimal threshold was 0.8. In other words, when one's fatality probability was larger than 0.8, the model predicted that a patient would die and the model had the best performance.

The model in Eq.(1) with the identified threshold was then used to predict the outcomes of 110 patients in the external test set that was not used to build the model. Although the data set was rather unbalanced (13 deaths and 97 survivals), the performance of the proposed logistic regression model was promising with 96 true negatives (1 false negative) and 12 true positives (1 false positive). The accuracy, f1-score, and AUC were $98.18\%$, $92.38\%$ and 0.996, respectively.

\subsection*{Multi-Day Ahead Forecasting}
Since there were multiple records of the three biomarkers for each patient, the model was used to forecast patient outcomes multi-days in advance. The samples were obtained by examining the records within each day for each patient. A total number of 909 records was obtained for all the 485 patients with 251 records for the 110 patients in the external test set. For multi-day ahead forecasting, we aimed to obtain the maximum days in advance. For example, if there are $N_{i}$ records in total for patient $i$ and the model was able to predict $m, m<=M$ days ahead, the predicted outcomes and the ground truth of the following maximum $n_{i}, 1<=n_{i}<=N_{i}$ consecutive records need to be the same. Here $m$/$M$ are the days between the dates of the $n_{i}$-th/first records and the date of the final outcome, respectively. For all the records, the model was able to predict 11.30 days (maximum = 34.91) ahead on average with a cumulative f1-score of 93.76\% (see Fig.1 (a) and (b)) and a cumulative accuracy value of 93.92\%. For the 251 records in the external test set, the model was able to predict 12.76 (maximum = 33.15) days ahead on average with a cumulative f1-score of 95.73\% (see Fig.1 (c) and (d)) and a cumulative accuracy value of 96.47\%. Thus, the proposed model can potentially give doctors the time needed to treat the patient accordingly.

\section*{Discussion}
Built on \cite{yan2020interpretable}, we proposed an explainable, intuitive, and yet accurate prediction model using logistic regression by incorporating two interaction items among the three most important biomarkers, including LDH, lymphocyte (\%) and hs-CRP, which can be easily measured in hospitals. Our model used no extra input information, compared to the decision tree model in \cite{yan2020interpretable}. Unlike the binary decisions produced by the decision tree, the logistic regression model produced a probability of fatality for each patient, which is more consistent with human-friendly explanations of machine learning models \cite{molnar2020interpretable}. For example, one of the rules in the decision tree in \cite{yan2020interpretable} was that IF LDH $> 365 UI^{-1}$, THEN death. Such a binary prediction may be not very intuitive without telling the likelihood of death. As a value of $364 UI^{-1}$ might result in a significantly different fatality probability from a value of $ 64 UI^{-1}$.

Our model, on the other hand, always gave a probability of death and also identified a threshold probability at 0.8, above which the model predicted that the patient would die. Furthermore, our model also outperformed the decision tree model in terms of average maximum days to the outcome and the cumulative f1-score and accuracy in Fig. 1.
Such a model can offer the clinicians time to identify high-risk patients before they become critical. Thus, an appropriate treatment strategy for COVID-19 patients depending on their likelihood of death can be made using corresponding medical resources. This can potentially alleviate the shortages of critical medical resources in hospitals in the current situation.

However, the model was built on a relatively small sample size. More research is needed to further test and optimize the model, taking both explainability and prediction performance into account.

\section*{Materials and Methods}
\subsection*{Samples} The data was originally from \cite{yan2020interpretable}. The model construction was based on the data of 375 (174 died) patients collected between January 10, 2020 and February 18, 2020 from Tongji Hospitals, Wuhan, China. Of them, 24 of the patients had missing data in the three biomarkers and thus were excluded from analysis. The external test data set was collected from another 110 patients (13 died) between February 19, 2020 and February 24, 2020 from the same hospital. We reported the performance of the model using metrics, including AUC \cite{walter2005partial}, micro-ave f1-score  \cite{zhou2011affect}, and accuracy \cite{zhou2020driver}.
\subsection*{Data Analysis}
Interaction items in logistic regression can potentially improve the performance of the model to a great extent \cite{levy2019don}. Hence, we first fitted a logistic regression model using the three most important biomarkers, i.e., LDH, lymphocyte (\%), and hs-CRP identified in \cite{yan2020interpretable} and further identified two interaction items that could be useful to improve the prediction model. Then, we added one item at a time to the logistic regression model with five-fold cross-validation for 100 rounds and verified the two identified interaction items did improve its prediction performance as shown in Table 1. In order to have a model with good generalizability, we used the median values of the coefficients produced from the 500 models when producing the results in Table 1. Finally, this model was used to predict the external test set and multi-day ahead forecasting.

\section*{Conclusions}
The proposed logistic regression model can effectively predict the outcomes of COVID-19 patients with fatality probabilities. The model is accurate, intuitive, and explainable with only three blood biomarkers and two of their interaction items as input, which can potentially help the doctors determine the best treatment route for COVID-19 patients with high risks and optimize the logistic planning in the medical systems around the world amid this COVID-19 pandemic.

\acknow{We thank the authors in \cite{yan2020interpretable} for
providing their data sets.}

\showacknow{}


\section{Introduction}

Control Barrier Functions (CBFs) have proven to be effective at enforcing safety in nonlinear systems. The goal of this paper is to introduce a new form of CBFs that are suitable for feedback systems where the controller is defined  by a differential equation \cite{Isidori90}, i.e., integral control \cite{freeman1995robust,khalil2000universal,jiang2001robust}.
The motivation comes from a dynamically defined tracking controller that has displayed effective tracking convergence \cite{Wardi19}, but also can exhibit large overshoots in the input controls at early transient phases. We investigate a CBF-based approach to limit these overshoots, i.e., satisfy input bounds, while also satisfying safety constraints on the state.  The results are presented in a more-general setting of pointwise constraints on input and state encoded by a novel form of CBF.

The general CBF framework (see \cite{Ames19} and references therein) requires solving an optimization problem at each time to enforce the CBF condition. If the dynamics of the state are control affine, this optimization problem is a convex quadratic program (QP), solvable in real time. Yet this real-time solvability hinges on the control affine nature of the dynamics. Moreover, when input constraints are also present, adding them to the QP can result in its infeasibility.  This paper will address both of these existing limitations in the context of integral control.

This paper introduces a new class of CBFs: \emph{integral control barrier functions (I-CBFs)}.  These barrier functions are developed in the context of systems with dynamically defined controllers, i.e., systems where the evolution of the state and input are described by an ordinary differential equation.  As a result, I-CBFs are defined on both the state and input, allowing for the inclusion of the input in the safety conditions encoded by this function. A related formulation is presented in \cite{huang2019guaranteed}, where control-dependent CBFs are defined: these, unlike I-CBFs, are considered for inputs with bounded time-derivative and, importantly, the integration with nominal tracking controllers was not considered. In this paper, given CBFs and I-CBFs, we present a resulting controller that guarantees safety in state and input, while minimally modifying a nominal dynamically defined controller. Additionally, we explicitly explore systems with both state and input constraints, and demonstrate how these can be unified through the I-CBF framework.

\section{Background Material}
\label{sec:background}

This section summarizes established results regarding safety-critical control via control barrier functions (CBFs), and the tracking-control technique that will inspire the main result of this paper.
Herein, we consider a nonlinear control system defined by the differential equation:
\begin{equation}
\label{eqn:controlsys}
    \dot{x}(t)=f(x(t),u(t)),
\end{equation}
where $t\geq 0$, the state is $x(t)\in \mathbb{R}^n$, the input is  $u(t)\in \mathbb{R}^m$, and an initial state $x_{0}:=x(0)\in \mathbb{R}^n$ is given.
Assume that the function $f:\mathbb{R}^n\times \mathbb{R}^m\rightarrow \mathbb{R}^n$ is continuously differentiable,
and it
satisfies sufficient conditions for the existence of a unique solution to \eqref{eqn:controlsys} over all $t\geq 0$ for every bounded, piecewise-continuous control $u(t)$ and $x_{0}\in \mathbb{R}^n$.

\subsection{Control Barrier Functions}
\label{ssec:cbfs}

Control Barrier Functions (CBFs) ensure the forward invariance of a set $\S\subset\mathbb{R}^n$, in which case the system is considered \emph{safe}. The CBF framework was introduced in \cite{Ames14,Ames17} where it was applied to adaptive cruise control, and has since been applied to a variety of application domains, including: automotive safety \cite{jankovic2018robust,xu2015robustness}, robotics \cite{nguyen2016exponential,wang2016multi} and multi-robot systems \cite{Wang17,lindemann2018control}.
See \cite{Ames19} for a recent survey.
Note that alternative approaches to CBFs include reference governors \cite{gilbert2002nonlinear,nicotra2018explicit} and integral barrier Lyapunov functions \cite{tee2012control}; yet, the former is computationally more expensive, whereas the latter only applies to systems in strict-feedback form and does not explicitly consider input bounds.

Let $h:\mathbb{R}^n\rightarrow \mathbb{R}$ be a continuously-differentiable function
such that $0$ is a regular value.  We will consider the set, $\S\subset\mathbb{R}^n$, given as the 0-superlevel set of $h$, i.e., $\S$ is defined by:
\begin{eqnarray}
\S & := & \{ x \in \mathbb{R}^n ~ : ~ h(x) \geq 0 \} , \nonumber\\
\partial \S & = & \{ x \in \mathbb{R}^n ~ : ~ h(x) = 0 \}, \nonumber\\
\mathrm{Int}(\S) & = & \{ x \in \mathbb{R}^n ~ : ~ h(x) > 0 \}. \nonumber
\end{eqnarray}
Consider a feedback control law applied to \eqref{eqn:controlsys} of the form:
\begin{equation}
\label{eqn:feedbackcont}
    u(t)=k(x(t)),
\end{equation}
with $k:\mathbb{R}^n\rightarrow \mathbb{R}^m$ continuous.  Under this controller, we say that $x(t)$ is {\it safe} if and only if
$x_0 \in \S$ implies that $x(t)\in \S$ for all $t \geq 0$, i.e., $\S$ is \emph{forward invariant}, i.e., $\S$ is safe.

The function $h$ is a {\it control barrier function (CBF)} if  there exists a continuous extended class-${\cal K}$ function $\gamma:\mathbb{R}\rightarrow \mathbb{R}$ (monotone increasing, $\gamma(0)=0$) together with a feedback controller \eqref{eqn:feedbackcont} such that along every trajectory of the closed-loop system,  the following CBF condition holds \cite{Ames17}:
\begin{equation}
\label{eqn:CBFcond}
    \frac{d}{dt}h(x(t))+\gamma(h(x(t)))\geq 0.
\end{equation}
Additionally, if Eq. \eqref{eqn:CBFcond} is satisfied for all $t\geq 0$, then the set $\S$ is forward invariant and asymptotically stable \cite{Ames17}.
Now by \eqref{eqn:controlsys}, it follows that $\dot{h}$ has the following form:
\begin{equation}
\label{eqn:cbfcondcontrol}
    \dot{h}(x(t),u(t)) =\frac{d}{dt}h(x(t))=\frac{\partial h}{\partial x}(x(t))f(x(t),u(t)).
\end{equation}
Therefore, the closed-loop system is safe if for every $t\geq 0$,
\begin{equation}
    \frac{\partial h}{\partial x}(x(t))f(x(t),u(t))+\gamma(h(x(t)))\geq 0.
\end{equation}
If this inequality is not satisfied for a particular $t\geq 0$, then the control law can be modified as follows to guarantee the safety of $\S$. Given $x\in \mathbb{R}^n$, define the set
$K_{x}\subset \mathbb{R}^m$ by:
\begin{align}
K_{x}= & \left\{u\in \mathbb{R}^m :  \frac{\partial h}{\partial x}(x)f(x,u)+\gamma(h(x))\geq 0\large \right\}.
\end{align}
Therefore, $u(t)$ given by:
\begin{equation}
\label{eqn:optimizationbasedcontrol}
    u(t)\in \underset{u \in \mathbb{R}^{m}}{\operatorname{argmin}} \{||u-k(x(t))||^2~:~u\in K_{x(t)}\}
\end{equation}
ensures the safety of the closed-loop  system while modifying the control law defined in \eqref{eqn:feedbackcont} in a minimally invasive fashion.

Finally, note that if the dynamics of the plant are control-affine, namely
$f(x,u)=f_{0}(x)+f_{1}(x)u$ for  functions $f_{0}:\mathbb{R}^n\rightarrow \mathbb{R}^n$
and $f_{1}:\mathbb{R}^n\rightarrow \mathbb{R}^{n\times m}$, then $u(t)$, defined by \eqref{eqn:optimizationbasedcontrol},  can be computed by a quadratic program (QP).  In particular, we obtain an example of the feedback controller $u = k_{\mathrm{s}}(x)$ that ``filters'' the controller in \eqref{eqn:feedbackcont} and renders the system safe:
\begin{align}
\label{eqn:CBFQP}
k_{\mathrm{s}}(x) = \underset{u \in \mathbb{R}^{m}}{\operatorname{argmin}} & ~ ||u-k(x)||^2 \\
\mathrm{s.t.} &  ~ \frac{\partial h}{\partial x}(x) f_0(x) +
\frac{\partial h}{\partial x}(x) f_1(x) u \geq - \gamma(h(x)), \nonumber
\end{align}
where we suppressed the dependence on $t$.
This QP admits a real-time implementation (see \cite{Ames14}).
Importantly, this paper will allow for quadratic program representations of safety-critical CBF-based controllers even when the plant is not control affine.

\subsection{Tracking Control by a Newton-Raphson Flow}
In this paper, we consider a tracking-control technique based on a flow version of the Newton-Raphson method for solving algebraic equations \cite{Wardi19}. This provides us with a convenient expression of a feedback control law defined by an ordinary differential equation, as it is recalled in this section (a more-detailed discussions and analyses can be found in \cite{Wardi19}). It is worth noting, however, that the main results presented in this paper do not depend on this specific method employed.  We consider a specific method to define an integral control law that achieves tracking for definiteness.

To explain the main idea behind the tracking controller utilized, consider first the simple case where the plant is a memoryless nonlinearity, input-output system of the form $y(t)=g(u(t))$,
where $u(t)\in \mathbb{R}^m$ is the input, $y(t)\in \mathbb{R}^m$ is the output, and $g:\mathbb{R}^m\rightarrow \mathbb{R}^m$ is a continuously-differentiable function.
Given a continuously-differentiable reference signal $r(t)\in \mathbb{R}^m$,
suppose that the objective is to design a control law such that $y(t)$ converges to $r(t)$ in a suitable sense as described below.  Consider the dynamically defined controller:
\begin{equation}
\label{eqn:udotnomem}
    \dot{u}(t)=\Big(\frac{\partial g}{\partial u}(u(t))\Big)^{-1}\big(r(t)-y(t)\big).
\end{equation}
This controller essentially implements a Newton-Raphson (NR) flow whose vector field at time $t$, $\dot{u}(t)$, is the direction defined by the classical NR method for
solving the  time-dependent algebraic equation: $r(t)-g(u)=0$.

Next, suppose that the plant is dynamic, i.e., defined by the state equation \eqref{eqn:controlsys}, and consider the output equation $y(t)=\zeta(x(t))$,
where the function $\zeta:\mathbb{R}^n\rightarrow \mathbb{R}^m$ is continuously differentiable. Now $x(t)$, and hence $y(t)$, are functions of the initial condition $x_0$ and past control actions: $u(\tau),~\tau\in[0,t)$.  Yet, instantaneously, $y(t)$ is not a function of $u(t)$.
Therefore the controller $u(t)$ cannot be defined by an equation like \eqref{eqn:udotnomem}.
However, given $t\geq 0$ and $T> 0$, $x(t+T)$ and hence $y(t+T)$ depend on $u(\tau),~\tau\in[t,t+T]$.  This observation forms the basis for our tracking controller.

Given a time $t \geq 0$ and $T > 0$, fixing $u_t := u(t)$ as a constant over $[t,t + T]$ and (approximately) forward integrating \eqref{eqn:controlsys}, i.e., $\dot{\hat{x}}(\tau) = f(\hat{x}(\tau),u_t)$, over this interval with initial condition $\hat{x}(t)= x(t)$ results in prediction of the state $\hat{x}(t + T)$ and hence a prediction of the output:
\begin{equation}
\label{eqn:predictor}
    \hat{y}(t+T)= \zeta(\hat{x}(t + T)) =: g(x(t),u(t))
\end{equation}
that, therefore, depends on $x(t)$ and $u(t)$.
In this case, a suitable extension of \eqref{eqn:udotnomem} has the following form:
\begin{equation}
\label{eqn:controleqn}
    \dot{u}(t)=\alpha\Big(\frac{\partial g}{\partial u}(x(t),u(t)\Big)^{-1}\big(r(t)-\hat{y}(t+T)\big);
\end{equation}
where $\alpha > 0$ is a controller gain.
In \cite{Wardi19}, the following convergence result was established:
\begin{equation}
\underset{t\rightarrow\infty}{\operatorname{limsup}}
||r(t)-y(t)||<\eta_{1}+\eta_{2}/\alpha,
\end{equation}
where $\eta_{1}:={\rm limsup}||y(t)-\hat{y}(t)||$, and $\eta_{2}:={\rm limsup}||\dot{r}(t)||$, where limsup means as $t\rightarrow\infty$. Thus, increasing the controller gain, $\alpha$, can reduce the error due variations in $r(t)$, but cannot attenuate the effects of prediction errors.

\section{Integral Control Barrier Functions}
\label{sec:mainresult}

This section presents the main construction and results of the paper.  Specifically, we introduce integral control barrier functions, and demonstrate that they can enforce safety for dynamically defined control laws.  This, in essence, creates a new paradigm for safety-critical integral control.

Consider a dynamical system defined by Eq. \eqref{eqn:controlsys} with an initial condition $x_{0}:=x(0)\in \mathbb{R}^n$. In a departure from the type of controller considered in Section II.A, which is algebraic, we define a general feedback law by the ordinary differential equation:
\begin{equation}
\label{eqn:dyncontrolequation}
  \dot{u}(t)=\phi(x(t),u(t),t),
\end{equation}
with an initial condition $u_{0}:=u(0)\in \mathbb{R}^m$.
We assume that the  function
$\phi:\mathbb{R}^n\times \mathbb{R}^m\times \mathbb{R}\rightarrow \mathbb{R}^m$ is  continuously differentiable. The closed-loop system is defined by
the state equation \eqref{eqn:controlsys} and the control equation \eqref{eqn:dyncontrolequation}.
We  write these two equations jointly in the following way:
\begin{equation}
\label{eqn:stateeqn}
  \left[\begin{array}{l}
  \dot{x}(t)\\
  \dot{u}(t)
  \end{array}
  \right]~=~
  \left[
  \begin{array}{l}
  f(x(t),u(t))\\
  \phi(x(t),u(t),t)
  \end{array}
  \right].
\end{equation}
Define $z(t):=(x(t)^{\top},u(t)^{\top})^{\top}\in \mathbb{R}^n\times \mathbb{R}^m$ with $z(t)$ the {\it augmented state}, which can be viewed as the state of the closed-loop system.  Note that the state equation \eqref{eqn:stateeqn} has no external input.

 Let $\S\subset \mathbb{R}^n\times \mathbb{R}^m$ be a set, the {\it safety set}, defined as the 0-superlevel set of a continuously differentiable function $h:\mathbb{R}^{n}\times \mathbb{R}^m\rightarrow \mathbb{R}$ with $0$ a regular value. In addition, suppose the existence of a continuous extended class-${\cal K}$ function $\gamma:\mathbb{R}\rightarrow \mathbb{R}$ such that the following is satisfied along every trajectory $z(t)$ of the closed-loop system:
 \begin{equation}
 \label{eqn:zCBFcond}
     \frac{d}{dt}h(z(t))+\gamma(h(z(t)))\geq 0.
 \end{equation}
Then, the set $\S$ is forward invariant and asymptotically stable.
 Throughout the remainder of this paper, for the sake of simplicity of exposition, we will suppress the dependence on time---this can easily be inferred from context.

 For notational simplicity, define the following vector valued $p(x,u) \in \mathbb{R}^m$ and scalar valued $d(x,u) \in \mathbb{R}$ functions:
  For given $t\geq 0$, define
  the scalar $d $ and vector $p \in \mathbb{R}^m$ by:
 \begin{align}
 \label{eqn:pxu}
      p(x,u)  :=  &  \Big(\frac{\partial h}{\partial u}(x ,u )\Big)^{\top} \\
      \label{eqn:pdxut}
     d(x,u,t)   := &  -\Big(\frac{\partial h}{\partial x}(x ,u )f(x ,u )  +\frac{\partial h}{\partial u}(x ,u )\phi(x ,u ,t) \nonumber\\
     & \qquad \qquad  +\gamma(h(x ,u ))\Big),
 \end{align}
Note that Eq. \eqref{eqn:zCBFcond} can thus be recast as $d(x,u,t) \leq 0$.

\newsec{Main Result.}
Traditional CBF methods cannot be directly applied to systems of the form given in \eqref{eqn:stateeqn} due to the dynamically defined control law, i.e., the control law $u$ is a result of integrating the augmented state equation in \eqref{eqn:stateeqn}.
Our approach is to add an auxiliary input to the state equation so that, if the inequality in \eqref{eqn:zCBFcond} is not satisfied for a particular $t\geq 0$, we are able to determine the minimal modification of the dynamically defined control law that will guarantee safety. To this end, we modify the augmented-state equation
\eqref{eqn:stateeqn} to include the input $v\in \mathbb{R}^m$ as follows:
\begin{equation}
\label{eqn:stateeqnv}
  \left[\begin{array}{l}
  \dot{x}\\
  \dot{u}
  \end{array}
  \right]~=~
  \left[
  \begin{array}{l}
  f(x,u)\\
  \phi(x ,u ,t)+v
  \end{array}
  \right].
\end{equation}

For every $z:=(x^{\top},u^{\top})^{\top}\in \mathbb{R}^n\times \mathbb{R}^m$
and $t\geq 0$, define:
\begin{align}
\label{eqn:Szset}
  K_{z,t}:= & \big\{v\in \mathbb{R}^m~
  :~\frac{\partial h}{\partial x}(x,u)f(x,u)    \\
   &+\frac{\partial h}{\partial u}(x,u)\big(\phi(x,u,t)+v\big) ~ +\gamma(h(x,u))~\geq~0\big\} \nonumber\\
  \label{eqn:Szsetshort}
   =& \left\{ v\in \mathbb{R}^m~
  :~ p(x,u)^{\top} v \geq d(x,u,t) \right\}.
\end{align}
Therefore, if $v(x,t) \in K_{z,t}$, it follows that Eq. \eqref{eqn:zCBFcond} is satisfied.  Additionally, if the dynamic control law $\dot{u} = \phi(x,u,t)$ is inherently safe, then $d(x,u,t) \leq 0$ and $v = 0$ imply that the system is safe.
This leads to the formulation of I-CBFs,
mirroring the classic definition of control Lyapunov functions \cite{sontag1989universal}) and similar in spirit to control-dependent CBFs \cite{huang2019guaranteed}.

\begin{definition}
For the system \eqref{eqn:stateeqnv}, with corresponding safe set $\S \subset \mathbb{R}^n \times \mathbb{R}^m$ defined as the 0-superlevel set of a function $h : \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}$ with $0$ a regular value: $\S = \{(x,u) \in \mathbb{R}^n \times \mathbb{R}^m ~ : ~ h(x,u) \geq 0\}$.  Then $h$ is an \emph{integral control barrier function (I-CBF)} if for any $(x,u) \in \mathbb{R}^n \times \mathbb{R}^m$ and $t \geq 0$:
\begin{eqnarray}
\label{eqn:relativedegree}
p(x,u) = 0 \qquad \Rightarrow \qquad d(x,u,t) \leq 0.
\end{eqnarray}
\end{definition}

We now have the necessary constructions to state the main result of this paper.

\begin{theorem}
\label{thm:main}
 Consider the control system $\dot{x} = f(x,u)$, with $x \in \mathbb{R}^n$ and $u \in \mathbb{R}^m$, and suppose that there is a corresponding dynamically defined controller: $\dot{u} = \phi(x,u,t)$.  If the safe set $\S \subset \mathbb{R}^n \times \mathbb{R}^m$ is defined by an integral control barrier function, $h : \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}$,
 then modifying the dynamically defined controller to be of the form:
 \begin{eqnarray}
 \label{eqn:umodifiedopt}
 \dot{u}  =  \phi(x,u,t) + v^*(x,u,t)
 \end{eqnarray}
 with $v^*$ the solution to the QP:
 \begin{align}
\label{eqn:CBFQPv}
 \qquad \qquad \quad  v^*(x,u,t)  =  \underset{v \in \mathbb{R}^{m}}{\operatorname{argmin}} & ~ ||v||^2 \\
\mathrm{s.t.} &  ~ p(x,u)^{\top} v \geq d(x,u,t) \nonumber
\end{align}
 results in safety, i.e., the control system $\dot{x} = f(x,u)$ with the dynamically defined controller \eqref{eqn:umodifiedopt} results in $\S$ being forward invariant: if $(x(0),u(0)) \in \S$ then $(x(t),u(t)) \in \S$ for all $t \geq 0$.
\end{theorem}

\begin{proof}
We only need to verify that Eq. \eqref{eqn:zCBFcond} is satisfied for $z(t) = (x(t),u(t))$ the solution to \eqref{eqn:stateeqnv} with $v = v^*(x,u,t)$.  Per \eqref{eqn:Szset} and \eqref{eqn:Szsetshort}, this is equivalent to:
\begin{eqnarray}
\label{eqn:proofmaineqn}
p(x,u)^{\top} v^*(x,u,t) \geq d(x,u,t).
\end{eqnarray}
This, in turn, will imply that $(x(t),u(t)) \in \S$ for all $t \geq 0$ if $(x(0),u(0)) \in \S$ per the main result of \cite{Ames19}.

Thus, to establish the result, we must show the solution to \eqref{eqn:CBFQPv} satisfies \eqref{eqn:proofmaineqn}.  This follows from the fact that we can obtain an explicit solution to \eqref{eqn:CBFQPv}.  In particular, the condition that $h$ is an integral control barrier function and, in particular, that it satisfies \eqref{eqn:relativedegree}, implies that the linear independent constraint qualification condition is satisfied \cite{boyd2004convex}.  Therefore, using the KKT optimality conditions (see \cite{xu2015robustness}), the solution to \eqref{eqn:CBFQPv} is given by the \emph{min-norm} controller:
\begin{align}
\label{eqn:minnormcon}
v^*(x,u,t) = \left\{
\begin{array}{lcr}
\frac{d(x,u,t)}{\| p(x,u) \|^2} p(x,u) & \mathrm{if~} d(x,u,t) > 0 \\
0 & \mathrm{if~} d(x,u,t) \leq 0
\end{array}
\right.
\end{align}
This controller is well-defined because, as $h$ is a control barrier function and thus satisfies \eqref{eqn:relativedegree}, it follows that: $d(x,u,t) > 0$ implies that $p(x,u) \neq 0$.  Additionally, from this explicit form, one can verify the Lipschitz continuity of this controller (assuming Lipschitz continuity of $p$ and $d$).

Finally, the explicit form of \eqref{eqn:minnormcon} makes it clear that it satisfies \eqref{eqn:proofmaineqn}.  If $d(x,u,t) \leq 0$ than it is naturally satisfied with $v^*(x,u,t) = 0$.  If $d(x,u,t) > 0$, then:
\begin{eqnarray}
\quad p(x,u)^{\top} v^*(x,u,t)  & = &   \underbrace{p(x,u)^{\top} p(x,u)}_{\| p(x,u) \|^2} \frac{d(x,u,t)}{\| p(x,u) \|^2} \nonumber\\
& = &  d(x,u,t). \hspace{2.6cm}  \hfill \qedhere \nonumber
\end{eqnarray}
\end{proof}

  \begin{remark}
  A natural consequence of this form of control barrier functions is, as its name suggests, its application to integral control.  In the case of ``pure'' integral control, we can consider a control system \eqref{eqn:stateeqnv} together with $\phi(x,u,t) \equiv 0$, i.e., the dynamic extension via the addition of an integrator.  In this case, the controller from Theorem \ref{thm:main} is just $u(t)~=~\int_{0}^t v^*(x(\tau),u(\tau),\tau) d \tau$.
  \end{remark}

\begin{remark}
The controller given in Theorem \ref{thm:main} can be viewed in the following fashion.  We began with a control system $\dot{x} = f(x,u)$ for which we synthesized a dynamic controller:
\begin{align*}
\dot{u} = \phi(x,u,t) +  \left\{
\begin{array}{lcr}
\frac{d(x,u,t)}{\| p(x,u) \|^2} p(x,u) & \mathrm{if} & d(x,u,t) > 0 \\
0 & \mathrm{if} &  d(x,u,t) \leq 0
\end{array}
\right.
\end{align*}
with $\phi(x,u,t)$ the ``feedforward'' integral controller that is modulated by the additional term to ensure safety through a minimal modification.  The controller, $\phi$, for example, can be given by \eqref{eqn:controleqn}.  The advantages of the integral instantiation of control barrier functions are: they allow CBFs to be applied to systems not in control affine form; the dynamic equation describing $u$ is integrated thus smoothing out the non-smooth nature of solutions to QPs; they can encode input bounds.
\end{remark}

\begin{remark}
\label{rmk:staticfeedback}
 We can also consider the case when we have a nominal controller $u = k(x)$ as in \eqref{eqn:feedbackcont}.  In this case, if we pick the following ``feedforward'' integral controller:
 $$
 \phi(x,u,t) = \frac{\partial k}{\partial x}(x) f(x,u) + \frac{\alpha}{2}\left(k(x,t)-u\right),
 $$
 for $\alpha > 0$, we provably get that $\| u - k(x) \| \to 0$ exponentially.  This can, therefore, be coupled with Theorem \ref{thm:main} to achieve provable tracking of the desired controller subject to safety.  Proving this is not within the scope of this paper, but will be the subject of future work---as it indicates the ability to safely track desired controllers even for non-affine control systems.
\end{remark}

\section{Applications and Extensions of I-CBFs}

In this section, we explore some extensions and applications of integral control barrier functions.  In particular, we demonstrate that I-CBFs can encode both input bounds and state constraints.  Importantly, we show that I-CBFs can simultaneously enforce both, thereby enabling input bounded provably safe integral control subject to feasibility.

\subsection{Application to Input Constrained Systems}
Consider the case when we have a system $\dot{x} = f(x,u)$ subject to input constraints: for $u_{\max} > 0$,
$$
\| u \|^2 \leq u_{\max} \qquad \Rightarrow \qquad h_u(u) := u_{\max} - u^{\top} u \geq 0.
$$
Naturally, this is a conservative method for enforcing input constraints due to its scalar instantiation.  Yet, it allows us to view input bounds as safety relative to the set $\S_u := \{ (x,u) \in \mathbb{R}^n \times \mathbb{R}^m ~ : ~ h_u(u) \geq 0\}$.
In this case, \eqref{eqn:pxu} and \eqref{eqn:pdxut} become:
 \begin{eqnarray}
 \label{eqn:pxu2}
      p_u(x,u)  &  =  &  - 2 u \nonumber\\
      \label{eqn:pdxut2}
     d_u(x,u,t)  &  =  &   2 u^{\top} \phi(x ,u ,t) - \gamma_u(h_u(u )). \nonumber
 \end{eqnarray}
Thus, following from Theorem \ref{thm:main}:

\begin{lemma}
\label{lem:torquebound}
For the system \eqref{eqn:stateeqnv}, $h_u$ is an integral control barrier function, i.e., $\S_u$ can be rendered safe (forward invariant), with the integral control law \eqref{eqn:umodifiedopt} where:
 \begin{align}
\label{eqn:CBFQPvubound}
  v^*(x,u,t)  =  \underset{v \in \mathbb{R}^{m}}{\operatorname{argmin}} & ~ ||v||^2 \\
\mathrm{s.t.} &  ~ 2 u^{\top} v \geq - 2 u^{\top} \phi(x ,u ,t) - \gamma_u(h_u(u ))
\nonumber
\end{align}
\end{lemma}

\begin{proof}
We only need to verify \eqref{eqn:relativedegree}.  If $p_u(x,u)= 2u = 0$, then $u = 0$, and therefore $d_u(x,u,t) = - \gamma_u(h_u(u))$.  But $u = 0$ implies that $(x,u) \in \S_u$, and therefore $h_u(u) \geq 0$ or $\gamma_u(h_u(u)) \geq 0$, and so $d_u(x,u,t) \leq 0$ as desired.
\end{proof}

\subsection{Application to State Constrained Control Affine Systems}
\label{ssec:affine}

Returning to the systems that are often considered in the context of control barrier functions---affine control systems---we wish to understand the relationship between classical and integral control barrier functions for systems of this form.
Consider the variation of \eqref{eqn:stateeqnv} in control affine form:
\begin{eqnarray}
\label{eqn:stateeqnaffine}
\dot{x} & = & f_0(x) + f_1(x) (\mu +  u) \\
\dot{u} & = & \phi(x,u,t) \nonumber
\end{eqnarray}
where here we added an auxiliary control law in the $x$ dynamics---the reason for this will become apparent soon.

For this system suppose we have a safe set defined in terms of state (and not input), i.e., by $h_x : \mathbb{R}^n \to \mathbb{R}$ with:
$$
\S_x := \{ (x,u) \in \mathbb{R}^n \times \mathbb{R}^m ~ : ~ h_x(x) \geq 0\}.
$$
Additionally, suppose that $h_x$ is a valid control barrier functions for the affine dynamics without the corresponding integrator: $\dot{h}_x(x,u) \geq -\gamma_x(h_x(x))$.
This can be quantified by a conditions analogous to \eqref{eqn:relativedegree}, specifically:
\begin{align}
\label{eqn:hxCBFcond}
\underbrace{\frac{\partial h_x}{\partial x}(x) f_1(x)}_{:= p_x(x)^{\top}} = 0
\Rightarrow
\underbrace{
\frac{\partial h_x}{\partial x}(x) f_0(x)  +  \gamma_x(h_x(x))   }_{ := -d_x(x)}  \geq 0
\end{align}
The following lemma gives a controller that will guarantee the forward invariance of $\S_x$ under the assumption that $h_x$ is a valid CBF. It will do so through a combined ``traditional'' CBF controller together with an integral controller.

\begin{lemma}
\label{lem:statebound}
Given the system \eqref{eqn:stateeqnaffine} with $h_x : \mathbb{R}^n \to \mathbb{R}$ a control barrier function for $\dot{x} = f_0(x) + f_1(x) u$.  Then the controller:
$$
k_{\mathrm{s}}(x,u,t) = \mu^*(x(t),u(t)) + \int_{0}^{t} \phi(x(\tau),u(\tau),\tau) d \tau
$$
where, for any nominal controller $u = k(x)$:
\begin{align}
\label{eqn:CBFQPint}
\mu^*(x,u) = \underset{\mu \in \mathbb{R}^{m}}{\operatorname{argmin}}&   ~ ||\mu +u -k(x)||^2 \\
\mathrm{s.t.}  ~ &  ~
p_x(x)^{\top} \mu  \geq d_x(x) - p_x(x)^{\top} u \nonumber
\end{align}
renders the set $\S_x$ forward invariant, i.e., safe.
\end{lemma}

\begin{proof}
By \eqref{eqn:hxCBFcond}, the QP \eqref{eqn:CBFQPint} is well-defined and has a solution.  This follows from the fact that $u$ not being a decision variable anymore does not affect \eqref{eqn:hxCBFcond}, i.e., $p_x(x) = 0$ still implies that $d_x(x) \leq 0$ since the term $-p_x(x) u$ vanishes.  Thus safety, i.e., forward invariance of $\S_x$, is guarenteed by the classic CBF result \cite{Ames17}.
\end{proof}

\subsection{Extension to Multiple Control Barrier Functions}

With the goal of now simultaneously enforcing state and input constraints in a holistic fashion,
we will begin by taking inspiration from Section \ref{ssec:affine}, wherein we must achieve similar results for systems with only integral control.
Returning to original state equations \eqref{eqn:stateeqnv}, but in affine form:
\begin{eqnarray}
\label{eqn:stateeqnaffinetwop}
\dot{x} & = & f_0(x) + f_1(x) u \\
\dot{u} & = & \phi(x,u,t) + v,  \nonumber
\end{eqnarray}
consider state constraints encoded by $\S_x$, i.e., by $h_x(x) \geq 0$. The input $v$ no longer appears in $\dot{h}_x$, but, the forward invariance of the set $\S_x$ is quantified by the condition $\dot{h}_x(x,u)  + \gamma_x(h_x(x)) \geq 0 $. Thus, in order to ensure the latter inequality is satisfied, we follow the approach in \cite{nguyen2016exponential} and let:
 \begin{eqnarray}
 \label{eqn:heeqn}
h_e(x,u)   :=  \dot{h}_x(x,u)  + \gamma_x(h_x(x)) =  p_x(x)^{\top} u - d_x(x).
\end{eqnarray}
The result is a function: $h_e : \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}$ with corresponding set: $\S_e = \{ (x,u) \in \mathbb{R}^n \times \mathbb{R}^m ~ : ~ h_e(x,u) \geq 0\}$.
Let $p_e(x,u)$ and $d_e(x,u,t)$ be the corresponding functions defined from $h_e$, for some $\gamma_e$, as in \eqref{eqn:pxu} and \eqref{eqn:pdxut}.

\begin{lemma}
\label{lem:cbftoicbf}
Given the system \eqref{eqn:stateeqnaffinetwop} with $h_x : \mathbb{R}^n \to \mathbb{R}$ a control barrier function with $p_x(x) \neq 0$, i.e., $h_x$ has relative degree 1, then $h_e(x,u) = \dot{h}_x(x,u)  + \gamma_x(h(x))$ is an integral control barrier function.  Additionally, for the I-CBF control law \eqref{eqn:CBFQPv}, using $p_e$ and $d_e$ in place of $p$ and $d$, if:
$$
\begin{array}{r}
x_0 \in \S_x  \\
(x_0,u_0) \in \S_e
\end{array}
\quad \Rightarrow \quad
\begin{array}{rr}
x(t) \in \S_x  & \quad \forall~ t \geq 0 \\
(x(t),u(t)) \in \S_e &  \quad \forall~ t \geq 0
\end{array}
$$
That is, $\S_x$ is safe subject to appropriate initial conditions.
\end{lemma}

\begin{proof}
It is easy to verify that $p_e(x,u) = p_x(x)$, thus the relative degree condition transfers from $h$ to $h_e$.  As a result, the controller in \eqref{eqn:CBFQPv} is well-defined, i.e., \eqref{eqn:relativedegree} is trivially satisfied.  As a result, $\S_e$ is forward invariant, i.e., assuming an initial condition $(x_0,u_0) \in \S_e$ then $\dot{h}_e(x(t),u(t)) = \dot{h}_x(x(t),u(t)) + \gamma(h_x(x(t))) \geq 0$.  Coupling this with the assumption that $x_0 \in \S_x$, i.e, that $h_x(x_0) \geq 0$, it follows that $\S_x$ is also forward invariant.
\end{proof}

\newsec{Combining State Constraints with Input Bounds.} With the I-CBF $h_e$, obtained from $h_x$ via Lemma \ref{lem:cbftoicbf}, we can synthesize controllers that
enforce state constraints with input bounds.

\begin{theorem}
\label{thm:maintwo}
Consider the system in \eqref{eqn:stateeqnaffinetwop}, together with corresponding CBF $h_x : \mathbb{R}^n \to \mathbb{R}$ and I-CBF $h_u : \mathbb{R}^m \to \mathbb{R}$, and safe sets $\S_x$ and $\S_u$ respectively, and let $h_e : \mathbb{R}^n \times \mathbb{R}^m$ be given as in \eqref{eqn:heeqn}.   If the following QP is feasible:
\begin{align}
\label{eqn:CBFQPintfinal}
v^*(x,u,t)
= &\underset{v\in \mathbb{R}^{m}}{\operatorname{argmin}}   ~  ||v||^2 \\
& \quad \mathrm{s.t.}  \quad
\begin{bmatrix}
p_x(x)^{\top}  \\
p_u(x,u)^{\top}
\end{bmatrix}
v \geq
\begin{bmatrix}
d_e(x,u,t) \\
d_u(x,u,t)
\end{bmatrix}
\nonumber
\end{align}
then the dynamically defined (integral) controller
$$\dot{u}~=~\phi(x,u,t) + v^*(x,u,t)$$
renders the set $\S = \S_x \cap \S_u$ forward invariant, i.e., safe, for appropriate initial conditions:
$$
\begin{array}{r}
x_0 \in \S_x  \\
(x_0,u_0) \in \S_e \cap \S_u
\end{array}
\quad \Rightarrow \quad
\begin{array}{r}
(x(t),u(t)) \in \S_x \cap \S_u \\
\forall ~ t \geq 0.
\end{array}
$$
\end{theorem}

\begin{proof}
The constraints in \eqref{eqn:CBFQPintfinal} are simply reformulations of:
\begin{eqnarray}
\hspace{-0.1cm}
\begin{array}{l}
 \dot{h}_e(x,u, v)  \geq  - \gamma_e(h_e(x,u))  \\
\dot{h}_u(x,u, v)  \geq   - \gamma_u(h_u(u))
\end{array}
\Leftrightarrow
\begin{array}{r}
 p_x(x)^{\top} v   \geq   d_e(x,u,t)  \\
p_u(x,u)^{\top} v  \geq  d_u(x,u,t)
\end{array}   \nonumber
\end{eqnarray}
Therefore, assuming a feasible QP implies that $\S_e \cap \S_u$ is forward invariant by Theorem \ref{thm:main} and Lemma \ref{lem:torquebound}.  The result follows by Lemma \ref{lem:cbftoicbf}.
\end{proof}

\begin{remark}
The requirement on the feasibility of the QP in Theorem \ref{thm:maintwo} encodes the fact that we are not assuming that $\S_x$ is a control invariant set.  Without this assumption, we do not know if there are feasible control inputs rendering $\S_x$ invariant---this is encoded by the feasibility of the QP. Moreover, if the QP is feasible, conditions for the Lipschitz continuity of its solution can be found (see, e.g., \cite{Ames17}). This and the relaxation of the feasibility assumption will be the subject of future work.
\end{remark}

\section{Simulation Results}

To demonstrate the results of the paper, we will return to an early motivating example considered for CBFs \cite{Ames14}: \emph{adaptive cruise control (ACC).}
The dynamics are given by:
\begin{eqnarray}
\label{eqn:accdynamics}
\dot{x} =
\begin{bmatrix}
x_2 \\
- \frac{1}{m} F_r(x) \\
v_0 - x_2
\end{bmatrix}
+
\begin{bmatrix}
0 \\
\frac{1}{m} \\
0
\end{bmatrix} u
\end{eqnarray}
where $(x_1,x_2)$ are the position and velocity ($x_2 = \dot{x}_1$) of the vehicle, $m$ its mass, $x_3$ is the distance between the vehicle and the lead vehicle traveling at a velocity of $v_0$, and $F_r(x) = c_0 + c_1 x_2 + c_2 x_2^2$ is the empirical form of rolling resistance. The control objective for the system is to drive the car to a desired speed ($x_2 \to v_d$).  This can be represented as an output: $y = \zeta(x) = x_2 - v_d$.  To obtain the predictor as in \eqref{eqn:predictor}, we can forward integrate \eqref{eqn:accdynamics} with $c_2 = 0$, i.e., the linear approximation, resulting in:
\[
\hat{y}(t + T) =
{\scriptstyle
-c_{1}^{-1}\left(c_{0}-u(t)+ m v_d-c_{1} e^{-\frac{c_{1} T }{m}} \left(x_2(t) +\frac{c_{0}-u(t)+m v_d}{c_{1}}\right)\right),
}
\]
yielding the dynamically defined control law per \eqref{eqn:controleqn}:
\[
\dot{u}(t) = \alpha c_1 \left( e^{\frac{-c_1}{m} T} - 1 \right)^{-1} \hat{y}(t + T)  =: \phi(x,u,t).
\]
The state safety constraints encode the ``half the speedometer rule'' which yields the CBF (since $p_x(x) = -1.8/m \neq 0$):
\[
h_x(x) = x_3 - 1.8 x_2 \geq 0 \qquad \mapsto \qquad
\S_x = \{h_x(x) \geq 0\} .
\]
The input constraint is given by the constraint that the wheel force is bounded by $|u| \leq m c_{a/d} g $ yielding an I-CBF:
\[
h_u(x) =  (m c_{a/d} g)^2 -  u^2 \geq 0
\quad \mapsto \quad
\S_u = \{h_u(u) \geq 0\} .
\]
\noindent where $c_{a/d}$ is the factor of $g$ for acceleration/deceleration.  In both cases, we pick $\gamma_x(r) = \gamma_u(r) = \gamma r$ for $\gamma > 0$.  From $h_x$ we get $h_e$ as in \eqref{eqn:heeqn}, and we pick $\gamma_e(r) = \frac{1}{2} \gamma r$.
\twocolumn[%
{\centering\includegraphics[width=\textwidth]{figures/final_fig4.pdf}
\captionof{figure}{Methods developed in the paper applied to adaptive cruise control: the goal is for a vehicle to reach a desired velocity ($v_d = 24 m/s$) while not colliding with the lead vehicle traveling slower (at $14 m/s$) and not exceeding a maximum wheel force (input bounds).}
\label{fig:finalfig}}
\vspace{0.5cm}
]
Utilizing the constructions leading to Theorem \ref{thm:maintwo} results in the forward invariance of $\S_x \cap \S_u$.  This is illustrated in Fig. \ref{fig:finalfig}, where the parameters and initial condition where chosen to match \cite{Ames14} with $\alpha = 10$ and $\gamma = 1$ above.
Three different controllers are plotted in Fig. \ref{fig:finalfig}, including: only input constraints and thus an I-CBF, $h_u$, with the controller from Lemma \ref{lem:torquebound} (red), only state constraints and thus a CBF, $h_x$, with the controller from Lemma \ref{lem:statebound} (green), and both input and state constraints, thus enforcing both $h_u$ and $h_e$, via the controller from Theorem \ref{thm:maintwo} (blue).  Therefore, via the presented unified controller, we are able to simultaneously satisfy the input and state constraints, and thus render the system safe for both constraints, i.e., render $\S_x \cap \S_u$ forward invariant.  Importantly, when compared against past uses of this example in \cite{Ames19,Ames14} we can do so holistically.

\section{Conclusions}

This paper introduced integral control barrier functions (I-CBFs).  By considering dynamically defined controllers, we were able to guarantee safety using I-CBFs defined in terms of both state and input.
This was applied in the context of dynamically defined tracking controllers for general nonlinear control systems (not necessarily control affine), wherein I-CBFs lead to minimal modification of these controllers via I-CBF based QPs (Theorem \ref{thm:main}).
We then considered the specific cases of input bounds (Lemma \ref{lem:torquebound}) and state constraints (Lemma \ref{lem:cbftoicbf}), wherein both can be independently enforced via I-CBFs.
Additionally, we proved the state constraints and input bounds can be simultaneously satisfied via the framework of I-CBFs assuming a feasible QP (Theorem \ref{thm:maintwo}).  This gives a holistic method for provably enforcing both through the use of safety-critical integral controllers.

Future work will be devoted to expanding this theoretic basis and demonstrating this theory experimentally.  From a theoretic perspective, understanding when the QP in Theorem \ref{thm:maintwo} is feasible is a rich problem that has important implications.
Additionally, we wish to improve the means in which nominal controllers can be enforced in conjunction with integral control barrier functions.
From an experimentally perspective, the ability of I-CBFs to bound inputs while achieving safety in state has important ramifications that we wish to explore on everything from walking robots to multi-robot systems to safey-critical autonomy.

\bibliographystyle{IEEEtran}