\section{Introduction}\label{sec:intro}

The purpose of this paper is to initiate the rigorous mathematical analysis of a model of the dynamics of disclination line defects in nematics proposed in \cite{ZANV}. Here, we focus on energetic aspects of the model. Combined with the ideas presented in \cite{AD} and the demonstrations provided in \cite{PAD, ZZAGW, ZANV}, which include static fields of straight $\pm \frac{1}{2}$  disclinations, their annihilation, the dissociation of closely bound pairs of straight disclinations, as well as static fields of disclination loops, the model can be considered as a thermodynamically consistent generalization of the Ericksen-Leslie (EL)  model to account for the dynamics of disclination lines, with total energy that remains bounded in finite bodies in the presence of these line defects.

The practical applications of liquid crystalline phases abound, from  liquid crystal displays for electronic devices and cell membranes in biology (use in the mechanically `soft' phase), to a vast variety of liquid crystal polymers including the mechanically `hard' Kevlar, for body armor, and in tires.  Equally, topological defects abound in liquid crystalline media, fundamentally due to microscopic structural symmetries related to the head-tail symmetry of the director. Due to this symmetry, a vector field assigned to a director field can undergo continuous changes in orientation around a non-unique surface terminating at a unique \textit{disclination line}, where the jump in orientation across the surface is quantized to be $\pi$ radians (see, e.g., De Gennes and Prost \cite[Sec.~4.2.1]{degennes-prost}). The line field representing the director and the corresponding vector field are both discontinuous at the disclination line, and the energy cost of such a line discontinuity can be sustained by the material. It is this fundamental insight, going back to the kinematics of singularities in linear elasticity theory due to Volterra and Weingarten (see, e.g., \cite{acharya2019weingarten} for a contemporary review), that forms the core idea of our model and, in fact, has been recently used to define an algorithm to detect line defects in molecular configurations of nematics produced by Molecular Dynamics simulations \cite{saha}, extending the notion of a director down to the level of a single nematic molecule.

The understanding of the energetics and \emph{dynamics} of topological defects and their interaction form an important part of the theoretical study of liquid crystalline media, and we are particularly interested in the universality of this behavior across the material behaviors of liquid crystals and crystalline solids. A primary justification of the type of model we consider is that it is inter-operable with a study of dislocation defects in crystalline solids with simply a change in interpretation of the field variables involved (this cannot be said of the Landau-DeGennes model \cite{Sonnet-Virga,schopohl1987defect,svenvsek2002hydrodynamics} which, nevertheless, employs a more adapted order parameter, the $Q$ tensor, to describe the head-tail symmetry of the nematic director, albeit at the cost of including a biaxial phase within the description as well). A distinguishing feature of our model is that, even in a `smooth' setting, defect cores can be identified as a locally calculable field variable, arguably a desirable feature as discussed in \cite{biscari2005field}.

Our model introduces an extra second-order tensor field, $B$, beyond the EL director field, $k$ (strictly speaking $k$ is a vector field representation of the director line field). This new field is to be physically thought of as a locally integrable realization, at the mesoscale, of the `singular' part of the director gradient field ($Dk$) in the presence of line defects, singular when viewed at the macroscale. Thus, at the mesoscale, both the director gradient field and the new field are integrable
\footnote{The length scales $\xi$ (core width) and $\e \xi$ (the layer width) that appear subsequently in Sec.~\ref{sec:main_res} and are at the heart of such a physical regularization, can be precisely defined in configurations of nematic molecules arising in Molecular Dynamics simulations, as shown in ongoing work \cite{saha}.}
- with this clear, we nevertheless refer to $B$ as the `singular part of the director distortion.' Notably, the field $B$ is not a gradient, and this allows it to encode information on the topological charge of line defects through its curl. The evolution of the director field $k$ continues to be obtained from the balance of angular momentum, as shown by Leslie \cite{L1992}, and the evolution of $B$ follows from a conceptually simple conservation law for the topological charge of the line defects, which is tautological before the introduction of constitutive assumption for the disclination velocity, the latter deduced from consistency with the second law of thermodynamics. The introduction of dynamics based on such a conservation law, rooted in the kinematics of defect lines, is a conceptual departure from what is done for dynamics with the Landau-DeGennes $Q$ tensor model (see, \cite{Sonnet-Virga, macmillan_1}), or Ericksen's variable degree of orientation model \cite{Ericksen_1991}. In doing so, the model also makes connections to the dynamics of dislocation line defects in elastic solids \cite{ZAWB,AZA}, as well as their statics \cite{Arora_AA_2020}.

At the length scales where individual line defects are resolved, partial differential equations-based dynamical models arising from continuum mechanical considerations involve Newtonian and thermodynamic driving forces that include nonlinear combinations of entities representing director distortions and the disclination density fields. This requires a minimum amount of regularity in these fields, and hence it is essential to have a formulation that utilizes at least locally integrable functions, and our model is designed to be consistent with this requirement (of course, this does not preclude the question of studying limiting situations of such models when such functions tend to singular limits, modeling fields that have discontinuities, and singularities in the macroscopic limit).

\subsection{Main Results}\label{sec:main_res}
We investigate the behavior of local minimizers for the previously discussed model for liquid crystals with disclination defects. Let $\Omega\subset \mathbb{R}^N$, where $N=2$ or $3$, be the domain occupied by a nematic liquid crystal. We consider the following energy for the director field $k \in W^{1,2}(\o;\mathbb{R}^2)$ and the singular part of the director distortion $B \in H_\text{curl}(\o;\mathbb{R}^{2\times 2})$,
 \begin{equation}\label{eq:energy}
 	E_{\e,\xi}[k,B] := \int_\o{\left[\frac{(|k|-1)^2}{\e\xi^2} + |\nabla k - B|^2 + \e \xi^2|\mathrm{curl}\, B|^2 + \frac{1}{\e \xi^2}W(\e \xi|B|)\right]dx},
 \end{equation}
where $W: [0,\infty) \to [0,\infty)$ is a nonconvex double-well continuous potential with wells at 0 and 2.

\begin{remark} \label{rmk:heuristics}
The heuristics behind the energy \eqref{eq:energy} above for the prediction $\pm \frac{1}{2}$ line defects are as follows: the nonconvex potential $W$ assigns vanishing energy cost when $|B| = 0$ or $ \frac{2}{\e \xi}$. This, along with the elastic energy term $| D k - B|^2$, assigns approximately vanishing elastic cost for pointwise values of the director gradient of the type $Dk \approx \frac{n - (-n)}{\e \xi} \otimes l$, where $0 < \e \leq 1$ and $n, l$ are unit vectors. Here, $n$ corresponds to the director field $k$, and $l$ represents the direction along which the jump of $n$ occurs.

In Fig. \ref{fig:gcr}, this scenario is presented for a rectangular transition layer for the director field.

 \begin{figure}
	\centering
	\includegraphics[width=0.5\linewidth]{GCR}
	\caption{A representation of the domain and the layer where the discontinuity is supported}
	\label{fig:gcr}
\end{figure}

If the transition layer in Fig. \ref{fig:gcr} did not terminate inside the domain, the energy cost would be minimal for a (diffuse) jump in the orientation of $k$ by $\pi$ radians across the layer. However, for a layer terminating inside the domain, $\mbox{curl} \, B$ is non-zero near the termination (or core), and if the width of the layer in the vertical direction was $0$ (i.e., $\e = 0$), $\mbox{curl} \, B$ would be singular (the classical defect solution results from the choice $\xi = 0$, when layer and core widths vanish). When the $\mbox{curl} \,$ does not vanish, the density $Dk$ cannot annihilate $B$ (regardless of $\xi = 0$, or not). To see this, the Euler Lagrange equation of a functional with just the energy density $|Dk - B|^2$ for admissible variations in $k$ with $B$ a specified field is, with $Dk - B  =: e$, $\mbox{div} \, e = 0$, and $\mbox{curl} \, e = -\mbox{curl} \, B$.

As an example, for $\mbox{curl} \, B$ a (mollified) Dirac supported at the layer termination, this produces the approximate elastic energy density field, given here by $|e|^2$, of some canonical line defects in 2 dimensions (screw dislocation in solids, the wedge disclination in nematics with unit vector constraint imposed, either exactly or approximately) \cite{ZZAGW,ZANV, nabarro, hirth_lothe, frank1958liquid, degennes-prost}. Since $e = Dk$ outside the layer, we recover the relevant director field (using the penalized unit vector constraint represented by the first term in \eqref{eq:energy} and a specified value of $k$ at one point of the domain). Within the layer, but outside core, the director field $k$ flips orientation by $\pi$ radians, with a somewhat more involved distribution in the core.

\end{remark}

 As we are interested in minimizers of the energy (\ref{eq:energy}), we first consider the relaxation at a fixed scale (see Theorem \ref{thm:relax} for complete details). The energetic relaxation provides a functional to which the direct method of the calculus of variations is amenable, and is the first step to understanding the structure of minimizers.

 \begin{thm*}
Let $\Omega\subset \mathbb{R}^N$, $N=2$ or $3$, be an open, bounded set with Lipschitz boundary. For $\e,\xi>0$ fixed, the lower semicontinuous relaxation of the energy (\ref{eq:energy}) is given by
	\begin{align}
		\bar E_{\e,\xi}[k,B] &:=
			\iO\left[\frac{(|k|-1)^2}{\e\xi^2}+|\nabla k - B|^2\right]\;dx \nonumber\\
			&\quad\quad + \iO\left[\e\xi^2|\mathrm{curl}\, B|^2 + \frac{1}{\e\xi^2}Q(W(|\cdot|))(\e\xi B)\right] \;dx.
	\end{align}
 for all $k \in W^{1,2}(\o;\mathbb{R}^N)$ and $B \in H_\mathrm{curl}\,(\o;\R^{N\times N})$.
	Here, $Qf$ denotes the quasiconvex envelope of $f$.
\end{thm*}

\begin{remark}	\label{rmk:quasiconvex}
	We show that $Q(W(|\cdot|))\equiv 0 $ in the ball $B(0,2)\subset \mathbb{R}^{N\times N}$. Clearly, $Q(W(|\cdot|))(p) \geq W^{**}(|p|)$ always. This shows that $Q(W(|\cdot|))(p) \geq 0$ everywhere. Further, we claim that each $p \in B(0,2)\subset \mathbb{R}^{N\times N}$ is given as the convex combination of two-elements in $2S^{N\times N-1}$ differing by a rank-one matrix. Indeed, note we can find $\lambda_+,\lambda_- \geq 0$ such that $$p_\pm:=p \pm \lambda_\pm e_1\otimes e_1 \in 2S^{N\times N-1}\quad\text{and}\quad \mathrm{rank}(p_+-p_-)=1.$$
		By taking $t= \frac{\lambda_-}{\lambda_++\lambda_-} \in [0,1]$, we can compute that
	$$t(p + \lambda_+ e_1\otimes e_1) + (1-t)(p - \lambda_- e_1\otimes e_1) = p.$$
	As $W(|\cdot|)\llcorner_{2S^{N\times N-1}} \equiv 0$ by (\ref{hyp:f2}), we apply rank-one convexity of the quasiconvex envelope \cite{Dacoragna} and the fact that the envelope is always below the original function to find that
 \begin{align*}
 Q(W(|\cdot|)) (p) \leq & tQ(W(|\cdot|)) (p_+) + (1-t)Q(W(|\cdot|)) (p_-) \\
 \leq & tW(|p_+|) +(1-t)W(|p_-|) = 0,
 \end{align*} thereby concluding that $Q(W(|\cdot|))$ vanishes in the ball $B(0,2)$.

 We further conjecture that $Q(W(|\cdot|))(p) = W^{**}(p)$ due to the radial symmetry, however characterization of the quasiconvex envelope poses challenges even in simple cases (see, e.g., \cite{QCenvelope} for one of the few nontrivial examples of a calculation of the envelope).
\end{remark}

To motivate the constraints we will place on the field $B$, we introduce a simple example. We now restrict our attention to dimension $N=2$ and consider the limit $\e \to 0$ with $\xi > 0$ fixed. Considering any $k \in C^{1}(\Omega;S^{1})$, we set $B := \nabla k$ to find that $\bar E_{\e,\xi}[k,\nabla k] = \int_\Omega\frac{1}{\e\xi^2}Q(W(|\cdot|))(\e\xi \nabla k) \,dx.$ Given Remark \ref{rmk:quasiconvex} above, the function $Q(W(|\cdot|))(\e\xi B)$ doesn't see the energy from $B$ if $|B|\leq \frac{2}{\e \xi}$. As a result, if $\epsilon \leq \frac{2}{\xi\|\nabla k\|_\infty +1}$, then $\bar E_{\e,\xi}[k,\nabla k]=0$. Such a result (though defect free) shows that further constraints on the field $B$ are required to gain physically meaningful insight in the limit as $\e\to 0$.

We consider the particular case of $B$ supported in a layer as in Figure \ref{fig:gcr}, a physically relevant geometric configuration (see, e.g., \cite{ZANV}).

To be precise, let $\Omega:= (-1,1)^2$ be the domain of a liquid crystal in the plane. We assume the defect is at the origin, and the surface of discontinuity is within a layer $L_{\varepsilon,\xi}$, with
\begin{equation}\label{eqn:BlayerConstraint}
{L_{\e,\xi}}:= (-\xi,1)\times\left(-\frac{\e\xi}{2},\frac{\e\xi}{2}\right) \quad \text{ and } \quad B \equiv 0 \text{ in }\Omega \setminus L_{\e,\xi}
\end{equation} with parameters $\e,\xi > 0$. In physical terms, $\xi$ is the core length of the crystalline defect and $\e$ is a parameter determining the thickness of the defect layer $L_{\e,\xi}$ (see Figure \ref{fig:gcr}).  In this paper, we are primarily concerned with $\pm\frac{1}{2}$ disclinations, which must satisfy the constraint
  \begin{equation}\label{eq:halfdisclin}
 \left|\int_\o{\mathrm{curl}\, B\;dx} \right| = 2.
 \end{equation}
 This is a model constraint requiring a disclination to exist in the domain. By Stokes' theorem, (\ref{eq:halfdisclin}) is consistent with a layer field in the form of $ B = \frac{n - (-n)}{\e \xi} \otimes l$ in a layer of width $\e \xi$ with normal in the direction $l$ and $n$ a unit vector corresponding to the director field $k$ (see Fig. \ref{fig:gcr}), and as described in Remark \ref{rmk:heuristics}.

After a change of variables analogous to typical dimension reduction problems \cite{DretRaoult}, we prove a compactness theorem for the rescaled fields, which is precisely stated in Theorem \ref{thm:compactness}. Next, we state a theorem that follows from Theorem \ref{thm:compactness} and which emphasizes the coupling of the physical quantities in the asymptotic limit.

\begin{theorem}\label{thm:compactCorollary}
    Let $k_\e$ and $B_\e$ have uniformly bounded energy as $\e \to 0$, that is, $\sup_{\e>0}E_{\e,\xi}[k_{\e},B_{\e}]<\infty$. Further, suppose $B_{\e}$ satisfies the geometric constraint (\ref{eqn:BlayerConstraint}) and corresponds to $\pm 1/2$ disclination by satisfying (\ref{eq:halfdisclin}). Then up to a subsequence (not relabeled), $k_\e \to k$ strongly in $L^2(\o;\mathbb{R}^2)$ where $k \in SBV(\o;S^1)\cap W^{1,1}(\Omega\setminus L_{0,\xi};S^1)$. Defining the jump of $k$ on $(-\xi,1)$ by $\jump{k}: = k^+ -k^-$, the compatibility condition
	$$\jump{k}(s) = \int_{-\xi}^{s}\int_{-\frac{\xi}{2}}^{\frac{\xi}{2}} \alpha \;dx \quad \text{ and }\quad |\jump{k}(1)| = 2 $$
	is satisfied for all $s\in (-\xi,1)$, where $\alpha\in L^2(L_{1,\xi})$ is the limit of the rescaled $\{\tilde B_\e := \e B(x_1,\e x_2)\}_{\e}.$ Furthermore, $\nabla k$, the part of $Dk$ absolutely continuous with respect to the Lebesgue measure has higher regularity, in the sense that
    $$\mathrm{curl}\,(\nabla k) = - \frac{d}{dx_1}\jump{k}\otimes e_2\mathcal{H}^{1}\llcorner_{(-\xi,1)\times\{0\}}.$$
\end{theorem}

We can also make a connection to the recent preprint \cite{GoldmanSerfaty}, where a SBV model for $\pm \frac{1}{2}$ disclinations is proposed and the constraint that $[k] = 2$ along the jump set is imposed. The energy we use can be viewed as an attempt to also relax the one used in \cite{GoldmanSerfaty} by being a Sobolev model allowing for a more general class of jumps in the SBV limit.

There are many open questions stemming from this work which we highlight in Section \ref{sec:conclusion}. Foremost, is the integral representation of a precise limiting energy for the case $\e \to 0$ with $\xi > 0$ fixed. It is also possible to consider the case of $\xi \to 0$ at various rates compared to $\e \to 0$. However, the limit $\xi \to 0$ will be complicated by the need to rescale the energy by $\log \xi$, which leads to a delicate Ginzbug-Landau type problem (see \cite{JerrardSoner}, \cite{AlicandroPonsig}).

\section{Mathematical Preliminaries}

Let $W^{1,2}(\o;\mathbb{R}^2)$ denote the usual Sobolev space, and we designate by $H_\text{curl}(\o;\mathbb{R}^{2\times 2})$ the space of $L^2$ matrix valued tensors, whose row-wise distributional curl is also in $L^2$. Under this setting, we consider the energy (\ref{eq:energy}) with a nonconvex continuous potential
$W: [0,\infty) \to [0,\infty)$ satisfying the following coercivity and growth properties,
\begin{equation}\label{hyp:f1}
	\frac{1}{C}|x|^2-C\leq W(x) \leq C(1+|x|^2) \quad \text{for }\; x \in [0,\infty)
\end{equation}
for some $C>0$, and
\begin{equation}\label{hyp:f2}
\{x \in [0,+\infty): W(x) = 0\} = \{0,2\}.
\end{equation}

The mathematical framework we use to study the convergence of the functional \eqref{eq:energy} is encapsulated by the notion of $\Gamma$-convergence, which we recall next.

\begin{definition}\label{def:Gamma_conv}
Given a metric space $(X,d)$, let $F_n: X \to [0,\infty]$ be a sequence of functionals for $n\in \mathbb{N}$. We say that $F_n$ $\Gamma$-converge to $F_0: X \to [0,\infty]$ with respect to the metric $d$ if the following two conditions hold:
\begin{enumerate}
    \item (Liminf Inequality) For every $u \in X $ and for every sequence $\{u_n\} \subset X $ such that $u_n \to u$ with respect to the metric $d$, we have $$F_0(u) \leq \underset{n \to \infty}{\liminf}\;{}{F_n(u_n)}.$$
    \item (Recovery Sequence) For every $u \in X$, there exists $\{u_n\} \subset X$ such that $u_n \to u$ with respect to the metric $d$, and the sequence recovers the energy, i.e., $$\underset{n \to \infty}{\limsup}\;{F_n(u_n)}= F_0(u).$$
\end{enumerate}
\end{definition}

For the relaxation of the energy (\ref{eq:energy}), which is by definition the $\Gamma$-limit of the constant sequence of functionals $F_n :=E_{\e,\xi}$, we will rely on the now classical notion of quasiconvexity introduced by Morrey \cite{Morrey}. Analogous to the characterization of convex functions via Jensen's inequality, quasiconvex functions satisfy a Jensen's type inequality for gradient fields. Specifically, a Borel measurable function $g: \mathbb{R}^{d\times N}\to [0,\infty]$, $d,N\geq 1$, is quasiconvex if
$$g(\xi) \leq \int_{(0,1)^N} g(\xi + \nabla \phi) \, dx \quad \quad \text{ for all }\xi \in \mathbb{R}^{d\times N} \text{ and }\phi \in C^1_c((0,1)^N;\mathbb{R}^d).$$
The quasiconvex envelope of a function $f$ is given by the greatest quasiconvex function beneath $f$, i.e., it is defined pointwise by
\begin{equation}\label{def:quasiconEnv}
Qf(x_0) : = \sup\{g(x_0):g \text{ is quasiconvex and }g\leq f\}.
\end{equation}
We refer the reader to \cite{Dacoragna} for further details on such functions.

As mentioned in the discussion preceding Theorem \ref{thm:compactCorollary}, we wish to model discontinuities across 2-d surfaces - when viewed at the macroscale -  in a vector field representation of the director field containing a line defect, while incorporating the fact that at the microscopic scale such a jump across the surface must necessarily be spread over a region roughly of the order of the spacing between adjacent mesogens (cf. \cite{saha}). Specifically, we impose the condition that $B$ vanishes outside of the layer, which is equivalent to the conditions
   \begin{align}
  	&B = 0 \quad\quad \mathrm{ in }\quad \o \setminus L_{\e,\xi}, \label{eq:B_zero}\\
  	&Bt = 0\quad\quad \mathrm{ on }\quad \partial {L_{\e,\xi}} \setminus \partial \o, \label{eq:B_tan}
  \end{align}
where $t$ is the tangent vector to the boundary point. The condition \eqref{eq:B_tan} comes for free as $B \in H_{\mathrm{curl}\,}(\o;\mathbb{R}^{2\times 2})$ and there is a well-defined tangential trace matching the condition \eqref{eq:B_zero} \cite{FabrieBoyer}.

Finally, we recall that a function $u$ belongs to $BV(\Omega)$ if its distributional gradient is given by a finite Radon measure. Informally, in the case that $u$ has only surface discontinuities, $u$ belongs to $SBV(\Omega)$, and for the sake Theorem \ref{thm:compactCorollary}, it suffices to know that if $u \in W^{1,1}(\Omega \setminus K)\cap L^\infty(\Omega)$, where $K$ is a closed set with finite surface measure, i.e., $\mathcal{H}^{N-1}(K)<\infty,$ then $u \in SBV(\Omega)$. For further details, we refer to \cite[Proposition 4.4]{AmbrosioFuscoPallara}.

\section{Relaxation for fixed $\e,\xi$}
We study properties of the energy (\ref{eq:energy}) with $\e,\xi>0$ fixed. A priori, it is not clear that minimizers to the problem exist nor is it clear what the value of the infimum is. In order to apply the direct method of the calculus of variations, we must consider the lower semicontinuous envelope of the functional. Specifically, we obtain an integral representation for the relaxation of energy (\ref{eq:energy}) in dimensions $N=2$ or $3$. This dimension constraint enables us to use the Helmholtz decomposition and the corresponding Sobolev spaces, as detailed in \cite{FabrieBoyer}. Here, $Qf$ denotes the quasiconvex envelope of $f$ as in (\ref{def:quasiconEnv}).
\begin{theorem}\label{thm:relax}
	Let $\Omega\subset \mathbb{R}^N$, with $N=2$ or $3$, be an open, bounded set with Lipschitz continuous boundary, and let $E$ be defined in \eqref{eq:energy}. For all $k \in W^{1,2}(\o;\mathbb{R}^N)$ and $B \in H_\mathrm{curl}\,(\o;\R^{N\times N})$, the relaxation of $E_{\e,\xi}$ is given by
	\begin{equation}\nonumber
		\bar E_{\e,\xi} [k,B] : = \inf \left\{\underset{n \to \infty}{\liminf}\; E_{\e,\xi}[k_n, B_n]: (k_n,B_n) \to (k,B) \right\},
	\end{equation}
	where  the convergence is such that $k_n \to k$ and $B_n \wk{} B$ in $L^2$.
	This relaxed energy has the integral representation
	\begin{align}\label{def:Ebar}
		\bar E_{\e,\xi}[k,B] &:=
			\iO\left[\frac{(|k|-1)^2}{\e\xi^2}+|\nabla k - B|^2\right]\;dx \nonumber\\
			&\quad\quad + \iO\left[\e\xi^2|\mathrm{curl}\, B|^2 + \frac{1}{\e\xi^2}Q(W(|\cdot|))(\e\xi B)\right] \;dx.
	\end{align}

\end{theorem}

We note that this result can equivalently be phrased as $E_{\e,\xi}$ $\Gamma$-converge to $\bar E_{\e,\xi}.$
To prove Theorem \ref{thm:relax}, we will introduce an intermediate functional, related to (\ref{def:Ebar}) through the Helmholtz decomposition of $B$.
We denote the space of the divergence-free fields with integrable $\mathrm{curl}\,$ as
\begin{equation}\label{def:Cspace}
	\mathcal{C} : = \{u \in H_{\mathrm{curl}\,}(\Omega;\R^{N\times N}) \; : \ {\rm div }\,  u = 0, \ u \cdot \nu = 0 \ \text{on } \partial \Omega \}.
\end{equation}
Define the functional $I: (L^2(\Omega;\mathbb{R}^N))^2 \times L^2(\Omega;\R^{N\times N}) \to [0,+\infty]$ by
\begin{equation}\label{def:I}
\begin{aligned}
    	I[\tilde{k},z,p] := &\int_\Omega\left[(|\tilde{k}+z|-1)^2+|\nabla \tilde{k}|^2\right]\;dx \\
		&\quad +\iO\left[|p|^2+|\mathrm{curl}\, p|^2 + W(|\nabla z + p|)\right]\;  dx,
\end{aligned}
\end{equation}
if $(\tilde{k},z,p) \in \mathcal{X}$, and $+\infty$ otherwise.
Here,
\begin{equation}\label{def:Xspace}
	\mathcal{X}:= W^{1,2}(\Omega;\mathbb{R}^N)\times \left(W^{1,2}(\Omega;\mathbb{R}^N) \cap \left\{\int_\Omega z \, dx = 0\right\}\right) \times \mathcal{C}.
\end{equation}

First, we investigate compactness of the functional $I$ in the following lemma.
\begin{lemma}[Compactness of $I$]\label{lem:compact}
	Let $\Omega \subset \mathbb{R}^N$, with $N=2$ or $3$, be an open, bounded set with Lipschitz continuous boundary. Consider a sequence $\{(\tilde{k}_n,z_n,p_n)\}$ such that $\sup_{n \in \mathbb{N} }{I[\tilde{k}_n,z_n,p_n]}\leq C$.
	Then there is $(\tilde{k},z,p) \in \mathcal{X} $ such that, up to a subsequence (not relabeled),
	\begin{align*}
		\tilde{k}_n &\wk{} \tilde{k} \quad \text{in } W^{1,2}(\Omega;\mathbb{R}^N),\\
		z_n &\wk{} z \quad \text{in } W^{1,2}(\Omega;\mathbb{R}^N),\\
		p_n &\wk{} p \quad \text{in } W^{1,2}(\Omega;\R^{N\times N}),
	\end{align*}
	and
	$$(\tilde{k}_n,z_n,p_n) \to (\tilde{k},z,p) \quad \text{strongly in } (L^2(\Omega;\mathbb{R}^N))^2 \times L^2(\Omega;\R^{N\times N}).$$
\end{lemma}
\begin{proof}
	The key is the inequality (see \cite{FabrieBoyer})
	\begin{equation}\label{bdd:div+curl}
		\|p_n\|_{W^{1,2}(\Omega; \mathbb{R}^N)} \leq C(\Omega)\left(\|p_n\|_{L^2(\Omega;\mathbb{R}^N)} + \|{\rm div }\, p_n \|_{L^2(\Omega)} + \|\mathrm{curl}\, p_n \|_{L^2(\Omega;\mathbb{R}^N)} \right).
	\end{equation}
	From the definition of $\mathcal{C}$ in (\ref{def:Cspace}), we have that ${\rm div }\,  p_n = 0$. As $I$ in (\ref{def:I}) controls $\|\mathrm{curl}\, p_n\|^2_{L^2(\Omega)}$, we apply the above inequality to conclude that $\|p_n\|_{W^{1,2}(\Omega,\R^{N\times N})}\leq C<\infty$. Convergence as in the statement of the lemma follows from weak compactness. By (\ref{hyp:f1}) and (\ref{def:I}), we have $\sup_n\|\nabla z_n\|_{L^2(\Omega)} \leq C<\infty,$ and the desired convergence follows from Poincar\'e's inequality because $\int_\Omega z_n \, dx = 0.$ Finally, the uniform bound of the energy (\ref{def:I}) implies a uniform bound on $\|\nabla \tilde k_n\|_{L^2(\Omega)}$. Combining this with control of $\|z_n\|_{L^2(\Omega;\mathbb{R}^N)}$ and $\|\tilde k_n\|_{L^2(\Omega;\mathbb{R}^N)}$, shows that $\sup_n\|\tilde k_n\|_{W^{1,2}(\Omega;\mathbb{R}^N)}\leq C<\infty$.

	To conclude strong convergence in $L^2,$ we apply the Rellich-Kondrachov compactness theorem.
\end{proof}

We will prove the relaxation of the functionals $I$ using techniques for \\
$\Gamma-$convergence (see Definition \ref{def:Gamma_conv}), i.e., that $I$ $\Gamma$-converges to $\bar{I}$, with
\begin{equation}\label{def:I_limit}
\begin{aligned}
    	\bar I[\tilde{k},z,p] :=& \int_\Omega \left[(|\tilde{k}+z|-1)^2+|\nabla \tilde{k}|^2\right]\;dx \\
		& \quad +\iO\left[|p|^2+|\mathrm{curl}\, p|^2 + Q(W(|\cdot|))( \nabla z + p)\right]\;  dx,
\end{aligned}
\end{equation}
if $(\tilde{k},z,p) \in \mathcal{X}$, and $+\infty$ otherwise.

We note that this functional is the same as the original functional with $W$ replaced by $Q(W(|\cdot|))$. In order to prove that this is indeed the correct limiting energy, we first show that the $\liminf$ inequality of $\Gamma-$convergence is satisfied (see Definition \ref{def:Gamma_conv}(1)).

 \begin{lemma}[Liminf of $I$]\label{lem:liminfI}
 	Let $\Omega\subset \mathbb{R}^N$, with $N=2$ or $3$, be an open, bounded set with Lipschitz continuous boundary, and assume that $W$ satisfies (\ref{hyp:f1}) and (\ref{hyp:f2}). For all sequences such that $$(\tilde{k}_n,z_n,p_n)\to (\tilde{k},z,p) \quad \text{strongly in } (L^2(\Omega;\mathbb{R}^N))^2 \times L^2(\Omega;\R^{N\times N}),$$ we have
 	$$
 	\bar{I}[\tilde{k},z,p]\leq \underset{n \to \infty}{\liminf}\; I[\tilde{k}_n,z_n,p_n].
 	$$

 \end{lemma}

 \begin{proof}
 	We define the function $h(s,\eta) : = W(|s + \eta|)$ and $Q h(s,\cdot)$ to be the greatest quasiconvex function below $h(s,\cdot)$ as in (\ref{def:quasiconEnv}).

    We claim that
 	\begin{multline}\label{bdd:hLSC}
 		\int_\Omega Q(W(|\cdot|))(p+\nabla z ) \, dx = \int_\Omega Q h(p,\nabla z) \, dx \\
 		\leq \underset{n \to \infty}{\liminf}\; \int_\Omega h(p_n,\nabla z_n) \, dx = \underset{n \to \infty}{\liminf}\; \int_\Omega W(|p_n +\nabla z_n|)\, dx.
 	\end{multline}
 	The first equality is easy to see because we note that there is a translational symmetry to $h$ which gives the equality
 	\begin{equation}\label{eqn:quasiTranslation}
 	Q h(s,\nabla z) = Q h(0,\nabla (sx+z(x))) = Q(W(|\cdot|))(s+\nabla z )
 	\end{equation}
 	for all $s \in \mathbb{R}$ and almost every $x \in \mathbb{R}^N$; the first inequality in (\ref{bdd:hLSC}) follows choosing $s = p(x)$ pointwise. To obtain the second inequality, note that $Qh(s,\eta)$ is $2$-Lipschitz in the second variable, so by (\ref{eqn:quasiTranslation}), it is continuous in the $s$-variable and hence $Qh$ is Carath\'eodory. Applying the lower semicontinuity result of Acerbi and Fusco for quasiconvex functions \cite{AcerbiFusco}, we find that
 	\begin{equation}\nonumber
 	\int_\Omega Q h(p,\nabla z) \, dx \leq \underset{n \to \infty}{\liminf}\; \int_\Omega Qh(p_n,\nabla z_n) \, dx	\leq \underset{n \to \infty}{\liminf}\; \int_\Omega h(p_n,\nabla z_n) \, dx,
 	\end{equation}
 	concluding (\ref{bdd:hLSC}).

 The lower semi-continuity of the $L^2$ norm shows that
 	\begin{equation}\label{bdd:LSCL2}
 		\underset{n \to \infty}{\liminf}\; \int_\Omega \left(|\nabla \tilde k_n|^2+|\mathrm{curl}\, p_n|^2 \right) \, dx \geq \int_\Omega \left(|\nabla \tilde k|^2 + |\mathrm{curl}\, p|^2\right) \, dx,
 	\end{equation}
 	with strong convergence of $z_n$, $k_n$, and $p_n$ in $L^2$ giving that
 	\begin{equation}\label{eqn:strConv}
 		\lim_{n \to \infty} \int_\Omega \left[(|\tilde{k}_n + z_n|-1)^2  + |p_n|^2\right] dx = \int_\Omega \left[(|\tilde{k}+z|-1)^2+ |p|^2 \right] dx.
 	\end{equation}
 	Combining (\ref{bdd:hLSC}), (\ref{bdd:LSCL2}), and (\ref{eqn:strConv}), the lemma is proven.
 \end{proof}

In order to complete the integral representation for the relaxation of $I$, we show the existence of a recovery sequence (see Def \ref{def:Gamma_conv}(2)).

\begin{lemma}[Recovery Sequence for $I$]\label{lem:limsupI}
	If $(\tilde{k},z,p) \in \mathcal{X}$, then there exists a sequence $\{(\tilde{k}_n,z_n,p_n)\}$ such that
	\begin{align}
		(\tilde{k}_n,z_n,p_n)\to (\tilde{k},z,p) \quad &\text{strongly in } (L^2(\Omega;\mathbb{R}^N))^2 \times L^2(\Omega;\R^{N\times N})\label{cv}\\
		\underset{n \to \infty}{\limsup}\; I[\tilde{k}_n,&z_n,p_n]\leq \bar{I}[\tilde{k},z,p]. \label{ls}
	\end{align}

\end{lemma}

\begin{proof}
    As $\tilde{k}$ and $p$ are admissible in the original energy, we take $p_n := p$ and $\tilde{k}_n := \tilde{k}$.
	Now define the function $g:\Omega \times \R^{N\times N} \to [0,\infty)$ by $g(x,\eta) = W(|p(x)+\eta|)$. Note that as $W$ is continuous and $p$ is measurable, $g$ is a Carath\'eodory function and has polynomial growth in $\eta$.
	By standard relaxation results \cite{Dacoragna}, we can find a sequence $ \{z_n\}\subset W^{1,2}(\Omega;\mathbb{R}^N)$ such that
	\begin{align}
		&z_n \to z \quad \text{strongly in } L^2(\Omega;\mathbb{R}^N), \label{str}\\
		&\underset{n \to \infty}{\limsup}\; \iO{g(x,\gr z_n(x))\; dx}\leq \iO{Qg(x,\gr z(x))\; dx}\label{qc}
	\end{align}
    where $Qg(x,\eta)$ is defined in (\ref{def:quasiconEnv}). Again, by a similar argument for (\ref{eqn:quasiTranslation}), we have that $$Qg(x,\gr z(x)) =  Q(W(|\cdot|))(p(x)+\nabla z(x)) $$ for $x\in \Omega$ almost everywhere.
	By the above relation, \eqref{str}, and \eqref{qc}, we obtain (\ref{ls}) as the other functions in $I$ (\ref{def:I}) are fixed.
\end{proof}

Combining the last two lemmas, we finish the proof of Theorem \ref{thm:relax}.

\begin{proof}[Proof of Theorem \ref{thm:relax}]
As $\e$ and $\xi$ are fixed, we will let $\e = \xi = 1$ without loss of generality.
Note that for any $B \in L^2(\Omega; \mathbb{R}^{N \times N})$ we can apply the Helmholtz decomposition (see \cite{FabrieBoyer}) row-wise to find $B = p + \nabla z$ for some $p \in \mathcal{C}$ that is divergence free (see (\ref{def:Cspace})) and $z \in W^{1,2}(\Omega;\mathbb{R}^N).$ As the $\mathrm{curl}\,$ vanishes on gradients and $p\in \mathcal{C}$ and $\nabla (z-k)$ are orthogonal in $L^2(\Omega;\R^{N\times N})$ by an integration by parts as ${\rm div }\,   p = 0$ and $p\cdot \nu = 0$ (see (\ref{def:Cspace})), we have
\begin{align*}
E_{1,1}[k,B] &= E_{1,1}[k,p+\nabla z]  \\
& = \int_\Omega \left[(|k|-1)^2 + |\nabla k - (p+\nabla z)|^2 + |\mathrm{curl}\, p |^2 + W(p + \nabla z)\right] dx\\
&= I[k-z,z,p],
\end{align*}
where we have used the fact that
\begin{align*} \int_\Omega |\nabla k - (p+\nabla z)|^2 \ , dx &= \int_\Omega \left[|\nabla k - \nabla z|^2 + 2 p \cdot\nabla( z-  k)   + |p|^2 \right] dx \\
&= \int_\Omega \left[|\nabla k - \nabla z|^2  + |p|^2\right] dx,
\end{align*}
and likewise
$$\bar E_{1,1}[k,B] = \bar I[k-z,z,p].$$
Using the above relations, one has that Lemma \ref{lem:liminfI} and Lemma \ref{lem:limsupI} translate to $\liminf$ and $\limsup$ relations for $E_{1,1}$, thereby proving $E_{1,1}$ $\Gamma$-converges to $\bar E_{1,1}$ and concluding the proof.
\end{proof}

\section{Constrained Minimizers}
In the following, we study low energy sequences for the energy $E_{\e,\xi}$ in dimension $N=2$ as $\e \to 0$, with $\xi>0$ fixed, when the fields $B$ satisfy geometric (\ref{eqn:BlayerConstraint}) and defect (\ref{eq:halfdisclin}) constraints.

Here the disclination layer ${L_{\e,\xi}}$ becomes thin in the limit. As it is typical in dimension reduction problems, we perform the change of variables
\begin{align}\label{eq:Change_of_Var}
	\tilde{k}(x_1,x_2):=k\left(x_1,\e x_2\right),\\
	\tilde{B}(x_1,x_2):= \e B\left(x_1,\e x_2\right).
\end{align}
We remark that we have rescaled the field $B$ as well because the quadratic coercivity of $W$ only gives control over $\e B$.
With $L_\xi:= (-\xi,1)\times(-\frac{\xi}{2},\frac{\xi}{2})$, we write the energy in terms of a bulk contribution and layer contribution as

 \begin{align}
	&E_{\e,\xi}[k,
	B] = \int_{\o\setminus{L_{\e,\xi}}}{\left[\frac{(|k|-1)^2}{\e\xi^2} + |\nabla k|^2 \right]\;dx}  \nonumber\\
	& + \e \int_{L_\xi} {\left[\frac{(|\tilde{k}|-1)^2}{\e\xi^2} + \left|\nabla_{\e} \tilde{k} - \frac{\tilde{B}}{\e}\right|^2+\frac{\xi^2}{\e}|\mathrm{curl}_{\e}\,\tilde{B}|^2 + \frac{1}{\e \xi^2}W(\xi|\tilde{B}|)\right]\; dx} \nonumber\\
	&= \int_{\o\setminus{L_{\e,\xi}}}{\left[\frac{(|k|-1)^2}{\e\xi^2} + |\nabla k|^2 \right]\;dx} \nonumber\\
	&\quad + \int_{L_{\xi}} {\left[\frac{(|\tilde{k}|-1)^2}{\xi^2} + \e\left|\nabla_{\e} \tilde{k} - \frac{\tilde{B}}{\e}\right|^2+\xi^2|\mathrm{curl}_{\e}\, \tilde{B}|^2 + \frac{1}{\xi^2}W(\xi|\tilde{B}|)\right]\; dx} \nonumber\\
	 & =: E^{bulk}_{\e,\xi}[k] + E^{layer}_{\e,\xi}[\tilde{k},\tilde{B}] \label{eq:bulk+layer}
\end{align}
where $\nabla_{\e} := [\partial_1,\frac{1}{\e}\partial_2]$ and the scaled curl operator is $\mathrm{curl}_{\e}\, g := \partial_1 g_2 - \frac{1}{\e}\partial_2 g_1$.
Furthermore, the $\mathrm{curl}\,$ constraint \eqref{eq:halfdisclin} becomes
 \begin{equation}\label{eq:curl_contraint_cov}
	\left|\int_{L_\xi}{\mathrm{curl}_{\e}\, \tilde{B}\;dx} \right| = 2.
\end{equation}
Finally, we define the limit layer
\begin{equation}\label{eqn:limitLayer}
L^0_\xi : = (-\xi,1)\times \{0\}.
\end{equation}

In the following, we take an arbitrary subsequence $\e_n \to 0$ and investigate sequences with uniformly bounded energy with an eye towards ultimately understanding the most physically relevant effective energies
\begin{equation} \label{eq:limit}
	E_\xi [k] : = \inf{\{\underset{n \to \infty}{\liminf}\; {E_{\e_n,\xi}[k_n,\tilde{B}_n]}: k_n \chi_{(L_{\e_n,\xi})^c} \to k \quad \mathrm{strongly \;in}\;\; L^2(\o;\mathbb{R}^2)\}}
\end{equation}
and
\begin{align*} \label{eq:limit}
	E_\xi [k,\alpha] : = \inf{\{\underset{n \to \infty}{\liminf}\; {E_{\e_n,\xi}[k_n,\tilde{B}_n]}: k_n \chi_{(L_{\e_n,\xi})^c} \overset{L^2}{\to} k, \quad \mathrm{curl}_{\e_n}\, \tilde{B}_n \wk{L^2} \alpha \}}.
\end{align*}

\subsection{Compactness}

Under the hypotheses (\ref{eqn:BlayerConstraint}) and (\ref{eq:halfdisclin}), we consider any sequence with uniformly bounded energy, and using (\ref{eq:Change_of_Var}) and (\ref{eq:bulk+layer}), we write it as the sum of the non-negative energies
\begin{equation}\label{eq:unifBddEnergy}
 E^{bulk}_{\e_n,\xi}[k_n] + E^{layer}_{\e_n,\xi}[\tilde{k}_n,\tilde{B}_n] = E_{\e_n,\xi}[k_n, B_n] \leq C<\infty.
 \end{equation}

We prove the following theorem.
\begin{theorem}\label{thm:compactness}
    Let $k_n,\tilde{B_n}$ have uniformly bounded energy  as in (\ref{eq:unifBddEnergy}). Then the following hold:
    \begin{itemize}
        \item Outside the layer, $k_n \chi_{(L_{\e_n,\xi})^c} \to k$ strongly in $L^2(\o;\mathbb{R}^2)$, where $k \in SBV(\o;S^1)$. Furthermore, the jump set of $k$ is contained in $\overline{L^0_\xi}.$
        \item In the layer, we can generate a rescaled $k_n$ which is denoted by $\tilde{k}_n$.
    We will have the convergences:
    \begin{align*}
    &\tilde{k}_n \wk{} \tilde{k} \quad \mathrm{weakly\;\; in\;}  L^2(L_\xi;\mathbb{R}^{2}), \\
    &\tilde{B}_n \wk{} B \quad  \mathrm{weakly\;\; in\;} L^2(L_\xi;\mathbb{R}^{2 \times 2}),\\
    &\mathrm{curl}_{\e_n}\, \tilde{B}_n \wk{} \alpha \quad \mathrm{weakly\;\; in\;}  L^2(L_\xi;\mathbb{R}^{2}),
    \end{align*}
    for some $\tilde{k} \in  L^2(L_\xi;\mathbb{R}^{2})$, $\alpha \in L^2(L_\xi;\mathbb{R}^{2})$, and $B := \begin{bmatrix}
	0 & \partial_2 \tilde{k}_1 \\
	0 & \partial_2 \tilde{k}_2
\end{bmatrix}.$
\item Further, for almost every $x_1 \in (-\xi,1)$, $\tilde k(x_1,\cdot) \in W^{1,2}(-\frac{\xi}{2},\frac{\xi}{2})$ and we may define $[\tilde{k}]: (-\xi,1) \to {\mathbb{R}^2}$ as
\begin{equation}\label{def:tildekjump}
[\tilde{k}](x_1) := \int_{-\frac{\xi}{2}}^{\frac{\xi}{2}}{\partial_2\tilde{k}\;dx_2}.
\end{equation} The following compatibility condition between $\tilde{k}$ and $\alpha$ holds:
	$$[\tilde{k}](s) = \int_{-\xi}^{s}\int_{-\frac{\xi}{2}}^{\frac{\xi}{2}} \alpha \;dx,\qquad\qquad |[\tilde{k}](1)| = 2. $$
    \end{itemize}
\end{theorem}

We note that even though $k$ and $\tilde{B}$ were independent in the beginning, these fields become coupled through $\tilde k$ in the limit. In Proposition \ref{prop:traceCoincidence}, we will show that this relation can be directly expressed without the intermediate field $\tilde{k}$.

\begin{proof}
\textbf{Step 1: Bulk Energy.}
To control the energy in the bulk, we observe that
\begin{equation}\label{eq:Uniform_bds_Bulk}
	\sup_n \int_{\o\setminus{L_{\e_n,\xi}}}{\left[\frac{(|k_n|-1)^2}{\e_n\xi^2} + |\nabla k_n|^2 \right]\; dx} \leq C .
\end{equation}
In particular, for any $U$ smooth open set which is compactly contained in the set $\o\setminus\overline{L^0_\xi}$ (recall (\ref{eqn:limitLayer})), we have that $$\sup_n\|k_n\|_ {W^{1,2}(U;\mathbb{R}^2)}\leq C.$$
Thus, up to a subsequence (not relabeled), we have that $k_n \to k$ strongly in $L^2(\o\setminus\overline{L^0_\xi};\mathbb{R}^2)$.
Because of the unit norm regularization, we have that $|k|=1$ almost everywhere. Furthermore, since $k \in W^{1,2}(\o \setminus\overline{L^0_\xi};\mathbb{R}^2)$, it is an integration by parts argument \cite[Proposition 4.4]{AmbrosioFuscoPallara} to show $k \in SBV(\o; \mathbb{R}^2)$ where the jump set of $k$ is contained in $\overline{L^0_\xi}$ up to a set of $\mathcal{H}^1$-measure zero.

\textbf{Step 2: Layer Energy.}
In this portion of the energy, we have
\begin{multline}\label{eq:Uniform_bds_Layer}
    	\int_{L_\xi} \left[\frac{(|\tilde{k}_n|-1)^2}{\xi^2}  + \e_n\left|\nabla_{\e_n} \tilde{k}_n - \frac{\tilde{B}_n}{\e_n}\right|^2\right]\;dx \\
    	+\int_{L_\xi}\left[\xi^2|\mathrm{curl}\,_{\e_n}\tilde{B}_n|^2 + \frac{1}{\xi^2}W(\xi |\tilde{B}_n|)\right]\; dx \leq C.
\end{multline}

 Using the quadratic coercivity of $W$ in \eqref{hyp:f1}, we have
\begin{align}\label{eq:bounds_layer}
	&\|\tilde{k}_n\|_{L^2} \leq C, \\
	&\|\e_n \nabla_{\e_n} \tilde{k}_n - \tilde{B}_n  \|_{L^2} \leq C \e_n^\frac{1}{2} \label{eq:k_b_relation},\\
	&\|\mathrm{curl}_{\e_n}\,\tilde{B}_n\|_{L^2} + \|\tilde{B}_n\|_{L^2} \leq C. \label{eq:B_bounds}
\end{align}
This implies that up to a subsequence, not relabeled, we obtain
\begin{align}
&\tilde{B}_n \wk{} B \quad  \mathrm{weakly\;\; in\;} L^2(L_\xi;\mathbb{R}^{2 \times 2}),\\
&\mathrm{curl}_{\e_n}\, \tilde{B}_n \wk{} \alpha \quad \mathrm{weakly\;\; in\;}  L^2(L_\xi;\mathbb{R}^{2}),
\end{align}
for some $B \in L^2(L_\xi;\mathbb{R}^{2 \times 2})$ and $\alpha \in L^2(L_\xi;\mathbb{R}^{2})$.

Furthermore, using the quadratic bounds on $\tilde{k}$ and \eqref{eq:k_b_relation}, we deduce that
\begin{align}
	&\tilde{k}_n \wk{} \tilde{k} \quad \mathrm{weakly\;\; in\;}  L^2(L_\xi;\mathbb{R}^{2}), \label{eq:wk_k_tilde}\\
	&\e_n \nabla_{\e_n} \tilde{k}_n \wk{} B \quad  \mathrm{weakly\;\; in\;} L^2(L_\xi;\mathbb{R}^{2 \times 2}), \label{eq:wk_dgrad}
\end{align}
for some $\tilde{k} \in  L^2(L_\xi;\mathbb{R}^{2})$.

In order to further characterize $B$, we can analyze component-wise for $i = 1,2$ using \eqref{eq:wk_dgrad}. To be precise, for $\phi \in C^\infty_c(L_\xi)$, it follows that
\begin{align*}
    &\int_{L_\xi} B^{i1} \phi\; dx= \underset{n \to \infty}{\lim}\;\int_{L_\xi} \e_n \partial_1 \tilde{k}_n^{i} \phi \;dx = - \underset{n \to \infty}{\lim}\;\e_n \int_{L_\xi} \tilde{k}^i_n \partial_1\phi\;dx = 0, \\
    &\int_{L_\xi} B^{i2} \phi\; dx= \underset{n \to \infty}{\lim}\;\int_{L_\xi} \partial_2 \tilde{k}_n^{i} \phi \;dx = - \underset{n \to \infty}{\lim}\; \int_{L_\xi} \tilde{k}^i_n \partial_2\phi\;dx = \int_{L_\xi} \partial_2 \tilde{k}^i \phi\;dx,
\end{align*}
where we have applied \eqref{eq:wk_k_tilde} after integrating by parts.

Thus, we conclude that
$$B = \begin{bmatrix}
	0 & \partial_2 \tilde{k}_1 \\
	0 & \partial_2 \tilde{k}_2
\end{bmatrix}.$$

Now we can get information on $\mathrm{curl}_{\e_n}\, \tilde{B}_n$ by integrating by parts. We will do it component-wise for $i = 1,2$. For any $\phi \in C^\infty(L_\xi)$ not necessarily compactly supported, we have
\begin{align*}
	&\int_{L_\xi}{\alpha^i \phi\; dx} = \underset{n \to \infty}{\lim}\; \int_{L_\xi}{[\mathrm{curl}_{\e_n}\, \tilde{B}_n]^i\;\phi\; dx},\\
	 & = \underset{n \to \infty}{\lim}\; - \int_{L_\xi}{\left[\tilde{B}^{i2}_n\partial_1\phi-\frac{1}{\e_n}\tilde{B}^{i1}_n\partial_2 \phi\right]\; dx} + \int_{\partial L_\xi} \phi [\tilde{B}^{i2}_n\nu_1 -\frac{1}{\e_n}\tilde{B}^{i1}_n\nu_2 ]\; d\mathcal{H}^1,
\end{align*}
where $\mathcal{H}^1$ is the one dimensional Hausdorff (surface) measure in $\mathbb{R}^2$.

Using the tangential relations in \eqref{eq:B_tan} and the weak convergence of $\tilde{B}_n$, we simplify
\begin{multline}\label{eqn:smooth_int_parts_alpha}
	\int_{L_\xi}{\alpha^i \phi\; dx} + \int_{L_\xi}{\partial_2 \tilde{k}^i \partial_1 \phi}\;dx \\= \underset{n \to \infty}{\lim}\; \left[\int_{L_\xi}{\frac{1}{\e_n}\tilde{B}_{i1}^n\partial_2 \phi\; dx} + \int_{-\frac{\xi}{2}}^{\frac{\xi}{2}} \phi(1,x_2) \tilde{B}^{i2}_n(1,x_2) dx_2\right].
\end{multline}
Taking $\phi \equiv 1$ leads to the relation
\begin{equation}\label{eqn:alpha_avg}
	\int_{L_\xi}{\alpha^i \; dx} = \underset{n \to \infty}{\lim}\;  \int_{-\frac{\xi}{2}}^{\frac{\xi}{2}} \tilde{B}^{i2}_n(1,x_2) dx_2 .
\end{equation}
Allowing $\phi\in C^\infty((-\xi,1))$, which means that $\partial_2 \phi = 0$, leads to the equation
\begin{align}\label{eqn:partial20}
	\int_{L_\xi}{\alpha^i \phi\; dx} + \int_{L_\xi}{\partial_2 \tilde{k}^i \partial_1 \phi\;dx}= \underset{n \to \infty}{\lim}\;  \phi(1)\int_{-\frac{\xi}{2}}^{\frac{\xi}{2}}\tilde{B}^{i2}_n(1,x_2) dx_2 .
\end{align}
Define $[\tilde{k}^i]: (-\xi,1) \to \mathbb{R}$ as $$[\tilde{k}^i](x_1) := \int_{-\frac{\xi}{2}}^{\frac{\xi}{2}}{\partial_2\tilde{k}^i\;dx_2}.$$
This is well-defined as an $L^2$ function since $\partial_2 \tilde{k}^i \in L^2(L_\xi)$.
Using \eqref{eqn:alpha_avg} and the fact that $\phi$ only depends on $x_1$, we can simplify the relation \eqref{eqn:partial20} further as
\begin{align}\label{eqn:partial21}
	\int_{L_\xi}{\alpha^i \phi\; dx} + \int_{-\xi}^1{[\tilde{k}^i]\partial_1 \phi\;dx_1}=  \phi(1)\int_{L_\xi}{\alpha^i \; dx} .
\end{align}
In particular, taking $\phi \in C_c^\infty((-\xi,1))$, we deduce that
$[\tilde{k}^i] \in W^{1,2}((-\xi,1))$ with
$$\frac{d}{dx_1}[\tilde{k}^i](x_1) = \int_{-\frac{\xi}{2}}^{\frac{\xi}{2}} \alpha^i \;dx_2.$$
Now, for generic $\phi\in C^\infty((-\xi,1))$, we can integrate by parts in \eqref{eqn:partial21} to get
\begin{align}
	&[\tilde{k}^i](1)\phi(1)-[\tilde{k}^i](-\xi)\phi(-\xi) = \phi(1)([\tilde{k}^i](1)-[\tilde{k}^i](-\xi)),\\
	&\text{and so  }\qquad [\tilde{k}^i](-\xi)(\phi(1)-\phi(-\xi)) = 0.
\end{align}
Since the equation has to hold for every such $\phi$, we have that $[\tilde{k}^i](-\xi) = 0.$
This gives us a complete characterization of the vertical jump as
\begin{equation}
	\label{eqn:k_jump}
	[\tilde{k}^i](s) = \int_{-\xi}^{s}\int_{-\frac{\xi}{2}}^{\frac{\xi}{2}} \alpha^i \;dx.
\end{equation}

By our convergences, we also have that \eqref{eq:curl_contraint_cov} passes to the limit. To be precise,
\begin{equation}\label{eq:curl_jump}
	2 = \left|\int_{L_\xi}{\mathrm{curl}_{\e}\, \tilde{B}\;dx} \right| \to \left|\int_{L_\xi}{\alpha \;dx} \right| = \left|[\tilde{k}](1)\right|  .
\end{equation}
\end{proof}

\subsection{Coincidence of traces}
We show that the compatibility relation from Theorem \ref{thm:compactness} relates the jump of $k$ directly to the limit defect field $\alpha.$
\begin{proposition}\label{prop:traceCoincidence}
Supposing that $k,\tilde k$ arise as in Theorem \ref{thm:compactness}, then the compatibility relation
$$\jump{k} = [\tilde{k}]  \quad \text{ for } \mathcal{H}^1\text{-almost every }x\in (-\xi,1)\times \{0\}$$
is satisfied, where $\jump{k}$ denotes the $BV$ jump, oriented as the trace from $\{x_2>0\}$ minus the trace from $\{x_2<0\}$.
\end{proposition}

\begin{proof}

We begin by noting that the function $k$ is such that $$Dk = \nabla k \mathcal{L}^2 + \jump{k}\otimes e_2 d\mathcal{H}^{1}\llcorner_{(-\xi,1)\times\{0\}}.$$
From this it follows that for $i=1,2$ and $\phi \in C_c^\infty (\Omega)$, we have
\begin{equation}\nonumber
\int_{\Omega} \partial_2 \phi k^i\, dx = - \left(\int_\Omega \phi \partial_2 k^i \, dx + \int_{(-\xi,1)\times\{0\}} \phi\jump{k^i} \, \mathcal{H}^{1} \right).
\end{equation}
Consequently, to prove the theorem, it suffices to show that for all $\phi \in C_c^\infty (\Omega \setminus \{(-\xi,0)\})$ it holds that
\begin{equation}\label{eqn:derivativeClaim}
\int_{\Omega} \partial_2 \phi k^i\, dx = - \left(\int_\Omega \phi \partial_2 k^i \, dx + \int_{(-\xi,1)\times\{0\}} \phi [\tilde{k}]^i \, \mathcal{H}^{1} \right).
\end{equation}

Let $\{k_n\}$ be the sequence from Theorem \ref{thm:compactness}. From here, we drop the superscript $i$ but continue to operate component-wise. Directly by the uniform $L^2(\Omega;\mathbb{R}^2)$ bound on $k_n$ and strong convergence away from $(-\xi,1)\times\{0\},$ we have that
\begin{equation}\label{eqn:limitBeforeIBP}
\int_{\Omega} \partial_2 \phi k_n\, dx \to \int_{\Omega} \partial_2 \phi k\, dx \quad\text{ as }n \to \infty.
\end{equation}
Performing an integration by parts, we also find
\begin{equation}\label{eqn:integralSplit}
\int_{\Omega} \partial_2 \phi k_n\, dx =  - \int_{\Omega} \phi  \partial_2 k_n\, dx = - \left(\int_{\Omega \setminus L_{\e_n,\xi}} \phi  \partial_2 k_n\, dx + \int_{L_{\e_n,\xi}}\phi  \partial_2 k_n\, dx\right).
\end{equation}
By (\ref{eq:Uniform_bds_Bulk}), $\nabla k_n\chi_{\Omega\setminus L_{\e_n,\xi}} \rightharpoonup \nabla k$ in $L^2(\Omega).$ The first term on right-hand side of (\ref{eqn:integralSplit}) converges with
\begin{equation}\label{eqn:bulkLimit}
\int_{\Omega \setminus L_{\e_n,\xi}} \phi  \partial_2 k_n\, dx \to \int_{\Omega } \phi  \partial_2 k\, dx \quad \text{ as }n \to \infty.
\end{equation}

For the second term, we recall that $\tilde{k}_n (y_1,y_2) = k_n(y_1,\e y_2),$ and perform a change of variables (recall $L_{\xi} : = L_{1,\xi}$)
\begin{equation}\label{eqn:layerBlowup}
\begin{aligned}
\int_{L_{\e_n,\xi}}\phi  \partial_2 k_n\, dx =& \int_{L_{\xi}}\phi(y_1,\e y_2) \partial_2 k_n(y_1,\e y_2) \e \, d(y_1,y_2)\\
=&  \int_{L_{\xi}}\phi(y_1,\e y_2) \partial_2 \tilde{k}_n(y)  \, d(y_1,y_2).\\
\end{aligned}
\end{equation}
Given the regularity of $\phi, $ the map $\phi_\e(y) : = \phi(y_1,\e y_2)$ converges strongly in $L^2(L_{\xi})$ to $\phi(y_1,0)$. By (\ref{eq:k_b_relation}) and (\ref{eq:B_bounds}), it follows that $\partial_2 \tilde{k}_n \rightharpoonup \partial_2 \tilde{k}$ in $L^2(L_{\xi})$. Passing to the limit in (\ref{eqn:layerBlowup}) as $n \to \infty$, applying Fubini's theorem, and recalling the definition (\ref{def:tildekjump}) for $[\tilde k]$, we find
\begin{equation}\label{eqn:layerLimit}
\int_{L_{\e_n,\xi}}\phi  \partial_2 k_n\, dx \to \int_{L_{\xi}}\phi(y_1,0)  \partial_2 \tilde{k}\, dy = \int_{(-\xi,1)\times\{0\}} \phi [\tilde{k}] \, d\mathcal{H}^{1}.
\end{equation}

Finally, to obtain (\ref{eqn:derivativeClaim}), we pass to the limit in (\ref{eqn:integralSplit}) using (\ref{eqn:limitBeforeIBP}), (\ref{eqn:bulkLimit}), and (\ref{eqn:layerLimit}), thereby concluding the proposition.
\end{proof}
\begin{corollary}\label{Cor: CurlNablak}
    Letting $k$ be as in Theorem \ref{thm:compactness}, $\mathrm{curl}\, {\nabla k}$ is a finite Radon measure.
\end{corollary}
\begin{proof}
We treat $\mathrm{curl}\, \nabla k$ row-wise.

We start from the relation $\mathrm{curl}\,(Dk^i)=0$ in the sense of distributions. Using the decomposition of derivatives for SBV functions, we can write for any $\phi \in C^\infty_c(\o)$:
\begin{equation}
    \langle{\mathrm{curl}\,(\nabla k^i \mathcal{L}^2)},\phi\rangle = - \langle{\mathrm{curl}\,(\jump{k^i} e_2 \mathcal{H}^{1}\llcorner_{(-\xi,1)\times\{0\}})},\phi\rangle.
\end{equation}
where the brackets denote the duality pairing of distributions.
On the left hand side, by definition,
$$\langle{\mathrm{curl}\,(\nabla k^i \mathcal{L}^2)},\phi\rangle = \langle\mathrm{curl}\,\nabla k^i,\phi\rangle.$$
On the right hand side, we also unwrap the duality and use the area formula to find
\begin{equation}\label{eqn:jump_curl}
    - \langle{\mathrm{curl}\,(\jump{k^i} e_2 \mathcal{H}^{1}\llcorner_{(-\xi,1)\times\{0\}})},\phi\rangle = \int_{-\xi}^1\jump{k^i}(x_1)\partial_1\phi(x_1,0)\;dx_1.
\end{equation}
But by Theorem \ref{thm:compactness} and Proposition \ref{prop:traceCoincidence}, we know that $\jump{k^i} \in W^{1,2}((-\xi,1)$ so we can integrate by parts and achieve
\begin{align}
    \left|\langle{\mathrm{curl}\,(\jump{k^i} e_2 \mathcal{H}^{1}\llcorner_{(-\xi,1)\times\{0\}})},\phi\rangle\right|&= \left|\int_{-\xi}^1\frac{d}{dx_1}\jump{k^i}(x_1)\phi(x_1,0)\;dx_1 \right|\nonumber\\
    &\leq \left\|\frac{d}{dx_1}\jump{k^i}\right\|_{L^2(-\xi,1)}\|\phi\|_{\infty}\leq C\|\phi\|_{\infty}.
\end{align}
Thus combining the previous equations, we conclude that
$$\left |\langle\mathrm{curl}\,\nabla k^i,\phi\rangle\right |\leq C\|\phi\|_{\infty}.$$
Therefore, $\mathrm{curl}\,\nabla k^i$ is a finite Radon measure, and furthermore the integration by parts argument shows that $\mathrm{curl}\,\nabla k^i = - \frac{d}{dx_1}\jump{k^i}e_2\mathcal{H}^{1}\llcorner_{(-\xi,1)\times\{0\}}.$
\end{proof}

\begin{proof}[Proof of Theorem \ref{thm:compactCorollary}]
This follows directly from Theorem \ref{thm:compactness}, Proposition \ref{prop:traceCoincidence}, and Corollary \ref{Cor: CurlNablak}.
\end{proof}

\section{Conjectures about the Limiting Energy}\label{sec:conclusion}
Though Theorems \ref{thm:relax} and \ref{thm:compactCorollary} address the behavior of the energy $E_{\e,\xi}$, they leave the $\Gamma$-limits of in the singularly perturbed regimes ($\e,\xi \to 0$) unresolved. In the constrained setting, the principal challenge to characterize the limiting energy is to understand the coupled term $$\int_{L_\xi} \e_n\left|\nabla_{\e_n} \tilde{k}_n - \frac{\tilde{B}_n}{\e_n}\right|^2 \;dx.$$

If we rewrite this term in a rescaled Helmholtz decomposition \\ $\tilde{B}_n = \nabla_{\e_n} z_n + p_n$, where $p_n$ will be rescaled divergence free, using orthogonality of $p_n$ with respect to $\nabla_{\e_n}z_n$, this coupled term will become
$$\int_{L_\xi} \left[\e_n\left|\nabla_{\e_n} \tilde{k}_n - \frac{\nabla_{\e_n} z_n}{\e_n}\right|^2 + \frac{|p_n|^2}{\e_n} \right] dx.$$
So we see that if this has bounded energy, then we will have $p_n \to 0$ strongly and the limit of $\tilde{B}_n$ comes purely from the rescaled gradient term.

Using this, one can heuristically argue for the structure of the limiting energy as follows. Considering the terms which depend on only $\tilde{B}_n$, we can view the relaxation in the double well-function as similar to the dimension reduction where we fix the so called bending vector $\frac{1}{\e}\partial_2 z_n \wk{} \partial_2 \tilde{k}$. Such a relaxation has been considered in the 3D-2D case in \cite{BouFonMasc}, and gives a cross convex-quasiconvex envelope (see also \cite{FonKindPed}). Using the convergences given in Theorem \ref{thm:compactness}, one could imagine leveraging the limiting structure of $B$ and weak lower semicontinuity of cross convex-quasiconvex envelopes (denoted here by $Q^*$) with respect to rescaled gradients
\begin{align*}
    &\int_{L_\xi}\left[\xi^2|\mathrm{curl}\,_{\e_n}\tilde{B}_n|^2 + \frac{1}{\xi^2}W(\xi |\tilde{B}_n|)\right]\; dx\\
    &\to \int_{L_\xi}\left[\xi^2|\alpha|^2 + \frac{1}{\xi^2}Q^*(W(|\cdot|))(\xi \partial_2\tilde{k}\otimes e_2)\right]\; dx.
\end{align*}
Since the envelope \eqref{eq:limit} we are considering should not depend on $\alpha,\tilde{k}$, we may optimize with respect to a lower bound achieved. In particular recalling definition \ref{def:tildekjump} and Proposition \ref{prop:traceCoincidence}, we see through Jensen's inequality and the definition of cross quasiconvexity a possible lower bound is
\begin{multline}\label{conj:jump}
    \int_{L_\xi}\left[\xi^2|\alpha|^2 + \frac{1}{\xi^2}Q(W(|\cdot|))(\xi \partial_2\tilde{k}\otimes e_2)\right]\; dx \\
    \geq \int_{-\xi}^1{\left[\xi\left|\frac{d}{dx_1}\jump{k}(x_1)\right|^2 + \frac{1}{\xi}Q(W(|\cdot|))(\jump{k})\right]}\; dx.
\end{multline}

In this setting, this equation is quite similar to a Modica-Mortola functional for the vertical jump of the director with a transition layer on the order of the core length $\xi$. This is the type of picture predicted by the numerical experiments in \cite{ZANV}.

However, if the correct energy is as in (\ref{conj:jump}), recalling Remark \ref{rmk:quasiconvex}, we see that the disclination layer $(-\xi,1)\times\{0\}$ allows for any jump of the unit-length director field. To counter this, most likely the energy should be modified so to obtain strong convergence (in at least $L^1$) of $\tilde{B}_n$, so that the quasiconvexification of $W(|\cdot|)$ doesn't occur. In principle, this could be done via the inclusion of a smaller order term $\int_{\Omega} \e^3|{\rm div} B|^2 \, dx$ within the energy. Given the smaller order of $\e$, it would most likely be negligible in the limit $\e \to 0$. However, at the $\e>0$ fixed level, it will require that $B=0$ in the sense of traces on the boundary of the layer. This has the effect that $\nabla k$ will still see some small energy contribution within the layer.

\section*{Acknowledgment}
This work was supported by the grant NSF OIA-DMR \#2021019.
L.G was partially supported by the Deutsche Forschungsgemeinschaft 320021702 GRK2326 Energy, Entropy, and Dissipative Dynamics (EDDy). K.S. also acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - GZ 2047/1, Projekt-ID 390685813.
\bibliographystyle{alphaurl}


\section*{NOMENCLATURE}

\section{Introduction}
\IEEEPARstart{I}{nter}-area oscillation (IAO) is common in modern interconnected power transmission networks. For example, 15 oscillation events ranging from 0.1 Hz to 0.8 Hz were observed in the China Southern Power Grid between 2008 and 2012 \cite{liu2018effect}. IAO is undesirable and needs to be effectively suppressed for normal power system operations as they limit power transfer capability and can result in system instability. Especially with the rapid growth of stochastic renewable energy resources, IAO damping control becomes more challenging in recent years \cite{eltigani2015challenges}.

Numerous efforts have been put into designing control strategies to dampen IAO. The Power System Stabilizer (PSS) is one of the most commonly used oscillation suppression devices in power systems, which controls the field excitation of a synchronous generator by providing a supplementary feedback signal \cite{devarapalli2022review}. Existing works often use model-based PSS control schemes such as sparsity-promoting linear quadratic control\cite{6740090}, second-order sliding mode-based damping control \cite{liao2016sliding}, delay-dependent state-feedback, dynamic output-feedback method \cite{li2016design} and $etc.$. More model-based IAO damping methods can be found in reviews \cite{zhang2016review, younis2013wide}. For the implementation of model-based IAO damping control, knowledge of the non-linear dynamics or small-signal model is required. However, the exact
power system model and corresponding parameters may not always be available due to the complex nature of modern power systems, and the time-varying and uncertain operation environment \cite{gupta2021coordinated}.

With the recent breakthroughs in deep learning and broad deployment of phasor measurement units (PMU), Deep Reinforcement Learning (DRL) shows its advanced performances in control and decision-making problems of power systems \cite{chen2022reinforcement}. DRL is a model-free method, only requiring the observable input and output data of the power system. Meanwhile, DRL is robust to the uncertainties of the system, $i.e.$, control signals latency, environment noise, and load variations. A Deep Deterministic Policy Gradient (DDPG) based method and the twin-delayed DDPG method are proposed to overcome various communication delays during damping control \cite{hashmy2020wide, cui2021twin}. To solve the high dimensionality problem of power systems, Mukherjee et al \cite{mukherjee2021scalable} introduce two model reduction approaches for scalable DRL wide-area damping control. For multi-mode oscillation control, a proximal policy optimization (PPO) algorithm for self-tuning of multi-band PSS is proposed \cite{zhang2021novel}. Accounting for $N-1$ and $N-2$ contingencies, a bounded exploratory control-based DDPG (BEC-DDPG) is proposed to coordinate the operation of different PSSs \cite{gupta2021coordinated}. Other related DRL applications include energy storage based oscillation damping via Soft Actor-Critic (SAC) approach \cite{li2021mechanism}, IAO damping using DRL controlled Thyristor Controlled Series Compensator (TCSC) \cite{huang2022damping}, preventing ultra-low-frequency oscillations with asynchronous advantage actor-critic (A3C) algorithm \cite{zhang2020deep}, $etc.$.

For existing DRL-based IAO damping control methods, two challenges remain to be solved. First, how to select the number and location of generators to be controlled are critical for determining the trade-off between control performance and control cost. Selecting all generators in the power system to be controlled under a coordinated manner improves transient response, while high control cost is required \cite{zenelis2018wide}. For a well-studied power system, we can select a specific generator to be controlled by expert experience \cite{9138480} to reduce the control cost. However, expert experience may not always be available. Meanwhile, more generators to be controlled by the DRL-based controller indicated a larger action set needs to be trained during the training process, which increases the computational burden. Second, many DRL-based controllers are designed under the linearized small-signal model of power systems, whose adaptability for practical non-linear power systems remains to be validated.

In this paper, a DRL method and related control strategies are introduced for IAO damping control. The main contributions of this paper are summarized as follows:
\begin{itemize}
	\item A DDPG-based DRL method is proposed to generate wide-area control signals for effectively suppressing IAOs.

	\item A subset of generators for the DRL-based control is selected by oscillation mode analysis and participation factors to balance the control performance and control cost.

	\item A switching control strategy between the PSSs and DRL controllers is proposed to improve the transient response.

	\item The robustness of the proposed DRL method against communication delay is illustrated. In addition, we demonstrate the transfer ability of the DRL model for the non-linear power systems model environment.
\end{itemize}

The remainder of the paper is organized as follows. Section \ref{sec1} introduces the power system model and problem setup of IAO damping control. Section \ref{sec2} presents the DRL method for IAO damping control. Numerical simulations of the IEEE-39 New England power grid model are provided in Section \ref{sec4}. Conclusions are given in Section \ref{sec5}.

\section{Power System Model and Problem Setup}\label{sec1}
In this section, we introduce the power system model for IAO damping control, then we set up the IAO damping control into an optimization problem.
\subsection{Power System Model}\label{sec1-1}
The power system model can be described in a nonlinear, and differential-algebraic form:
\begin{equation}
\boldsymbol{\dot{x}}\left( t \right) =f\left( \boldsymbol{x}\left( t \right) ,\boldsymbol{z}\left( t \right) ,\boldsymbol{u}\left( t \right) ,\boldsymbol{\eta }\left( t \right) \right)
 \label{eq1}
\end{equation}
\begin{equation}
0=g\left( \boldsymbol{x}\left( t \right) ,\boldsymbol{z}\left( t \right) ,\boldsymbol{u}\left( t \right) ,\boldsymbol{\eta }\left( t \right) \right)
 \label{eq2}
\end{equation}
where the dynamic and algebraic variables $\boldsymbol{x}(t) \in \mathbb R ^n $ and $\boldsymbol{z}(t) \in \mathbb R ^s$ represent the state, $\boldsymbol{u}(t) \in \mathbb R ^p$ constitutes the control action and $\boldsymbol{\eta}(t) \in \mathbb R ^q$ is white noise. (\ref{eq1}) accounts for the synchronous generators' electromechanical dynamics and the excitation control equipment. (\ref{eq2}) accounts for load flow, generator stator, and power electronic circuit equations.

Linearizing the power system model in (\ref{eq1}) and (\ref{eq2}) at a stationary operating point and solve the algebraic equations for $\boldsymbol{z}(t)$ will arrive at the linear state-space model:
\begin{equation}
\boldsymbol{\dot{x}}\left( t \right) =\boldsymbol{Ax}\left( t \right)  +\boldsymbol{B}_1\boldsymbol{u}\left( t \right)+\boldsymbol{B}_2\boldsymbol{\eta}\left( t \right)
\label{eq3}
\end{equation}
where $\boldsymbol{A} \in \mathbb R^{n \times n}$, $\boldsymbol{B_1} \in \mathbb R^{n \times p}$, and $\boldsymbol{B_2} \in \mathbb R^{n \times q}$. The system's state variable $\boldsymbol{x}\left( t \right) $ at time $t$ is defined as
\begin{equation}
\begin{cases}
	\boldsymbol{x}\left( t \right) =\left[ \boldsymbol{\theta }\left( t \right) ,\boldsymbol{\omega }\left( t \right) ,\boldsymbol{x}\left( t \right) ^{rem} \right] ^T\\
	\boldsymbol{\theta }\left( t \right) =\left[ \theta _1\left( t \right) ,\theta _2\left( t \right) ,\cdots ,\theta _{n_g}\left( t \right) \right]\\
	\boldsymbol{\omega }\left( t \right) =\left[ \omega _1\left( t \right) ,\omega _2\left( t \right) ,\cdots \omega _{n_g}\left( t \right) \right]\\
	\boldsymbol{x}^{rem}\left( t \right) =\left[ x_{1}^{rem}\left( t \right) ,x_{2}^{rem}\left( t \right) ,\cdots ,x\left( t \right) _{n-n_g}^{rem} \right]\\
\end{cases}
\label{eq4}
\end{equation}
where $n_g$ is the number of generators, $\boldsymbol{\theta }^T,\boldsymbol{\omega }^T\in \mathbb{R}^{n_g}$ are the rotor angles and frequency deviations concerning the normal frequency, and $\boldsymbol{x}_{rem}\in \mathbb{R}^{n-n_g}$ are the remaining state variables.

For the non-linear power system model (\ref{eq1}-\ref{eq2}) or linear power system model (\ref{eq3}), the IAO damping controller can be described as \cite{6740090}:
\begin{equation}
\boldsymbol{u}\left( t \right) =\boldsymbol{u}_{loc}\left( t \right) +\boldsymbol{u}_{wac}\left( t \right)
\label{eq5}
\end{equation}

In (\ref{eq5}), the first term is the local control signal $\boldsymbol{u}_{loc}(t)$ from the PSS on each generator, which is designed based on locally measured states. The second term is the wide-area control signal $\boldsymbol{u}_{wac}(t)$ from the DRL controller, which is an additional control signal designed based on integrating information from all generators. In practice, signals from the controller (\ref{eq5}) will adjust the voltages applied at the field windings of the controlled generators. For generator $i$, the local control signal of the PSS in the Laplace domain denotes \cite{2009.0669}:
\begin{equation}
u_{loc,i}(s)=k_i \cdot \frac{T_{w,i}s}{1+T_{w,i}s} \cdot \frac{1+T_{n1,i}s}{1+T_{d1,i}s} \cdot \frac{1+T_{n2,i}s}{1+T_{d2,i}s} \cdot \dot \theta_i(s)
\label{eq6}
\end{equation}
where $\theta_i(s)$ is the rotor angle of generator $i$ and $k_i$, $T_{w,i}$, $T_{n1,i}$, $T_{d1,i}$, $T_{n2,i}$, and $T_{d2,i}$ are controller gains of the PSS.

For wide-area control $\boldsymbol{u}_{wac}(t)$, a feedback control mechanism is introduced:
\begin{equation}
\boldsymbol{u}_{wac}\left( t \right) =-\boldsymbol{Kx}\left( t \right)
\label{eq7}
\end{equation}
where $\boldsymbol{K}\in \mathbb{R}^{p\times n}$ is the feedback gain matrix.
\subsection{Problem Setup}\label{sec1-3}
The physical nature of power system oscillations is determined by the spectrum of close-loop system matrix $\boldsymbol{A}-\boldsymbol{B}_1\boldsymbol{K}$, $i.e.$, the eigenvalues of $\boldsymbol{A}-\boldsymbol{B}_1\boldsymbol{K}$ reveal the natural frequency and damping ratio of each oscillation mode:
\begin{equation}
f_i=\sqrt{\mathrm{Re}\left( \lambda _i \right) ^2+\mathrm{Im}\left( \lambda _i \right) ^2}
\label{naturalf}
\end{equation}
\begin{equation}
\zeta _i=-\frac{\mathrm{Re}\left( \lambda _i \right)}{f_i}
\label{dampingr}
\end{equation}
where $\lambda _i$ is the $i$th pair eigenvalues of matrix $\boldsymbol{A}-\boldsymbol{B}_1\boldsymbol{K}$, $\mathrm{Re}\left( \cdot \right)$ and $\mathrm{Im}\left( \cdot \right)$ denote the real and imaginary parts of $\lambda _i$, respectively. (\ref{naturalf}) is the natural frequency of $i$th oscillation mode, describing how fast the motion oscillates. (\ref{dampingr}) is the damping ratio of $i$th oscillation mode, describing how much amplitude decays per oscillation. Accordingly, we set the problem of IAO damping control into a framework for the optimization of eigenvalues distributions:

\begin{subequations}
\begin{align}
\begin{aligned}
\mathrm{minimize}\,\,J\left( \boldsymbol{K}(t) \right) &= \alpha\sum_i{\left[ \mathrm{Re}\left( \lambda _i \right) ^2-\mathrm{Re}\left( \hat{\lambda}_i \right) ^2 \right]}\\
&+\beta\sum_i{\mathrm{Im}\left( \lambda _i \right) ^2}
\label{eq8a}
\end{aligned}
\end{align}
\,\,\,\,\,\,\,    \text{subject to}
\begin{align}
\text{dynamics:}\,\,\,\,\, &\boldsymbol{\dot{x}}(t)=\boldsymbol{Ax}(t)+\boldsymbol{B}_1\boldsymbol{u}(t)+\boldsymbol{B}_2\boldsymbol{\eta}\left( t \right)
\label{eq8b}
\\
\text{stability:}\,\,\,\,\,\,\,\, &(\boldsymbol{A}-\boldsymbol{B}_1\boldsymbol{K}(t)) \, \text{Hurwitz}
\label{eq8e}
\end{align}

\label{eq8}
\end{subequations}
where $\hat{\lambda}_i$ is the $i$th pair eigenvalues of the open-loop matrix $\boldsymbol{A}$, $\alpha$ and $\beta$ are weight constants. As shown in the objective (\ref{eq8a}), the goal is to minimize the real part differences between eigenvalues of close-loop matrix $\boldsymbol{A}-\boldsymbol{B}_1\boldsymbol{K}$ and eigenvalues of open-loop matrix $\boldsymbol{A}$, as well as the imaginary part of eigenvalues of close-loop matrix $\boldsymbol{A}-\boldsymbol{B}_1\boldsymbol{K}$. Since from (\ref{naturalf}) and (\ref{dampingr}) we can find that when the real part of eigenvalues does not change significantly and the imaginary part of eigenvalues are nearly zeros, the natural frequencies will decrease and the damping ratios will increase. Smaller natural frequencies and larger damping ratios make the oscillations of the close-loop system become slower and amplitude decay per oscillation becomes larger. In this paper, the control gain matrix $\boldsymbol{K}(t)$ serves as the variable to be optimized with the constraints of power system dynamics (\ref{eq8b}), stability criterion (\ref{eq8e}) and implicit constrain of linear feedback control (\ref{eq5}) and (\ref{eq6}). Note that $\boldsymbol{K}(t)$ is time-varying during the damping control.

The optimization problem in (\ref{eq8}) can be further formatted as a Markov Decision Process (MDP), which converts the damping control of IAOs into sequence decision making. The MDP can be represented as a five-tuple$(S, A, P, R, \gamma)$, among which $S$ is the state of the system, $A$ is the action set, $P$ is the state transition probability, $R$ represents the reward for each transition, and $\gamma$ is a discount factor \cite{sutton2018reinforcement}. At time step $t$, the MDP yields:

\textbf{System State:}
The system state corresponds to the power system state variable $\boldsymbol{s}_t=\boldsymbol{x}\left( t \right)$.

\textbf{Action:}
The immediate action is denoted as $\boldsymbol{a}_t = \boldsymbol{K}(t)$.

\textbf{Reward:} The reward takes the negative value of the objective (\ref{eq8a}), $r_t=-J\left( \boldsymbol{K}\left( t \right) \right)$. The goal of the MDP is to find the optimal action $\boldsymbol{a}_t$ at each step, based on which the largest expect reward can be obtained accordingly.

\textbf{State Transition:}
The state transition at time $t$ is a tuple $(\boldsymbol{s}_t,\boldsymbol{a}_t,\boldsymbol{s}_{t+1},r_t)$. The next system state $\boldsymbol{s}_{t+1}$ is determined by the action $\boldsymbol{a}_t$, system state $\boldsymbol{s}_t$, and the state transition probability $\boldsymbol{s}_{t+1}=P(\boldsymbol{s}_t,\boldsymbol{a}_t)$.

\section{DRL Method}\label{sec2}
To solve the MDP introduced in the former section, a model-free DRL method based on the Deep Deterministic Policy Gradient (DDPG) model is used. In the following, we introduce the architecture, performance improvement techniques, and training of DDPG.
\begin{figure*}[h]
  \centering
  \includegraphics[width=6.6in]{fig/DDPG.eps}
  \caption{The architecture of DDPG for IAO damping control. The circles in green, red, and blue represent neurons in the input layers, hidden layers, and output layers, respectively. }
  \label{DDPG}
  \end{figure*}

\subsection{Architecture}
Fig. \ref{DDPG} illustrates the architecture of the DDPG. The DDPG includes two main components interacting with each other: the DRL controller and environment. The DRL controller is composed of the actor network, critic network, and target network. The environment corresponds to the power system. Details of each component are given as follows.

\textbf{Actor network:}
The actor-network is a multi-layer fully-connected neural network mapping the system state $\boldsymbol{s}_t$ to the action $\boldsymbol{a}_t$ by a non-linear actor function $\mu$:
\begin{equation}
	\boldsymbol{a}_t=\mu(\boldsymbol{s}_t)
	\label{eq10}
\end{equation}

Specifically, the actor-network includes three layers: the input layer, the hidden layer, and the output layer. The input of the actor-network is system state $s_t$. The hidden layer $\boldsymbol{h_{a}}$ can be described as:
\begin{equation}
\boldsymbol{h_{a}}=relu(\boldsymbol{W_{ha}}\cdot \boldsymbol{s}_t^T+\boldsymbol{b_{ha}})
\label{eq12}
\end{equation}
where $\boldsymbol{W_{ha}}$, $\boldsymbol{b_{ha}}$, and $relu$ function are weights, biases, and the activation function of the hidden layer, respectively. The $relu$ function reads:
\begin{equation}
relu(x)
\begin{cases}
\ 0 &\ \ \ \ if\ \ x < 0\\
\ x &\ \ \ \ if\ \ x\geqslant0
\end{cases}
\label{eq13}
\end{equation}

The output of the actor-network is action $\boldsymbol{a}_t$:
\begin{equation}
\boldsymbol{a}_t=tanh(\boldsymbol{W_{oa}}\cdot \boldsymbol{h_a}+\boldsymbol{b_{oa}})
\label{eq14}
\end{equation}
where $\boldsymbol{W_{oa}}$, $\boldsymbol{b_{oa}}$, and hyperbolic tangent function $tanh$ are weights, biases, and the activation function of the output layer.

\textbf{Critic network:}
The critic network is a multi-layer fully-connected neural network aiming to estimate the expected reward of the current action and state:
\begin{equation}
	Q^\mu(\boldsymbol{s}_t,\boldsymbol{a}_t)=E_\mu\left[\sum_{k=0}^\infty\gamma^k\cdot r_{t+k} \Big|\boldsymbol{s}_t,\boldsymbol{a}_t\right]
	\label{eq11}
\end{equation}
where $Q^\mu$ is the critic function, $0<\gamma<1$ is the discount factor accounting for the future and immediate rewards, and $E_{\mu}\left[ \cdot \right]$ represents the expected value.

The critic network has three layers: the input layer, the hidden layer, and the output layer.
The input of the critic network includes the system state $\boldsymbol{s}_t$ and the action $\boldsymbol{a}_t$. The hidden layer $\boldsymbol{h_c}$ can be described as:
\begin{equation}
\boldsymbol{h_c}=relu(\boldsymbol{W_{hc}}\cdot (\boldsymbol{s}_t,\boldsymbol{a}_t)^T+\boldsymbol{b_{hc}})
\end{equation}
where $\boldsymbol{W_{hc}}$, $\boldsymbol{b_{hc}}$, and the $relu$ function are weights, biases, and the activation function of the hidden layer.

The output of the critic network is the expected reward $Q_t$ according to the value from the output layer:
\begin{equation}
Q_t=\boldsymbol{W_{oc}}\cdot \boldsymbol{h_c}+\boldsymbol{b_{oc}}
\end{equation}
where $\boldsymbol{W_{oc}}$ and $\boldsymbol{b_{oc}}$ are weights and biases of the output layer. Note that no activation function is used in the output layer of the critic network.

\textbf{Target network:}
Since we also use the critic network to estimate the expected reward, the updating process is prone to divergence, which is undesired. Inspired by \cite{mnih2015human}, one possible solution is creating the target actor and critic networks which are expressed as $\mu'$ and $Q'$ respectively. The parameters of actor and critic networks are donated as $\Theta$, and the parameters of target actor network and target critic networks are donated as $\Theta'$. The two target networks are updated by tracking the learned networks $\Theta'\leftarrow \tau \Theta+(1-\tau) \Theta'$ while $\tau \ll 1$. Thus, the target networks are updated more slowly, which provides steady support for actor and critic networks and then improves the stability of updating significantly \cite{8882242}.

\textbf{Environment:} The power system acts as the environment for DDPG architecture. Phasor Measurement Units (PMUs) are assumed deployed at each generator bus, making the data of rotor angle and frequency of each generator available. The power system provides the data of rotor angles and frequencies as the system states $\boldsymbol{s}_t$ for the DRL controller. In turn, the DRL controller sends the action values $\boldsymbol{a}_t$ as wide-area control signals to the power system.

\begin{figure}[h]
  \centering
  \subfigure[Generator frequency components]{
    \includegraphics[width=0.22\textwidth]{fig/53c.eps}\label{fig_c5a}
  }
  \subfigure[Participation factors]{
    \includegraphics[width=0.22\textwidth]{fig/53p.eps}\label{fig_c5b}
  }
  \caption{Example of the IAO mode and corresponding participation factors.}
  \label{fig_c5}
  \end{figure}

\subsection{Performance Improvement Techniques}\label{sec2-2}

\textbf{Reward:} To ensure the stability constraint in (\ref{eq8e}), we check the maximum eigenvalue $\lambda_{\max}(\boldsymbol{A}-\boldsymbol{B}_1\boldsymbol{K}(t))$ of the close-loop system matrix for the linear state space model (\ref{eq3}) with the IAO damping controller (\ref{eq4} to \ref{eq5}), $i.e.$, the system is unstable when $\lambda_{\max}(\boldsymbol{A}-\boldsymbol{B}_1\boldsymbol{K}(t)) > 0$; otherwise, the system is stable. During the training of DRL, we penalize the instability by a large negative reward $r_t = -300$.

It should be noted that, for the DRL controller, the power system is a black box model, $i.e.$, the eigenvalues shown in objective \ref{eq8a} and the maximum eigenvalue used for check the stability constraint \ref{eq8e}, can not be directly calculated by the DRL controller, while the system states, $i.e.$, the generator frequency deviations and rotor angles, can be observed by the DRL controller. Hence, a data-driven algorithm called dynamic
mode decomposition (DMD) \cite{schmid2010dynamic} is applied for eigenvalues estimation based on the observable system states. The DMD algorithm can be described as follows \cite{tu2013dynamic}:
\begin{itemize}
	\item \textbf{Step 1:} Build the information matrices from the PMU measurements:
	\begin{subequations}
		\begin{align}
			&\boldsymbol{X}\triangleq \left[ \boldsymbol{x}(0),\boldsymbol{x}(1)\cdots ,\boldsymbol{x}(T-1) \right]
			\\
			&\boldsymbol{Y}\triangleq \left[ \boldsymbol{x}(1),\boldsymbol{x}(2)\cdots ,\boldsymbol{x}(T) \right]
		\end{align}
	\end{subequations}
	where vectors $\boldsymbol{x}(t)$, $t=0,1,\cdots T$ denote the observable input and output system states defined in (\ref{eq4}).
	\item \textbf{Step 2:} Compute the singular value decomposition of matrix $\boldsymbol{X}$ as
	\begin{equation}
		\boldsymbol{X}=\boldsymbol{U\varSigma V}^*
	\end{equation}
	where $\boldsymbol{U}$, $\boldsymbol{V}$, and  $\boldsymbol{\varSigma}$ is the matrix of left singular vectors, right singular vectors, and singular values, respectively, and $\left( \cdot \right) ^*$ represents the transposition operation.
	\item \textbf{Step 3:} Define the matrix
	\begin{equation}
		\boldsymbol{F}\triangleq \boldsymbol{U}^*\boldsymbol{YV\varSigma }^{-1}
	\end{equation}
	\item \textbf{Step 4:} Compute the approximate eigenvalues and eigenvectors by
	\begin{subequations}
		\begin{align}
			&\boldsymbol{F\phi }_i=\tilde{\lambda}_i\boldsymbol{\phi }_i
			\\
			&\boldsymbol{\tilde{\varphi}}_i=\boldsymbol{U\phi }_i
		\end{align}
	\end{subequations}
	where $\tilde{\lambda}_i$ and $\boldsymbol{\tilde{\varphi}}_i$ represent the estimated $i$th eigenvalue and the corresponding eigenvector.
\end{itemize}

As shown in Fig. \ref{fig_estimated}, the small Absolute Errors (AEs) between exact eigenvalues and estimated eigenvalues of the IEEE-39 New England power grid model. The AEs is calculated during an oscillation damping control simulation. The average AEs for the real part and the imaginary part of eigenvalues are below 0.10 and 0.13, respectively.

\begin{figure}[h]
	\centering
	\subfigure[Real part AEs]{
		\includegraphics[width=0.22\textwidth]{fig/real.eps}\label{fig_real}
	}
	\subfigure[Imaginary part AEs]{
		\includegraphics[width=0.22\textwidth]{fig/imag.eps}\label{fig_imag}
	}
	\caption{Absolute Errors (AEs) between exact eigenvalues and estimated eigenvalues. (a) shows AEs of the real part and (b) shows AEs of the imaginary part. The lines in blue and red denote the AEs every time step and the average AEs, respectively.}
	\label{fig_estimated}
\end{figure}

\textbf{Training Reward:}
During the training of DRL, usually experience transitions were uniformly
sampled from the experience replay buffer. However, this approach simply replays transitions at the same frequency that they were originally experienced regardless of their significance, which leads to inefficient training. Hence, the Prioritized Experience Replay (PER) \cite{schaul2015prioritized} is used in this paper. Since it is the unideal experience that can bring more significant adjustment to the agent during training, the PER defines the unideal experience as the prioritized target. The prioritized target is automatically obtained every $M$ episodes from the proposed test, so that the agent could try repeatedly on the unideal circumstance with its updating DNNs until the transition become ideal, or in other words, the agent masters the task that it once failed.

\section{Control Strategies}\label{sec3}
\subsection{Generator Selection}\label{sec3-1}
To save the control cost, not all generators, but a subset of generators are selected to be controlled by the DRL controller. Participation factors are used to select the subset of generators. Participation factors are derived from eigenvalues and eigenvectors of the open-loop matrix $\boldsymbol{A}$, illustrating the impact of each generator on a given oscillation mode.

As an example, Fig. \ref{fig_c5} shows an oscillation mode and the corresponding participation factors of the IEEE-39 New England power grid model with ten generators. The IEEE-39 New England power grid model is displayed in Fig. \ref{DDPG}. We denote the generators at bus 30 to 39 as generator 1 to generator 10. In Fig. \ref{fig_c5a}, the oscillation mode is an IAO mode, where generator 10 is swinging against other generators. Fig. \ref{fig_c5b} reflects that generators 4 to generator 7 have higher impact on the IAO mode due to their larger participation factors, and the participation factors of the rest generators are relatively smaller. Hence, for the given IAO mode shown in Fig. \ref{fig_c5a}, generators 4 to generator 7 are selected to be controlled by the DRL controller.

\subsection{Switching Control Strategy}\label{sec3-2}

Recalling (\ref{eq5}), the IAO damping controller has two types of control signals. One is the local control signal from the PSS on each generator, damping the oscillation with information from the individual generator. The other is the wide-area control signal from the DRL controller, damping the oscillation with integrated information from all generators, $i.e.$, the DRL controller generates wide-area control signals based on the rotor angles and frequency deviations of all generators in the power system. In this section, a Switching Control Strategy (SCS) is proposed to make the DRL controller participate in the IAO damping control Intermittently. As shown in Fig. \ref{DDPG}, the SCS is used to decide whether the DRL controller participates in the IAO damping control according to the reward $r_t$. The rule of SCS is shown below:

\begin{equation}
	\left\{
	\begin{array}{lc}
		\boldsymbol{u}(t)= \boldsymbol{u}_{loc}(t)+\boldsymbol{u}_{wac}(t),& P(t)> \hat{r} \\
		\boldsymbol{u}(t)=\boldsymbol{u}_{loc}(t),& P(t)\le \hat{r}\\
	\end{array}
	\right.
	\label{scsrule}
\end{equation}

\begin{equation}
\begin{aligned}
P\left( t \right) &=\kappa _1\sum_{i=1}^{n_g}{\sum_{j=1}^{n_g}{\left[ \omega _i\left( t \right) -\omega _j\left( t \right) \right] ^2}}+\kappa _2\sum_{i=1}^{n_g-1}{\omega _{i}^{2}\left( t \right)}
\\
&+\kappa _3\sum_{i=1}^{n_g-1}{\left[ \theta _i\left( t \right) -\theta _{ref}\left( t \right) \right] ^2}
\end{aligned}
\label{energy}
\end{equation}
where $\hat{r} > 0$ is a given threshold for switching the DRL controller on or off, $\kappa _1, \kappa _2$and $\kappa _3$ are weights, $\theta _{ref}$ is the phase of the reference generator. $P(t)$ is a system energy-like framework, by using which we quantify the control performances.  In (\ref{energy}), the first term quantifies the amplitude of IAOs between each pair of generators, the second and third terms quantify the kinetic and potential energy of generators, respectively. The smaller $P(t)$ is, the better the the performance of transient responses. Specifically, at the beginning of IAO, $P(t)> \hat{r}$ will be satisfied concerning relatively large frequency violation and angular difference. The wide-area controller is switched on for effective suppression of the oscillation, $i.e.$, the control signals from the trained DRL model provide auxiliary effort on fast oscillation damping. However, the wide-area controller requires additional cost for remote communication, since the DRL needs to obtain the global system states. Considering communication cost, when the reward $P(t) \le \hat{r}$, which means the oscillation has been considerably suppressed, the wide-area controller is switched off and no remote communication cost is needed.

\section{Numerical Simulations}\label{sec4}
In this section, the IEEE-39 New England power grid model is used to validate the proposed method. The diagram of the IEEE-39 New England power grid model is illustrated on the right side of Fig. \ref{DDPG}. The power grid model includes 39 buses and 10 generators. The corresponding electrical parameters are obtained from the power system toolbox \cite{207380}. We denote the generators at bus 30 to 39 as generator 1 to generator 10.

The DRL model introduced in Section \ref{sec2} is trained on the Intel(R) Xeon(R) CPU E5-2620 v3 CPU platform with the TensorFlow software environment. During the training, the minibatch size is 32, and the capacity of the replay buffer is 1000. The max episodes and max steps in one episode are set as 5000 and 500, respectively. The discount factor $\gamma$ is 0.95. The learning rate for both actor and critic networks is 0.0001. The white noise $\boldsymbol{\eta}(t)$ is sampled from Gaussian distribution with zero mean and standard deviation 0.01.

For SCS in Section \ref{sec3-2}, we made use of the cumulative system energy-like framework to determine suitable the threshold $\hat{r}$. The cumulative system energy-like framework $\bar{P}(t)$, $a.k.a$, the return, is is the sum of $P(t)$ values at each time point in a control process and the goal of control is to minimize its value. Different threshold $\hat{r}$ gives different switching points between the wide-area DRL controller and the local PSS controller, and in turn, leads to different $\bar{P}(t)$. Fig. \ref{SCS} shows the normalized $\bar{P}(t)$ changes with the different threshold $\hat{r}$. The normalized $\bar{P}(t)$ is obtained by the max-min normalization to characterize the $\bar{P}(t)$ concerning different thresholds on an equal footing. Each blue point represents a normalized $\bar{P}(t)$ of one control process with a given $\hat{r}$. At each threshold, we do controls with 20 different initial states, and the average normalized $\bar{P}(t)$ of the 20 times control is shown by the red curve. The average normalized $\bar{P}(t)$ has a minimum at $\hat{r}\approx 0.06$ and then increases monotonically. Since the control of the proposed model aims to minimize the $\bar{P}(t)$, $\hat{r}=0.06$ is selected as the reward threshold for the SCS rule in (\ref{scsrule}).

\begin{figure}[htb]
	\centering
	\includegraphics[width=3in]{fig/threshold.eps}
	\caption{The normalized cumulative system energy-like framework $\bar{P}(t)$ with different threshold $\hat{r}$. Each blue point represents the normalized $\bar{P}(t)$ for one control. The red curve is the average normalized $\bar{P}(t)$ of 20 different controls using the same threshold $\hat{r}$.}
	\label{SCS}
\end{figure}

\subsection{Performance Evaluation}
As shown in Section \ref{sec3-1}, the performance of the trained DRL model is tested under the IAO mode that generator 1 to generator 9 oscillates against generator 10. According to Section \ref{sec3-1}, the DRL controllers are deployed on a subset of generators including generator 4 to generator 7. The initial conditions are aligned with the eigenvector of the corresponding IAO mode. Generator 10 is set as a reference generator without the damping controller. We compared the proposed DRL method with the conventional local PSSs method. For the local PSS in (\ref{eq6}), the control gains are chosen according to the tuning strategy \cite{2009.0669}.

Fig. \ref{traincurve} shows the learning curve of the DRL model during training. The line in blue shows the average reward for each episode. The line in red shows the average reward for every 20 episodes. The rewards in blue approaching -300 indicate the penalty for unstable states during training. After 4000 episodes, nearly no unstable states exist. The average reward in red starts at a low value and gradually increases. Finally, the average reward keeps a high value close to zero after 4000 episodes, validating the effectiveness of the training process.
\begin{figure}[h]
	\centering
	\includegraphics[width=3.8in]{fig/traincurve.eps}

	\caption{The learning curve of the DRL model during training. The line in blue shows the reward of each episode. The line in red shows the average reward for every 20 episodes.}
	\label{traincurve}
\end{figure}

Fig. \ref{linear} illustrates the time-domain control results of the linearized state-space model (\ref{eq3}) of the IEEE-39 New England power grid model. The first and second row shows the generator frequency deviations and angle differences, respectively. The angle differences reflect the inter-area power transfer, $i.e.$, larger angle differences correspond to larger inter-area power transfer. The first and second column is control results under local PSS control and the proposed DRL method, respectively. For generator frequency deviations, The lines in blue are frequency deviations of generator 1 to generator 9, and the lines in red are frequency deviations of generator 10. The IAO is shown by generator 10 oscillates against other generators. It can be found that compared with local PSS control, the frequency differences between generator 10 and the other generators are smaller, and the frequency deviation of all ten generators decays faster under the control of the proposed method. Although the angle differences of the proposed method are slightly larger around 2s to 4s, they are more effectively suppressed after 4s of the proposed method.

\begin{figure}[h]
	\centering
	\includegraphics[width=3.3in]{fig/results.eps}

	\caption{Time-domain simulation of the linearized state-space model (\ref{eq3}) of the IEEE-39 New England power grid model. The subfigures show the generator frequency deviations and the angle differences resulting from the use of the local PSS control and the proposed method.}
	\label{linear}
\end{figure}

\begin{figure}[b]
	\centering
	\includegraphics[width=3.3in]{fig/P_t.eps}
	\caption{The quantified control performance $P(t)$ of different control methods for the linear power system environment (\ref{eq3}).}
	\label{per}
\end{figure}

Fig. \ref{per} displays the $P(t)$ during the IAO damping control. The solid line in red and the dashed line in green are the $P(t)$ with local PSS control and the proposed DRL method, respectively. At $t=0$ of initial perturbation of IAO, $P(t)$ is the largest. Then the $P(t)$ gradually decays due to the frequency deviations and angle differences are suppressed by the IAO damping control. However, for local PSS control and the proposed DRL method, the decay speed and amplitude per oscillation are different. It can be found that the $P(t)$ under the proposed method decays faster with a larger decay amplitude per oscillation. Hence, from the quantified performance given by the $P(t)$, the proposed method shows advanced control performance.

\subsection{Transfer ability for the nonlinear power system}
\begin{figure}[t]
	\centering

	\includegraphics[width=3.3in]{fig/nonlinear.eps}

	\caption{Time-domain simulation of the full nonlinear differential-algebraic power system environment (\ref{eq1} to \ref{eq2}) of the IEEE-39 New England power grid model. The subfigures show the generator frequency deviations and active power outputs using the local PSS control and our proposed method both without delay and with an 800 ms delay. Initially, the system is at a steady state, and a three-phase fault is triggered at line \{3,4\} at 0.1 s, which is cleared at 0.2 s and the remote end is cleared at 0.5s.}
	\label{nonlinear}
\end{figure}

\begin{figure}[b]
	\centering
	\includegraphics[width=3.3in]{fig/quant.eps}
	\caption{The quantified control performance $P(t)$ of different control methods for the full nonlinear differential-algebraic power system environment (\ref{eq1}-\ref{eq2}).}
	\label{quant}
\end{figure}

In this section, we demonstrate the transfer ability of the proposed DRL method for the more practical nonlinear power system environment. We directly apply the proposed DRL model trained in the linear power system environment (\ref{eq3}) to control the IAO in the full nonlinear differential-algebraic power system environment (\ref{eq1}-\ref{eq2}). It is meaningful to validate the transfer ability of the DRL model applied in the non-linear power system environment. Since in practical power systems, the dynamics are always non-linear. Fig. \ref{nonlinear} shows the control results for the corresponding nonlinear power system. The first and second column is the control result under local PSS control and our proposed method, respectively. The first and second row is the frequency deviations and power outputs, respectively. For frequency deviations, the curves in red correspond to generator 10 and the curves in blue correspond to other generators. A three-phase fault is triggered at line \{3,4\} at 0.1s, which is cleared at 0.2s and the remote end is cleared at 0.5s. It can be found from Fig. \ref{nonlinear} when the DRL model transfers to the non-linear power system environment, it is also effective for damping the IAO. Specifically, compared with the local PSS control, the proposed DRL method suppress the IAO faster and allows the power system to reach a new steady state earlier.

\subsection{Performance Under time Delays}
\begin{figure}[H]
	\centering
	\includegraphics[width=3.3in]{fig/delay.eps}
	\caption{The control performance of different control methods under communication delay. We set the time delay as 350 and 800 ms respectively.}
	\label{delay}
\end{figure}
In this section, we further investigate the robustness of the proposed method against communication delays. In practical, wide-area control is inevitable to face the challenge of time delays. The time delays may arise from different resources such as the long-distance communication channels, latencies and multiple data rates in the SCADA (supervisory control and data acquisition) network, asynchronous measurements, and local processing and operating times \cite{6740090}. Typically, time delays of communication in power systems range from 100ms (fiber-optic cables) to 700ms (satellite link) \cite{9138480}. We set a moderate time delay as 350ms and a  severe time delay as 800ms to test the robustness of the proposed method. Fig. \ref{delay} shows the control performance of the proposed method. Comparing Fig. \ref{nonlinear} (b) and (d) with Fig. \ref{delay}, we observe that the close-loop performance of the proposed method is not affected significantly by 350 and 800ms delays, which ensures the robustness of the proposed method against time delays.

\section{Conclusions}\label{sec5}
A transfer deep reinforcement learning approach with switch control strategy is developed for the optimal inter-area oscillation damping control. Compared with the local PSS control, the faster decay speed of inter-area oscillations and the larger decay amplitude per oscillation can be realized by the control of proposed DRL method. The switch control strategy between the local PSS and DRL controller helps to improve the performance as well as reduce the communication cost. We illustrate that not all generators, but a subset of generators selected by the participation factors needs to be controlled, mitigating the computing burden of the DRL controller.
Considering the practical application, we show that the proposed DRL method is robustness against the communication delays. The performance is slightly drop even for the serve 800ms time delays. In addition, we demonstrate the transfer ability of the proposed method for full nonlinear differential-algebraic power system environment. The DRL model trained in the linear power system environment can be directly transferred to suppress the inter-area oscillation in nonlinear power system environment.

\section*{ACKNOWLEDGMENT}

The authors would like to thank Prof. Florian Dorfler and Dr. Xiaofan Wu for helpful discussions regarding the Power Systems Toolbox.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi
\bibliographystyle{IEEEtran}