\section{Introduction}
There has been a  growing interest in unmanned aerial vehicle (UAV) communications, seen as a promising enabler for the {fifth generation (5G) and beyond-5G networks} \cite{Bor2019}--\cite{Li2019}. The deployment of UAVs, driven by the emergence of new applications, such as aerial security inspection, smart agriculture, and aerial delivery, is expected to continue shaping future breakthrough services. UAVs can also serve as aerial base stations (UAV-BS) in order to support ubiquitous connectivity, particularly in rural areas and disaster areas, and high data rates, when deployed in dense areas. In this case, the UAV-BS will be a source of/or experience interference as in conventional cellular networks. Fortunately, such a problem can be addressed by the additional degrees of freedom provided solely by a UAV, i.e., the design of {three-dimensional (3D)} location and high probability of line-of-sight (LoS) to ground users \cite{Alzenad2018}. Nevertheless, these additional degrees of freedom are not sufficient to provide the best experience to a constantly growing number of ground users. In this respect, downlink multiple access techniques play an important role in realizing the data rate requirements, low latency, and connectivity, without added resources.

{Recently, rate-splitting multiple access (RSMA) has been identified as a highly-reliable and spectrum-efficient multiple access scheme, that is capable of outperforming both non-orthogonal multiple access (NOMA) and space division multiple access (SDMA) \cite{Mao2018}. RSMA relies on the implementation of a linear precoder at the transmitter and successive interference cancellation (SIC) at the receiver. The process starts by dividing user's messages into common and private parts at the transmitter. The common parts of (all and/or subsets of) users are combined together and encoded into a single common stream, while the private parts are encoded into distinct private streams. These streams are superimposed and sent over a multiple-input multiple-output (MIMO) channel. Then, each user decodes the first common stream and recovers its own data. At the receiver, the interference from the decoded common stream is removed using the SIC. This is followed by successive decoding of the next common parts (of users' subsets) and removing them, then by private part decoding, while treating other (common and) private signals as noise.}

{In \cite{Mao2018}, Mao \textit{et al.} showed that RSMA outperforms NOMA and SDMA systems over a wide range of network loads and user deployment scenarios, and that RSMA has a lower-complex receiver design than NOMA. In \cite{Mao_conf2018}--\nocite{Mao2018_2}\cite{Mao2018_3}, it was shown that RSMA provides better energy and spectral efficiency than NOMA and SDMA for several user deployments, and unicast and multicast transmissions. Extensions of these works have been made to other system models, e.g., downlink coordinated multi-point transmission \cite{Mao_conf2019}, cloud-radio access network (C-RAN) \cite{Yu_2019}, multi-user multi-antenna wireless information and power transfer (SWIPT) \cite{Mao2019spawc}, hybrid radar-communication system \cite{xu2020ratesplitting}, and multi-beam satellite communications \cite{yin2020ratesplitting}.
 {In \cite{Rahmati2019}, downlink RSMA was considered from one ground BS to serve several aerial users only, whereas in} \cite{Ahmad2019}, a UAV-assisted C-RAN system, where RSMA parameters were optimized for joint
transmissions from ground and aerial BSs to terrestrial users, was proposed.

}
{UAV-BS placement has been previously considered with multiple access, thus the UAV-BS' altitude and its horizontal location can be either separately or jointly optimized. In \cite{Yin2019}, the authors investigated joint UAV-BS placement, power and time duration allocations for time-division and frequency-division multiple access (TDMA and FDMA). Xia \textit{et al.} found in \cite{Xiao2016} optimal SDMA user grouping and precoding for ground-to-air (user to UAV-BS) communications in order to maximize the achieved sum data rate.
Also, UAV-BS placement and NOMA power allocation were separately solved in \cite{Liu2019} to maximize the system's sum rate, whereas an extension to joint NOMA power allocation, user pairing, and UAV-BS placement has been studied in \cite{Nguyen2019_inrs}.}

{However, the aforementioned works either focused on simple and well-investigated multiple access techniques such as TDMA, FDMA, SDMA, and single-input-single-output (SISO) NOMA, which are not flexible, and are typically designed for specific scenarios (e.g., single-antenna BS, small number of users, etc.) \cite{Yin2019}--\cite{Nguyen2019_inrs}, or ignored UAV-BS placement optimization when using RSMA on the air-to-ground channels \cite{Ahmad2019}. It is obvious that given the recent research interest into RSMA, the latter has not been sufficiently researched in the context of aerial communications, which motivates us to combine the UAV and RSMA technologies to improve wireless data transmissions to ground users. Specifically, we investigate the joint optimization of UAV placement, RSMA precoding, and rate splitting, in order to maximize the weighted sum rate (WSR) of the downlink communication. Since the optimization problem is non-convex, we propose an alternating optimization (AO) method, where the main problem is divided into two sub-problems, namely UAV placement and RSMA parameters optimization,  respectively solved for fixed RSMA configuration (using successive convex optimization) and UAV location (using the weighted minimum mean square error method), then combined into a global algorithm, where the optimization approaches for the sub-problems are iteratively executed.}

The rest of the paper is organized as follows. In Section II, the system model is presented and the problem is formulated. Section III details the proposed solutions {and discusses their complexity}, while Section IV presents the results. Finally, Section V concludes the paper.

\textit{Notations:} In the remainder of the paper, boldface uppercase and lowercase represent matrices and vectors, respectively. $(\cdot )^{T}$ and $(\cdot )^{H}$ denote the transpose and Hermitian operations respectively. $\mathbb{E}(\cdot)$ is the expectation, {$\Re(\cdot)$ is the real part of a complex number}, $||\cdot||$ is the Euclidean norm, $|\cdot|$ is the absolute value, \textbf{\textit{0}} is the zero matrix, {\textbf{\textit{I}} is the identity matrix}, {$\mathbbm{1}_{a \times b}$ is $a \times b$ all-ones matrix, and} tr$(\cdot)$ is the trace of a matrix.

\section{System Model and Problem Formulation}

\subsection{System Model}
{We consider a downlink one layer rate-splitting (RS) transmission where a hovering UAV-BS, equipped with $N_t$ antennas, is communicating with $K$ single-antenna users, using a dedicated frequency band. At the UAV side, the distinct messages of all users denoted by $M_k$ ($k=1,2,...,K$) are initially divided into common and private parts as $M_k=\{M_k^{p}, M_k^{c} \}$. Then, all common messages from all users are encoded together into a single stream $s_{0}$, for the purpose of reducing interference. This stream $s_0$ will be eventually decoded by all users. Meanwhile, the private message of each user is encoded into a separate private stream $s_k$ ($k=1,2,...,K$) that will be decoded by the corresponding user only. Hence, the vector of the $K+1$ streams to be transmitted is denoted by $\textbf{s}=[s_{0},s_{1},s_{2},...,s_K]^T$.
It is to be noted that a generalized RS model can be obtained by following the same steps as in \cite{Mao2018}. For the sake of simplicity, we keep such design out of the scope of this paper, and focus on the simplified version of one common signal only. In order to reduce multi-user interference, precoding $\textbf{P}=[\textbf{p}_{0},\textbf{p}_{1},\textbf{p}_{2},...,\textbf{p}_{K}]_{N_t \times (K+1)}$ is utilized to give the streams appropriate weights at each transmitting antenna. Finally, signals are superimposed and broadcast as $\textbf{x}=[x_1,\ldots,x_{N_t}]^T=\textbf{P}\;\textbf{s}$.}
Hence, the received signal at user $k$ can be written as
\begin{equation}
\small
\label{eq:receiver}
    y_k=\mathbf{h}_k^H \mathbf{P} \mathbf{s} + n_k,\qquad \forall k=\{1,\ldots,K\},
\end{equation}
where $n_k$ is the circularly symmetric complex additive white Gaussian noise (AWGN) with zero-mean and variance $\sigma^2$, and $\mathbf{h}_k=[h_k^1,\ldots, h_k^{N_t}]^T \in {\mathbb{R}^{N_t \times 1}}$ is the air-to-ground channel from the UAV to user $k$.
{We assume that the air-to-ground communication channels are dominated by line-of-sight (LoS) links, thus,} the channel coefficient $h_k^j$ ($j=1, \ldots, N_t$) follows {the free-space path loss model}
expressed by {\cite{Zeng2016}}
\begin{equation}
\small
    \label{eq:Ricechannel}
    h_k^{j}={{d_k}^{-{{\frac{\beta}{2}}}}}

    , \qquad \forall j =\{1,\ldots,N_t\},

\end{equation}
where $d_k=||\textbf{q}-\textbf{q}_k||_2$ is the distance between the UAV and user $k$, $\textbf{q}=[x,y,z]$ and $\textbf{q}_k=[x_k,y_k,z_k]$ are the UAV and user $k$ 3D locations, respectively, {and $\beta=2$ is the free space path-loss factor \cite{Zeng2016}.}

{Following \cite{Mao2018_3,Ahmad2019}, channels are assumed perfectly known at the transmitter and receivers.}
Also, we assume that {$\mathbb{E}\{\textbf{s} \textbf{s}^H\}=\textbf{\textit{I}}$} and that the transmit power budget at the UAV is constrained by $\text{tr}(\textbf{P}\textbf{P}^H)\leq P_t$, where $P_t$ is the UAV's maximal transmit power.

{
The decoding procedure at the $k^{th}$ user is described as follows. Since the  common stream $s_c$ is allocated the highest power, it will be decoded first by treating the rest of the received signal as noise. Then, the $k^{th}$ user will extract its intended information from the common stream. In order to improve the detection of the private stream, each user will apply SIC to remove the effect of the common stream. Finally, the $k^{th}$ user will decode its intended private stream $s_k$ by considering the rest of the other users' private streams as noise. Hence, the received {signal-to-interference-plus-noise ratios (SINRs)} of the common and private streams $s_{0}$ and $s_k$ at $k^{th}$ user, can be given by}
{
\begin{equation}
\small
    \label{eq:SINR_c}
    \gamma_{k}^{c}=\frac{|\mathbf{h}_k^H \mathbf{p}_{0}|^2}{\sum \limits_{i=1}^K |\mathbf{h}_k^H \mathbf{p}_i|^2+\sigma^2}=\frac{d_k^{-2}|\mathbbm{1} \mathbf{p}_0|^2}{\sum \limits_{i=1}^K d_k^{-2}|\mathbbm{1} \mathbf{p}_i|^2|+\sigma^2}, \; \forall k \in \{1,\ldots,K\}
\end{equation}
and
\begin{equation}
\small
    \label{eq:SINRk}
    \gamma_k^{{k}}=\frac{|\mathbf{h}_k^H \mathbf{p}_{k}|^2}{\sum \limits_{\substack{i=1\\i \neq k}}^K |\mathbf{h}_k^H \mathbf{p}_i|^2+\sigma^2}=\frac{d_k^{-2}|\mathbbm{1} \mathbf{p}_k|^2}{\sum \limits_{\substack{i=1\\ i \neq k}}^K d_k^{-2}|\mathbbm{1} \mathbf{p}_i|^2|+\sigma^2}, \; \forall k \in \{1,\ldots,K\},
\end{equation}
where $\mathbbm{1}=\mathbbm{1}_{1 \times N_t}$.}
Let $R_k^j=B \; \text{log}_2\left(1+\gamma_k^j\right)$ be the data rate corresponding to the received SINR $\gamma_k^j$, $\forall k\in \{1,2,...,K\}$, $\forall j\in \{c,1, 2,...,K\}$, where $B$ is the used bandwidth.
Then, in order to ensure that the common message {$s_{0}$} is successfully decoded by all users, the achievable common data rate should not exceed $R_{c}=\text{min}\left\{R_{1}^{c},R_2^{c},...,R_{K}^{c}\right\}$.
We denote {by} $R_{k,\rm{com}}$ the portion of the common rate allocated to user $k$. Then, we have
{$R_{c}=\sum \limits_{k=1}^K R_{k,\rm{com}}$}.
Finally, the overall achievable data rate of user $k$ can be expressed as
\begin{equation}
R_{k,\rm{ov}}=R_{k,\rm{com}}+R_k^{{k}}. \; \; \forall k \in \{1,...,K\}.
\end{equation}

\subsection{Problem Formulation}
In this paper, we jointly optimize the {precoding} matrix $\textbf{P}$, the common rate vector $\textbf{r}=[R_{1,\rm{com}},R_{2,\rm{com}},...,R_{K,\rm{com}}]$, and the UAV {3D} location, with {the aim of} maximizing the weighted sum of overall achievable data rates (WSR), defined as $R_{\rm{ov}}=\sum_{k=1}^K w_k R_{k,\rm{ov}}$, where $w_k$ ($k=1,\ldots,K$) is the weight reflecting $k$th user {traffic} priority\footnote{{For simplicity, we assume that each user is demanding a specific service with a given priority. If the priority is the same for all traffic demands, then WSR becomes the sum data rate ($\forall w_k=1$) or average data rate ($\forall w_k=\frac{1}{K}$).}}. For a given weight vector $\textbf{w}=[w_1,\ldots,w_K]$,
the optimization problem can be formulated as follows:

\begin{subequations}
\small
	\begin{align}
	\max_{\mathbf{P},\mathbf{r},\mathbf{q}} & \quad
	{\text{WSR}=R_{\rm{ov}}=\sum_{k=1}^K w_k R_{k,\rm{ov}}} \tag{P1} \\
	\label{c1}
	\text{s.t.}\quad & \sum_{k=1}^K R_{k,\rm{com}}\leq R_{c} \tag{P1.a} \\
	\label{c2_1} & R_{k,\rm{ov}}\geq R_{k,\rm{th}},\;\forall k \in \{1,\ldots,K\} \tag{P1.b} \\
	\label{c3} & \text{tr}(\mathbf{P}\mathbf{P}^H)\leq P_t  \tag{P1.c}\\
		\label{c4} & \mathbf{r}\geq \mathbf{0}  \tag{P1.d}\\
	\label{c5} & {x_{\rm{min}}\leq x \leq x_{\rm{max}}, y_{\rm{min}}\leq y \leq y_{\rm{max}}, z_{\rm{min}}\leq z \leq z_{\rm{max}}},  \tag{P1.e}
	\end{align}
\end{subequations}
where
$R_{k,\rm{th}}$ is the minimum required rate at user $k$ to ensure respect of QoS, {and ($x_{\min}$,$x_{\max},y_{\min}$,$y_{\max}$,$z_{\min}$,$z_{\max}$) are the minimum and maximum 3D placement coordinates.}
Problem (P1) is highly non-convex. In fact, for a given UAV location, the problem reduces to WSR maximization by optimizing the precoding matrix and rate-splitting for one layer RSMA, as in \cite{Joudeh2016}. It has been shown that the latter problem is non-convex and non-trivial due to the appearance of the precoder weights in the denominator of the SINR equations \cite{Christensen2008}. Thus, by reduction, we deduce that (P1) is non-convex.

\section{Proposed {Solution}}
In order to solve the joint problem in (P1), {we propose a low-complexity iterative algorithm, based on {the alternate optimization of the UAV placement and RSMA parameters.}
}

\subsection{UAV Placement Optimization}
{Assuming} that $\textbf{P}$ and $\textbf{r}$ are given, then (P1) reduces to the following UAV placement problem (P2):

\begin{subequations}
\small
	\begin{align}
	\max_{\mathbf{q},\boldsymbol{\eta}} & \quad {
	{\sum_{k=1}^K w_k \; \eta_k } } \tag{P2}
	\label{c2_21}\\
		\label{c2_0} &
	{\eta_k \leq \text{log}_2 \left(1+{|\mathbbm{1} \mathbf{p}_k|^2}/\left({\sum \limits_{\substack{i=1\\i\neq k }}^K |\mathbbm{1} \mathbf{p}_i|^2 + \sigma^2 ||{\mathbf{q}-\mathbf{q}_k}||_2^{2}}\right)  \right), \; \forall k,} \tag{P2.a}
	\end{align}
\end{subequations}
\begin{subequations}
\small
	\begin{align}
	\label{c2_2} & {||{\mathbf{q}-\mathbf{q}_k}||_2^2 \leq {\frac{1}{\sigma^2} \left( \frac{|\mathbbm{1}\mathbf{p}_k|^2}{\Lambda_k}-\sum_{\substack{i=1\\ i \neq k}}^K|\mathbbm{1} \mathbf{p}_i|^2 \right), \; \forall k}, \tag{P2.b}}\\
	\label{c2_4} & {x_{\rm{min}}\leq x \leq x_{\rm{max}}, y_{\rm{min}}\leq y \leq y_{\rm{max}}, z_{\rm{min}}\leq z \leq z_{\rm{max}}, \tag{P2.c}}
	\end{align}
\end{subequations}
where {$\eta_k \geq 0$ is a SINR slack variable, which respects constraint (\ref{c2_0}), $\Lambda_k=2^{A_k/B}-1$, $A_k=R_{k,\rm{th}}- R_{k,\rm{com}}$, and constraint (\ref{c2_2}) is equivalent to (\ref{c2_1}).
Problem (P2) is non-convex due to the non-convexity of constraint (\ref{c2_0}) with respect to $\textbf{q}$. To handle this issue, we opt for successive convex approximation (SCA) technique, where in each iteration, the right-hand side of (\ref{c2_0}) is replaced by its concave lower bound at a given UAV location denoted by $\textbf{q}^{(l)}$, with $l$ designating the $l^{th}$ iteration. Recalling that any convex
function is globally lower-bounded by its first-order Taylor
expansion at any point, and by following a similar approach as in \cite{Wu2017}, (P2) can be approximated in iteration $l$ by}

{
\begin{subequations}
\small
	\begin{align}
	\max_{\mathbf{q},\boldsymbol{\eta}} & \quad {
	\sum_{k=1}^K w_k \eta_k } \tag{P3}
	\label{c3_}\\
	\label{c3_0} & \eta_k \leq R_{k}^{k,(l)}, \; \forall k \in \{1,\ldots,K \} \tag{P3.a}\\
	\label{c3_1} & {\text{(\ref{c2_2})}-\text{(\ref{c2_4})}, \tag{P3.b}}
	\end{align}
\end{subequations}
where $R_{k}^{k,(l)}$ is the lower bound of the private signal data rate at the $l^{th}$ iteration, expressed by
\begin{equation}
\small
    \label{eq:LB}
    R_{k}^{k,(l)}=-A_k^{(l)} \left( d_k^{2} - \left(d_k^{(l)}\right)^{2} \right) + B_k^{(l)},
\end{equation}
with
\begin{equation}
\small
    \label{eq:Akl}
    A_{k}^{(l)}=\frac{\log_2(e) |\mathbbm{1} \mathbf{p}_k|^2 }{\left( \sum \limits_{\substack{i=1\\i \neq k}}^K |\mathbbm{1} \mathbf{p}_i|^2 + \sigma^2 \left( d_k^{(l)} \right)^{2}\right) \left( \sum \limits_{i=1}^K |\mathbbm{1} \mathbf{p}_i|^2 +  \sigma^2 \left( d_k^{(l)} \right)^{2} \right) },
\end{equation}
\begin{equation}
\small
    \label{eq:Bkl}
    B_{k}^{(l)}=\log_2\left( 1+{|\mathbbm{1} \mathbf{p}_k|^2}/\left({\sum \limits_{\substack{i=1 \\ i \neq k}}^K |\mathbbm{1}\mathbf{p}_i|^2 + \sigma^2 \left( d_k^{(l)} \right)^{2}}\right) \right),
\end{equation}
and $d_k^{(l)}=||\mathbf{q}^{(l)}-\mathbf{q}_k||_2$, $\forall k=1,\ldots,K$.
Hence, problem (P3) is a convex quadratically constrained quadratic programming problem, which can be solved efficiently using existing software tools such as CVX \cite{CVX}.
}

\subsection{RSMA {Precoding} and Rate-Splitting}
Given {a} UAV location $\textbf{q}$,
problem (P1) can be reduced to

\begin{subequations}
\small
	\begin{align}
	\small
	\max_{\mathbf{P}, \mathbf{r}} & \quad
	{\text{WSR}=R_{\rm{ov}}} \tag{P4} \\
	\label{c_P31}
	\text{s.t.}\quad & \text{(\ref{c1})--(\ref{c4})}.  \nonumber \tag{P4.a}
	\end{align}
\end{subequations}
Similar to (P1), problem {(P4)} is also non-convex. To solve it, we adopt the same weighted minimum mean square error (WMMSE) approach as in {\cite{Mao2018,Joudeh2016}}, {where (P4) is transformed into an augmented weighted mean square error (AWMSE) problem (P5) as follows.}
{First, any user $k$ ($k=1,\ldots,K$) detects and estimates $s_{0}$ as $\hat{s}_{0}=e_k^{0}y_k$, where $e_k^{0}$ is the equalizer. After successfully decoding $s_{0}$ and subtracting it from the received signal, $s_k$ can be detected and estimated as
\begin{equation}
\small
\hat{s}_k=e_k^k (y_k-\textbf{h}_k^H \textbf{p}_{0}s_{0})=e_k^k (y_k-d_k^{-1} \mathbbm{1}  \textbf{p}_{0}s_{0}).
\end{equation}
The mean square error (MSE) of each stream can be defined as $\varepsilon_k=\mathbb{E}\{|s_k-\hat{s}_k|^2\}$, calculated as
\begin{equation}
\small
\label{eq:equalizer}
\varepsilon_k^{0}=|e_k^{0}|^2 T_k^{0}- 2 \Re(e_k^{0} d_k^{-1} \mathbbm{1} \textbf{p}_{0})+1
\end{equation}
and
\begin{equation}
\small
\label{eq:equalizer2}
\varepsilon_k^{k}=|e_k^{k}|^2 T_k^{k}- 2 \Re(e_k^{k} d_k^{-1} \mathbbm{1} \textbf{p}_{u})+1,
\end{equation}
where
$T_k^{0}=\sum_{i=0}^K d_k^{-2}|\mathbbm{1} \textbf{p}_i|^2+\sigma^2$
and
$T_k^k=T_k^{0}-|d_k^{-1} \mathbbm{1} \textbf{p}_{0}|^2+\sigma^2$
are the received power at user $k$ ($k=1,\ldots,K$) to decode signals $s_{0}$ and $s_k$, respectively. The optimal MMSE equalizers can then be written as \cite{Joudeh2016}
\begin{equation}
\small
\label{eq:equalizeropt}
(e_k^{0})^{\rm MMSE}= \textbf{p}_{0}^H \textbf{h}_k (T_k^{0})^{-1}\;\text{and}\;(e_k^{k})^{\rm MMSE}= \textbf{p}_{k}^H \textbf{h}_k (T_k^{k})^{-1}.
\end{equation}
By substituting (\ref{eq:equalizeropt}) into (\ref{eq:equalizer})--(\ref{eq:equalizer2}), the MMSEs are written as
\begin{equation}
\small
(\varepsilon_k^{0})^{\rm MMSE}=\min_{f_k^{0}} \varepsilon_k^{0}=(T_k^{0})^{-1} I_k^{0} \; \text{and}\;(\varepsilon_k^{k})^{\rm MMSE}=(T_k^{k})^{-1} I_k^{k},
\end{equation}
where
$I_k^{0}=T_k^{0}-d_k^{-2} |\mathbbm{1} \textbf{p}_{0}|^2$
and
$I_k^{k}=T_k^{k}-d_k^{-2}|\mathbbm{1} \textbf{p}_{k}|^2$
are the interference terms when decoding $s_{0}$ and $s_k$, respectively. Subsequently, the SINRs can be expressed by
\begin{equation}
\small
 \gamma_k^{0}=1/((\varepsilon_k^{0})^{\rm MMSE})-1 \; \text{and}\; \gamma_k^{k}=1/((\varepsilon_k^{k})^{\rm MMSE})-1,
\end{equation}
and the common and private data rates by
\begin{equation}
\small
 R_k^{0}=-\log_2\left( (\varepsilon_k^{0})^{\rm MMSE}\right) \;\text{and}\;   R_k^{k}=-\log_2\left( (\varepsilon_k^{k})^{\rm MMSE}\right),
\end{equation}
respectively. Consequently, the AWMSEs are given by
\begin{equation}
\small
    \label{eq:awmse}
    \zeta_k^{0}= u_k^{0} \varepsilon_k^{0} -\log_2\left( u_k^{0} \right)\; \text{and}\;     \zeta_k^{k}= u_k^{k} \varepsilon_k^{k} -\log_2\left( u_k^{k} \right),
    \end{equation}

where $u_k^{j}>0$ ($j=0,k$), are weights associated with the MSEs of user $k$ ($k=1,\ldots,K$). By setting the optimization variables as the equalizers and weights, the relation data rates--WMMSEs can be given as \cite{Joudeh2016}
\begin{equation}
\small
\label{eq:wm1}
(\zeta_k^{0})^{\rm MMSE}= \min_{\mu_k^{0},f_k^{0}} \zeta_k^{0} =1-R_k^{0}\;\text{and}\;(\zeta_k^{k})^{\rm MMSE}= \min_{\mu_k^{k},f_k^{k}} \zeta_k^{k}=1-R_k^{k}.
\end{equation}
By substituting (\ref{eq:wm1}) into (\ref{eq:awmse}), and after some manipulations
\begin{equation}
\small
\label{eq:mu1}
(u_k^{0})^*=\left((\varepsilon_k^{0})^{\rm MMSE}\right)^{-1}\; \text{and}\;(u_k^{k})^*=\left((\varepsilon_k^{k})^{\rm MMSE}\right)^{-1}.
\end{equation}
}

{Motivated by the data rates-WMMSE relations in (\ref{eq:wm1}), the optimization problem (P4) can be reformulated as
\begin{subequations}
	\begin{align}
	\small
	\min_{\mathbf{P}, \mathbf{v}, \mathbf{u}, \mathbf{e}} & \quad
	\sum_{u=1}^K w_k \zeta_{k}^{\rm tot} \tag{P5} \\
	\text{s.t.}\quad &
	\sum_{k=1}^K v_k^0+1 \geq \zeta_{0},  \nonumber \tag{P5.a} \\
	\label{c2_P5} & \zeta_{k}^{\rm tot}\leq 1-R_{k,\rm th},\;\forall k \in \{1,\ldots,K\} \tag{P5.b}\\
	\label{c3_P5} & \text{tr}(\mathbf{P}\mathbf{P}^H)\leq P_t  \tag{P5.c}\\
		\label{c4_P5} & \mathbf{v}\leq \mathbf{0}  \tag{P5.d}
	\end{align}
\end{subequations}
where $\textbf{v}=[v_1^{0},\ldots,v_K^{0}]=-\textbf{r}$, $\textbf{u}=[u_1^{0},\ldots,u_K^{0},u_1^1,\ldots,u_K^K]$, $\textbf{e}=[e_1^{0},\ldots, e_K^{0}, e_1^1,\ldots,e_K^K]$, $\zeta_k^{\rm tot}=v_k^{0}+\zeta_k^k$ ($k=1,\ldots,K$), and $\zeta_{0}=\max \{\zeta_1^{0},\ldots,\zeta_K^{0}\}$.
}
{When minimizing the objective in (P5) for $\textbf{u}$ and $\textbf{e}$ (fixed $\textbf{P}$ and $\textbf{v}$), the optimal MMSE  ($\textbf{u}^{\rm MMSE}$,$\textbf{e}^{\rm MMSE}$) is obtained
according to (\ref{eq:equalizeropt}) and (\ref{eq:mu1}). The obtained values satisfy the Karush Kuhn Tucker (KKT) optimality conditions in (P5) for $\textbf{P}$. Thus, given (\ref{eq:wm1}) and the common rate transformation $\textbf{v}=-\textbf{r}$, (P5) can be transformed into (P4). Similarly, for any solution ($\textbf{P}^*$,$\textbf{v}^*$,$\textbf{u}^*$,$\textbf{e}^*$) that satisfies the KKT conditions in (P5), the solution ($\textbf{r}^*=-\textbf{v}^*$,$\textbf{P}^*$) satisfies the KKT conditions in (P4). Consequently, (P4) can be transformed into (P5). However, (P5) is also non-convex for joint parameters optimization. To solve it, Mao \textit{et al.} proposed in \cite{Mao2018} to adopt an AO method as follows. In the $l^{th}$ iteration of the AO algorithm, the equalizers and weights are updated using the precoders obtained in the $(l-1)^{th}$ iteration, i.e.,
\[
(\textbf{u},\textbf{e})=(\textbf{u}^{\rm MMSE}(\textbf{P}[l-1]),\textbf{e}^{\rm MMSE}(\textbf{P}[l-1])).
\]
Then, ($\textbf{v}$,$\textbf{P}$) is updated by solving (P5) for ($\mathbf{u}$,$\mathbf{e}$). Hence, ($\textbf{u}$, $\textbf{e}$) and ($\textbf{v}$, $\textbf{P}$) are iteratively updated until convergence of the WSR. The details of the procedure are presented in Algorithm \ref{Algo2}, where $[l]$ is the iteration index, $\mathbf{u}$ is the stream's weight vector, $\mathbf{e}$ is the equalizer vector, $\mathbf{v}$ is the transformation of $\mathbf{r}$, and $\varepsilon \ll 1$ is the convergence condition \cite{Mao2018}.
}

\begin{algorithm}[t]
\small{
\caption{RSMA {Precoding}-Rate Splitting Algorithm}
\label{Algo2}
\begin{algorithmic}[1]
\State {Initialize $l \xleftarrow{}0$, $\mathbf{P}^{(l)}$, ${R}_{\rm{ov}}^{(l)}$}
\Repeat
\State Set $l \xleftarrow{} l+1$; $\mathbf{P}^{(l-1)} \xleftarrow{} \mathbf{P}$
\State Set $\mathbf{u} \xleftarrow{} \mathbf{u}^{\rm{MMSE}}(\mathbf{P}^{(l-1)})$; $\mathbf{e} \xleftarrow{} \mathbf{e}^{\rm{MMSE}}(\mathbf{P}^{(l-1)})$
\State Solve (P5) for updated ($\mathbf{u}, \mathbf{e}$), then update ($\mathbf{P}$, $\mathbf{v}$)
\Until {$|{R}_{\rm{ov}}^{(l)}-{R}_{\rm{ov}}^{(l-1)}|\leq \varepsilon$}
\State Return ($\textbf{P},-\textbf{v}$) and $R_{\rm{ov}}^{(l)}$ \% since $\textbf{v}=-\textbf{r}$ \cite{Mao2018}
\end{algorithmic}}
\end{algorithm}

\subsection{{Joint UAV Placement and RSMA Parameters Optimization}}
{Given the solutions to the independent UAV placement and RSMA precoding-rate splitting sub-problems, we propose here a joint optimization algorithm that solves (P1) iteratively. Indeed, our solution, shown in Algorithm \ref{Algo3}, alternates between solving (P2) for fixed $\textbf{P}$ and $\textbf{R}$ and solving (P4) for a fixed UAV location. This procedure continues until convergence of the WSR, i.e., the difference between WSR performances of the current and previous iterations is below $\varepsilon$.}

\begin{algorithm}[t]
\small{{
\caption{{Joint UAV Placement and RSMA Parameters Optimization Algorithm}}
\label{Algo3}
\begin{algorithmic}[1]
\State {Initialize $l \xleftarrow{} 0$, ($\textbf{q}^{(l)},\textbf{P}^{(l)}, \textbf{u}^{(l)}, \textbf{e}^{(l)}, \textbf{v}^{(l)}, {R}_{\rm ov}^{(l)}$)}
\State Set $l \xleftarrow{} l+1$
\Repeat
\State \multiline{Given ($\textbf{P}^{(l-1)}, \textbf{u}^{(l-1)}, \textbf{e}^{(l-1)}, \textbf{v}^{(l-1)}$), solve (P3) using the SCA method and get the solution $\textbf{q}^{(l)}$}
\State \multiline{Given $\textbf{q}^{(l)}$, solve (P4) using the AO method in Algorithm \ref{Algo2} and get the solution ($\textbf{P}^{(l)}, \textbf{u}^{(l)}, \textbf{e}^{(l)}, \textbf{v}^{(l)}, {R}_{\rm ov}^{(l)}$)}
\State Set $l \xleftarrow{} l+1$
\Until {$|{R}_{\rm{ov}}^{(l-1)}-{R}_{\rm{ov}}^{(l-2)}|\leq \varepsilon$}
\State Return $(\textbf{q}^*, \textbf{P}^*, \textbf{r}^*, R_{\rm{ov}}^*)$
\end{algorithmic}}}
\end{algorithm}

\subsection{{Complexity Analysis}}
{In Algorithm \ref{Algo3}, the complexity of solving problem (P1) resides mainly in the complexity of solving problems (P2) and (P4). The UAV placement problem (P2) is solved using the SCA approach. Since we have $(2K+3)$ constraints in (P3), the required number of iterations by SCA is $\mathcal{O}\left( \sqrt{2K+3} \log_2(1/\varepsilon) \right)$, where $\varepsilon$ is the accuracy of SCA \cite{CVX}. At each iteration, (P3) is solved with complexity $\mathcal{O}\left(Y_1^2 Y_2\right)$, where $Y_1=K+3$ and $Y_2=2K+3$ are the number of variables and constraints, respectively \cite{LOBO1998}. Thus, the overall SCA complexity to solve (P2) is $\mathcal{O}\left(K^{1.5}\log_2\left(1/\varepsilon\right)\right)$. As presented in Algorithm \ref{Algo2}, (P5) can be solved using alternating optimization (AO). For each iteration in Algorithm \ref{Algo2}, the complexity is dominated by step 5, which solves (P5) using the interior-point method. Given $(N_t+1)K$ variables, step 5 has complexity $\mathcal{O}\left( (N_t K)^3 \right)$ \cite{Boyd2004}, and subsequently, the complexity of Algorithm \ref{Algo2} is $\mathcal{O}\left(Y_3 K^3 N_t^3  \right)$, where $Y_3$ is the number of iterations in Algorithm \ref{Algo2}.}
{Finally, the complexity of Algorithm \ref{Algo3} that solves problem (P1) is $\mathcal{O}( Y_4 K^{1.5} \log_2(1/\varepsilon)+ Y_4 Y_3 K^3 N_t^3)$, where $Y_4$ is the number of iterations in Algorithm \ref{Algo3}.}

\section{Numerical Results}
We consider an RSMA-based UAV system, where one $N_t$-antenna UAV-BS is deployed to serve randomly located {$K$ users in an area of $300 \times 300$ m$^2$}.
For the sake of simplicity, we assume that users are on the ground, i.e., $z_k=0$, {$\forall k\in \{1,\ldots,K\}$}, {the noise power $\sigma^2=1$},
{$\textbf{w}=[1,\ldots,1]$ (corresponds to calculating the sum rate)},
{the} bandwidth $B=20$ MHz, and {the data rate threshold} $R_{k,\rm{th}}=0$ Mbps.

\begin{figure}[t]

	\centering
	\includegraphics[width=0.78\linewidth]{Fig1.eps}
	\caption{WSR vs. number of iterations ($K=2$, $N_t=2$, SNR$=\frac{P_t}{\sigma^2}=20$ dB, $z_\min=80$ m, $z_\max=120$ m).}
	\label{Fig:00}
\end{figure}
\begin{figure}[t]

	\centering
	\includegraphics[width=0.88\linewidth]{Fig2.eps}
	\caption{Convergence of the UAV-BS to the best location ($K=2$, $N_t=2$, SNR$=\frac{P_t}{\sigma^2}=20$ dB, $z_\min=80$ m, $z_\max=120$ m).}
	\label{Fig:01}
\end{figure}

{In Figs. \ref{Fig:00}--\ref{Fig:01}, we illustrate the iteration convergence behavior of Algorithm \ref{Algo3} in terms of WSR and UAV-BS location respectively. In Fig. \ref{Fig:00}, we see  that RSMA and SDMA achieve the best performance, with SDMA converging the fastest. Indeed, since the system is underloaded, the number of antennas $N_t$ is sufficient to efficiently handle the multi-user interference. RSMA converges slower than SDMA since it requires more time to adapt its behavior to act like it. Finally, NOMA presents the worst performance since it does not mitigate as efficiently the multi-user interference.
}

{From Fig. \ref{Fig:01}, we see that the UAV-BS converges to a different location, depending on the used multiple access technique.
Indeed, the convergence locations are either on or are very close to the Y-Z plan (where the users are located), and are almost at the same minimum allowed altitude $z_\min$ \cite{cherif2020optimal}.
Nevertheless, RSMA and SDMA favor locations close to one of the users, while the UAV-BS for NOMA is closer to the middle point between them. In fact, RSMA and SDMA look for UAV-BS locations that provide non-degraded channels, whereas NOMA prefers a location with degraded air-to-ground channels to perform best \cite{Mao2018,Chen2016}.
}

\begin{figure}[t]

	\centering
	\includegraphics[width=0.78\linewidth]{Fig3.eps}
	\caption{WSR vs. SNR$=\frac{P_t}{\sigma^2}$ ($K=4$, $N_t=4$, $z_\min=80$m, $z_\max=120$m).}
	\label{Fig:03}
\end{figure}

\begin{figure}[t]

	\centering
	\includegraphics[width=0.78\linewidth]{Fig4.eps}
	\caption{WSR vs. SNR$=\frac{P_t}{\sigma^2}$ ($K=4$, $N_t=4$, $z_\min=80$m, $z_\max=120$m).}
	\label{Fig:02}
\end{figure}

{Fig. \ref{Fig:03} illustrates the WSR performance as a function of SNR for a system of $K=4$ users, with their locations selected as $\textbf{q}_1=[0, 0, 0]$, $\textbf{q}_2=[0, 100, 0]$, $\textbf{q}_3=[150, 150, 0]$, and $\textbf{q}_4=[200, 50, 0]$. The proposed Algorithm \ref{Algo3} is compared to a baseline method, called ``Avg. UAV-BS loc.'', which separates the multiple access problem from the UAV-BS placement one, as in \cite{Liu2019}. Subsequently, we notice that Algorithm \ref{Algo3} outperforms ``Avg. UAV-BS loc.'' for any multiple access scheme. Moreover, RSMA is improved over SDMA and NOMA using our approach, whereas the same performance is achieved by the latter techniques for ``Avg. UAV-BS loc.''. This is due for the optimized UAV-BS locations using Algorithm \ref{Algo3}.}

{Finally, in order to explicitly emphasize the superiority of RSMA over SDMA and NOMA, we adopt for Fig. \ref{Fig:02} the same scenario as Fig. \ref{Fig:03} but with a Rician channel model as follows.}
{The channel coefficient $h_k^j$, ($j=1$ to $N_t$), follows a Rician distribution with parameter $\alpha(\theta_k)$ \cite{Azari2018} and is expressed by
\begin{equation}
\small
    \label{eq:Ricechannel1}
    h_k^{j}={d_k^{-{{\frac{\beta}{2}}}}} {g_k^j(\alpha(\theta_k))}, \qquad \forall j =1,\ldots,N_t,

\end{equation}
where $\theta_k$ is the elevation angle between user $k$ and the UAV, $g_k^j(\alpha(\theta_k))$ is the small-scale fading, which follows a Rician distribution with K-factor $\alpha(\theta_k)$ and $\mathbb{E}\left\{ ||g_k^j||^2 \right\}=1$.}
{Also, $\alpha(\theta_k)=a_1 e^{\theta_k b_1}$,
and ($a_1, b_1$)=(10$^{0.5}$, 10$^{1.5}$)
are constants \cite{Azari2018}.}
{For the sake of simplicity, we solve (P4) for the generated Rician channel coefficients, while (P3) is solved based on the large-scale channels only.
The Rician small-scale coefficients are selected such that the resulting channels are neither aligned, nor orthogonal \cite{Mao2018}.}
{According to Fig. \ref{Fig:02}, RSMA is superior to NOMA and SDMA. Indeed, rate-splitting and precoding, combined to adequate UAV-BS placement, provides non-degraded channels that enable significant WSR gains.
At low SNR, NOMA performs slightly better than SDMA due to the limited power, which performs well over degraded channels. In contrast, SDMA outperforms NOMA at high SNR due to its capability to mitigate interference efficiently.}

\section{Conclusion}
In this paper, we presented {an} in-depth look into the integration of RSMA in UAV{-based} networks. We formulated {the joint UAV placement, beamforming, and rate-splitting} problem, which is non-convex for one UAV-BS. {An alternating optimization } solution was then proposed, where {the UAV placement and the RSMA parameters are optimized iteratively in order to maximize the WSR performance.}
The obtained results illustrate the efficiency and robustness of {the proposed approach} compared to other conventional schemes.
Finally, it is worth noting that {we validate the efficiency of RSMA over NOMA for air-to-ground communications,}
which makes it a promising technology for {beyond-5G non-terrestrial networks}.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\balance

\bibliographystyle{IEEEtran}


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

Recent detections of high-energy, TeV-PeV, neutrinos realized by the IceCube experiment \citep{aartsen2014}
have made more exciting electromagnetic identification of the neutrino sources.
Such high-energy neutrinos are produced by decay of charged pions which are created
through cosmic-ray interactions with radiation ($p\gamma$; \cite{winter2013})
or gas ($pp$; \cite{murase2013}).
Therefore,
detection of high-energy neutrino is a smoking-gun evidence of existence of high-energy protons (cosmic rays).
Observations of high-energy neutrinos
provide unique information about cosmic-ray acceleration mechanisms
and their acceleration sites, if their origins are identified.

In \citet{icecube2015}, detections of 54~neutrino events by the IceCube experiment in total
are reported from the 4-year data.
Arrival directions of the 54~IceCube events are consistent with being
isotropic and do not show any clustering in the Galactic plane \citep{icecube2015},
indicating the neutrino sources would be extragalactic \citep{aartsen2014}.
Neutrinos can travel cosmological distances without being deflected by cosmic magnetic field
nor absorbed by photon field.
On the other hand, contribution of high-redshift sources represents
only a small fraction of the total observed flux due to the redshift dilution.
Thus, the competition between the neutrino's penetrating power and
the redshift dilution makes emissions from sources
at redshift of $z\sim1-2$ dominant to the IceCube neutrino
for single-neutrino ({\it singlet})
events (\cite{kotera2010}; \cite{icecube2017_multiplet}).
The IceCube collaboration started issuing
real-time alerts in April 2016 \citep{aartsen2017realtimealert} \sout{and}
and the number of the neutrino detection is increasing.

Many hypotheses for the origin of the high energy neutrinos have been proposed,
including
blazars \citep{mannheim1995,mucke2003,becker2005},
starburst galaxies \citep{loeb2006,bechtol2017},
Type~IIn supernovae \citep{murase2011,aartsen2015sn2n},
choked-jet supernovae \citep{razzaque2004,senno2016},
gamma-ray bursts (GRBs; \cite{waxman1997,aartsen2017grb}),
tidal disruption events (TDEs; \cite{senno2017}),
active galactic nuclei (AGN) core \citep{eichler1979,inoue2019}.

These theories have been tried to be assessed observationally and some observational
results set constraints on the theories as described in the
this paragraph.
\citet{kadler2016} showed that
a high fluence GeV blazar, PKS~1424-418, is a possible origin for
one of the high-energy starting events (HESEs), HESE-35,
based on the temporal and positional coincidence
between the neutrino detection and
the $\gamma$-ray flare of the blazar.
This blazar shows broad CIV, CIII], and Mg~II lines in its optical spectrum and
is classified as a flat-spectrum radio quasar (FSRQ) at $z=1.522$ \citep{white1988}.
On the other hand, in general,
there is not a good spatial
correlation between the neutrino detections and Fermi $\gamma$-ray blazars,
indicating that contribution of the Fermi $\gamma$-ray blazars
to the diffuse neutrino flux is
as small as $<27$\% \citep{aartsen2017blazar}.
No other origin candidates have been strongly supported
as a significantly contributing source to the neutrino background
\citep{bechtol2017,senno2016,senno2017}.
Very recently,
a radio-emitting TDE at $z=0.051$
discovered on April 1, 2019,
is claimed to be an associated source to
the neutrino event IceCube-191001A \citep{stein2020}.

On 22 September 2017 at 20:54:30.43 (MJD=58018.871),
the IceCube experiment
detected a track event of a high-energy neutrino
($\sim290$~TeV; \cite{icecube2018_science1}),
IceCube-170922A, and
\citet{kopper2017icecubegcn} reported it
via
The Gamma-ray Coordinates Network (GCN) Circular.
The direction is constrained to be
R.A.$=77.43^{+0.95}_{-0.65}$ and Dec.$=+5.72^{+0.50}_{-0.30}$ degrees
in J2000 equinox (90\% confidence region).
TXS~0506+056, a BL~Lac blazar within the error region,
was pointed out to be a good candidate of the counterpart
(\cite{tanaka2017_IC170922a}; see also \S\ref{sec:sec_discnirvari})
and observed with many telescopes over a wide wavelength range
after that \citep{icecube2018_science1}.
During the intensive follow-up observations,
the redshift of TXS~0506+056\ was successfully determined to be $z=0.3365$ \citep{paiano2018}.
In addition, TXS~0506+056\ was detected
in high-energy $\gamma$-ray with
Large Area Telescope (LAT) on the Fermi satellite (\cite{tanaka2017_IC170922a}; \cite{icecube2018_science1}),
MAGIC telescope (\cite{icecube2018_science1}; \cite{ansoldi2018}),
and AGILE $\gamma$-ray telescope \citep{lucarelli2019}, independently,
which strength the coincidence between TXS~0506+056\ and the neutrino source.

TXS~0506+056\ is a blazar
registered in a catalog of
Texas Survey of Radio Sources
\citep{douglas1996} and
one of the highest energy $\gamma$-ray emitting blazars
among blazars detected by the Energetic Gamma Ray Experiment Telescope~(EGRET)
$\gamma$-ray (30~MeV-30~GeV) satellite \citep{dingus2001}.
The radio-to-gamma-ray spectral energy distribution (SED; \cite{icecube2018_science1}) in combination with
its featureless spectra (\cite{paiano2018}; \cite{icecube2018_science1})
and peak frequency of $10^{14.5\pm0.25}$~Hz \citep{padovani2019}
indicates that TXS~0506+056\ belongs to
intermediate synchrotron peaked
($10^{14}<\nu_{S}<10^{15}$~Hz) BL~Lac objects
(ISPs or IBLs; \cite{padovani2018})
although \citet{padovani2019} claim that TXS~0506+056\ is a masquerading BL~Lac.
The bolometric luminosity
$L_{\rm{bol}}$
is estimated to be a few $\times10^{45}$~erg~s$^{-1}$ \citep{padovani2019},
which is roughly consistent with being an ISP
by following the so-called blazar sequence
\citep{fossati1998,kubo1998,ghisellini2017}.

In SEDs of BL~Lacs,
non-thermal emission from a relativistic jet dominates thermal emission
from an accretion disk
in rest-frame UV-optical wavelengths and
from a dusty torus in rest-frame near-infrared (NIR) wavelengths,
as well as its host galaxy.
Temporal flux (luminosity) variability of blazars is
sometimes
explained by a shock-in-jet model \citep{marscher1985} and
useful to see a possible link between
neutrino emission, probably from a relativistic jet,
and electromagnetic emission activities.
To assess if or not TXS~0506+056\ is a special blazar and September 22, 2017 is a special timing,
it is worth examining variability properties of TXS~0506+056.
Blazars, in general, show large and rapid variability in optical wavelengths \citep{bauer2009blazar},
which could be a key to understanding relationship with neutrino emission.
Although it is dependent on blazar types,
intranight variability is significantly detected for
several tens of percent of blazars
and its duty cycles are also several tens of percent
(\cite{paliya2017}; \cite{bachev2017}; \cite{gaur2017}; \cite{rani2011}).
Thus, short time scale variability is also worth being examined.

Strong polarization in optical wavelengths
is one of the unique characteristics of BL~Lacs
(e.g., \cite{mead1990,falomo2014,ikejiri2011,itoh2016}),
which supports an idea that synchrotron emission dominates optical emission
of BL~Lacs.
Temporal changes of polarization degrees and angles have been observed for many blazars \citep{itoh2016,jermak2016}
and these observables are keys to understanding structure and physics in relativistic jets of blazars
\citep{brand1985,visvanathan1998}.

The structure of this paper is as follows.
In section \ref{sec:sec_followup},
we describe our follow-up observations after the alert of the IceCube-170922A,
including imaging and spectroscopic observations in the optical and near-infrared wavelengths.
We describe the observational results and its discussion in section \ref{sec:sec_resultsdiscussion}.
The neutrino direction and TXS~0506+056\ are
inside the Pan-STARRS1 (PS1; \cite{chambers2016}) footprint
and
outside the Sloan Digital Sky Survey (SDSS; \cite{york2000}) footprint.
Cosmological parameters used in this paper are
$\Omega_{M}=0.3, \Omega_{\Lambda}=0.7, H_{0}=70$~km~s$^{-1}$~Mpc$^{-1}$.
All the observing times are specified in UT.

\section{Follow-Up Observations and Data Reduction}\label{sec:sec_followup}

\begin{figure*}
 \begin{center}
	\includegraphics[width=145mm,bb=0 0 981 757]{kwfc_g.pdf}
 \end{center}
\caption{
Coadded $r$-band image of the field taken with KWFC.
North is up and east is left.
The IceCube 90\% containment region \citep{icecube2018_science1} is shown in blue thick ellipse and
its center is indicated as a blue cross.
The seven flat-spectrum radio sources in the preliminary BROS catalog
which we observed first are indicated in blue small circles
and TXS~0506+056\ region is marked as a red rectangle, close to the blue cross.
Zoom-up $1\times1$~arcmin$^2$ view around TXS~0506+056\ is also shown in the top right inset.
The field-of-views of the optical and NIR instruments used in this paper are shown in brown.
For Tomo-e~Gozen, the field-of-view of one sensor is shown.
}\label{fig:fig_imgfield}
\end{figure*}

We started intensive follow-up observations in optical and NIR wavelengths
right after receiving the realtime alert of the event.
In this section,
we summarize our strategy in general to identify an electromagnetic counterpart
of an IceCube neutrino-emitting source (\S\ref{sec:sec_searchstrategy}),
quick observations to search for a counterpart candidate
of the IceCube-170922A\
(\S\ref{sec:sec_optimg}),
and follow-up observations after a likely counterpart, TXS~0506+056, was identified,
including monitoring imaging
and spectroscopic observations
(\S\ref{sec:sec_optimg},\S\ref{sec:sec_optspec}).
Polarimetric observations for TXS~0506+056\ are also described (\S\ref{sec:sec_optnirpol}).

\subsection{Strategy for IceCube-170922A Electromagnetic Counterpart}\label{sec:sec_searchstrategy}

Since the realtime alert system of IceCube was started in April 2016 \citep{aartsen2017realtimealert},
we have organized a strategic optical and NIR follow-up observing group
utilizing the Optical and Infrared Synergetic Telescopes for Education and Research (OISTER; Yamanaka et al. in prep.)
and other Japanese facilities.
Considering multiple possibilities for
transient high-energy neutrino sources,
we adopt multiple observing strategies using our optical and NIR facilities.

In order to test the blazar scenario,
we first select blazar candidates
with flat radio spectra
catalogued in the BROS catalog
\citep{itoh2018bros}
The BROS catalog
collects flat-spectrum radio sources selected using
a combined catalog \citep{degasperin2018} of
NRAO Very Large Array Sky survey (NVSS; \cite{condon1998}; 1.4~GHz) and
TIFR GMRT Sky Survey (TGSS; \cite{intema2017}; 150~MHz) catalogs.
The final criterion on the radio spectral slope to select flat-spectrum objects is
$\alpha>-0.6$ where $f_{\nu}\propto\nu^{\alpha}$
\citep{itoh2018bros}
while the criterion was $\alpha>-0.5$
at the observation times in this paper.
We carry out optical and NIR imaging observations to detect any brightness change
possibly from a neutrino-emitting blazar within an IceCube localization.
For this purpose, 1-2~m class telescopes in the OISTER collaboration
are mainly used.
A fraction of the BROS sources are apparently faint in optical and not detected in the PS1 data.
We search for variabilities of such faint blazar candidates with Hyper Suprime-Cam (HSC; \cite{miyazaki2012})
on the 8.2-m Subaru telescope.

A supernova is also thought to be
one of the candidates
to generate high-energy neutrinos.
To find a supernova which neutrino may originate from,
we also carry out wide-field imaging surveys in optical to
cover a significant fraction of IceCube localization
using wide-field optical imagers
such as Kiso Wide Field Camera (KWFC; \cite{sako2012})\footnote{KWFC was decommissioned in August 2018. }
and Tomo-e~Gozen\ \citep{sako2016,sako2018} on the 1.05-m Kiso Schimidt telescope
and HSC on the 8.2-m Subaru telescope.

For the IceCube-170922A\ event,
we performed optical and NIR observations
with both the strategies.
Our blind survey for supernovae is described in a separate paper
(Morokuma et al., in prep.).

\subsection{Optical and NIR Imaging}\label{sec:sec_optimg}

All the imaging observations are summarized below and in Table~\ref{tab:tab_obsimg}.
Detailed characteristics of the telescopes and instruments are summarized in the table.

\subsubsection{Initial Response to the Alert: Search for a rapidly variable blazar}\label{sec:sec_searchblazars}
We started follow-up imaging observations
$0.8$~days after the alert of IceCube-170922A
Hiroshima Optical and Near-Infrared camera (HONIR; \cite{akitaya2014})
on the 1.5-m Kanata telescope at the Higashi Hiroshima Observatory ($R_C$ and $J$-bands) and
KWFC on the 1.05-m Kiso Schmidt telescope ($g,r$ and $i$-bands).
In subsequent nights,
HSC on the 8.2-m Subaru telescope ($z$-band)
and 0.5-m MITSuME Akeno telescope
(\cite{yatsu2007}; \cite{shimokawabe2008}; \cite{kotani2005};
$g, R_C$, and $I_C$-bands) were also used.
Field-of-views of
KWFC
and HSC
are well suited for effectively covering the localization area given by IceCube as shown in Figure~\ref{fig:fig_imgfield}.

Brightness of a counterpart was unknown and difficult to predict
because of unknown nature of a counterpart and unknown distance to a counterpart.
Then, we adopt multiple exposure times to be
3-180 seconds as summarized in Table~\ref{tab:tab_obsimg}.

There were 7 flat-spectrum radio sources
in a preliminary version of the BROS catalog
within or right outside the 90\% localization of IceCube-170922A\
as shown in Figure~\ref{fig:fig_imgfield} and
as listed in Table~\ref{tab:tab_brossources}.
Five sources among the seven are detected
in the optical archival PS1 DR1 data.
Each of the 7 BROS sources was observed basically one by one
with Kanata/HONIR and MITSuME while
all the seven BROS sources were covered by two KWFC pointings
and two HSC pointings.
As described in \S\ref{sec:sec_discnirvari},
we made a quick data reduction for the data
and a Kanata/HONIR difference image for TXS~0506+056\
revealed that it was fading in a time scale of a day.

\subsubsection{Monitoring for TXS~0506+056}\label{sec:sec_monitoring}
After detecting a rapid NIR variability of TXS~0506+056\ with Kanata/HONIR (described in \S~\ref{sec:sec_discnirvari}),
we continued monitoring TXS~0506+056\ with telescopes described
in \S\ref{sec:sec_searchblazars}.
In addition, we carried out
$J$, $H$, and $K_s$-band simultaneous imaging with
SIRIUS \citep{nagayama2003,nagayama2012} on the 1.4-m Infrared Survey Facility (IRSF) and
$V$-band imaging with the 0.4-m Kyoto telescope.
Exposure times of the single images are 10 and 60~seconds, respectively.
We also took non-filter CMOS imaging data
with Tomo-e~Gozen\ \citep{sako2016,sako2018} on the Kiso Schmidt telescope.
Two~fps (consecutive 0.5-sec exposure with negligible readout time) full-frame readout mode
\citep{sako2018} was adopted.
No filters were used
and
the sensors of Tomo-e~Gozen\ are sensitive
in $\sim350-900$~nm \citep{kojima2018}.

\begin{table*}[!htbp]
  \tbl{Telescopes and Instruments Used for Follow-up Imaging Observations.}{
  \begin{tabular}{cllllll}
      \hline
      Mode & Telescope & Aperture [m] & Instrument & FoV & Filter & $t_{\rm{exp}}$ [sec]\\
      \hline
      M  & Kyoto           & 0.4 & - & 18~arcmin (rectangle)  & $V$         & 60 \\
      SM & MITSuME (Akeno) & 0.5 & - & 28~arcmin (rectangle)  & $g,R_C,I_C$ & 60 \\
      SM & Kiso   & 1.05   & KWFC    & 2.2~deg (rectangle)    & $g,r,i$ & 10,30,60,180 \\
      M  & Kiso   & 1.05   & Tomo-e Gozen  & ($39.7\times22.4$~arcmin$^2$)$\times$4$^{*}$ & No      & 0.5 \\
      M  & IRSF   & 1.4    & SIRIUS  & 7.7~arcmin (rectangle) & $J,H,K_s$ & 10 \\
      SM & Kanata & 1.5    & HONIR   & 6~arcmin (diameter)    & $R_C,J$ & 25-95 ($R_C$), 10-80 ($J$) \\
      SM & Subaru & 8.2    & HSC     & 1.5~deg (diameter)     & $z$     & 3-5\\
	\hline
  \end{tabular}}
  \label{tab:tab_obsimg}
\begin{tabnote}
In the first column, S and M denote
``survey'' and ``monitoring'', respectively.
(*) indicates that
when we took Tomo-e~Gozen\ data,
the instrument
was operated with the limited number (4) of the sensors
before the completion of the instrument with 84~sensors in 2019.
\end{tabnote}
\end{table*}

\begin{table*}[!htbp]
  \tbl{Seven flat-spectrum radio sources observed within or right outside the IceCube 90\%\ error region.
  }{%
  \begin{tabular}{lllllrrrrl}
      \hline
           &    RA  & Dec    & RA    & Dec   & sep.     & $f_{\rm{TGSS}}$ & $f_{\rm{NVSS}}$ &          & \\
      Name & (NVSS) & (NVSS) & (PS1) & (PS1) & [arcmin] & [mJy]           & [mJy]           & $r(PS1)$ & note\\
      \hline
      \input{brossource_kanata_202009a}
	\hline
    \end{tabular}}\label{tab:tab_brossources}
\begin{tabnote}
  These are registered in the preliminary version of the BROS catalog when we started the follow-up observations.
  Coordinates of the sources in the NVSS and PS1 catalogs,
  separation from the IceCube detection center in arcmin unit,
  radio fluxes in 150~MHz (TGSS) and 1.4~GHz (NVSS) in mJy unit,
  PS1 $r$-band Kron magnitude, and some notes are listed.
\end{tabnote}
\end{table*}

\begin{table*}[!htbp]
  \tbl{Telescopes and Instruments Used for Follow-up Spectroscopic Observations.}{
  \begin{tabular}{llllllrlll}
      \hline
      Telescope & Aperture [m] & Instrument & Grism & Filter & $\lambda$~[\AA] & $R$ & $t_{\rm{exp}}$ [sec] & Date (UT) & MJD \\
      \hline
      Nayuta       & 2.0 & MALLS  & 150~l/mm & WG320 & 4700--8600  & $600$ & $900\times6$ & 2017-09-29 & 58025.7 \\
      Kanata       & 1.5 & HOWPol & grism    & NONE  & 4700--9200  & $400$ & $900\times3$ & 2017-09-29 & 58025.7 \\
      Subaru       & 8.2 & FOCAS  & 300B     & SY47  & 4600--8200  & $400$ & $100\times11$& 2017-09-30 & 58026.6 \\
      Subaru       & 8.2 & FOCAS  & VPH850   & SO58  & 7500--10500 & $1200$& $100\times7$ & 2017-10-01 & 58027.6 \\
      Gemini-North & 8.2 & GMOS   & B1200    & NONE  & 3400--4960  & $3700$& $600\times(3+8)+240\times1$ & 2017-11-11,15 & 58071.5\\
	\hline
    \end{tabular}}\label{tab:tab_obsspec}
\begin{tabnote}
$\lambda$ indicates usable wavelength ranges rather than observed wavelength ranges.
$R$ indicates spectral resolutions.
\end{tabnote}
\end{table*}

\subsubsection{Data Reduction of Imaging Data}
\label{sec:sec_dataredimg}
All the data are reduced in a standard manner,
including
bias, overscan, and dark current subtractions if necessary and possible,
flat-fielding,
astrometry with
Optimistic Pattern Matching (OPM; \cite{tabur2007})
or
Astrometry.net \citep{lang2010}
although some differences exist;
e.g., distortion correction is applied only for HSC data.
Details of the data reductions are described
in \citet{morokuma2014} for KWFC,
in
\citet{tachibana2018} for MITSuME.
The HSC data are analyzed with {\it hscPipe} v4.0.5, which is a standard reduction pipeline
of HSC \citep{bosch2017}.
Tomo-e~Gozen\ CMOS sensor data are also reduced with a dedicated pipeline described in \citet{ohsawa2016}.

To see intranight variability of objects,
we measured apparent magnitudes of TXS~0506+056\ and neighboring stars in individual frames.
With these photometry results, we obtain daily-averaged magnitudes to obtain daily light curves of TXS~0506+056\ and neighboring stars.

Photometry was done with SExtractor \citep{bertin1996} and MAG\_AUTO is used.
All magnitudes in optical and NIR wavelengths are measured in AB system.
Optical magnitudes are calibrated
relative to the PS1 Data Release~2 (DR2) catalog (\cite{magnier2016}; \cite{flewelling2016}).
Johnson-Cousins filter data are also calibrated to PS1 data in filters
whose bandpasses are similar
(i.e., PS1~$r$ for $R_C$ and PS1~$i$ for $I_C$).
No filter Tomo-e~Gozen\ data are also calibrated relative to $r$-band PS1 data.
Magnitudes of field stars in NIR wavelengths
are derived from the 2MASS database~\citep{skrutskie2006}
and converted to those in the AB system by following \citet{tokunaga2005}.

To
search for a counterpart (\S\ref{sec:sec_searchblazars}),
we apply an image subtraction method (\cite{alard1998}; \cite{alard2000})
for the imaging data
with reference images of PS1 in optical and 2MASS in NIR, respectively.
Except for Subaru/HSC, all of our optical imaging data are shallower than PS1 images and
the depths of the search are limited by the depths of our data.
On the other hand,
in NIR, 2MASS data are not so deep and
we also performed another subtraction in which
our first (reference) images are subtracted from our new data.

For photometry for the normal images without image subtractions,
we add addtional 3\% errors to the errors measured with SExtractor.
This is because photometry errors based on photon statistics usually
underestimate the errors which should follow Gaussian distribution around the real values for non-variable objects,
for example, due to imperfect flat-fielding procedure.

\begin{figure}
 \begin{center}
  \includegraphics[width=96mm,bb=0 0 466 792]{brosall_honir_astropy.pdf}
 \end{center}
\caption{
HONIR $J$-band images and subtracted images for the seven
flat-spectrum radio sources catalogued in the preliminary BROS
when we started the follow-up observations
(Table~\ref{tab:tab_brossources}).
From left to right,
HONIR $J$-band image on Sep. 23, 2017,
HONIR $J$-band image on Sep. 24, 2017,
subtracted HONIR $J$-band image (Sep. 23 - Sep. 24),
and
PS1 $r$-band image are shown for each BROS source.
Red circles on the PS1 images are NVSS radio locations with radii of 2~arcsec,
which are typical positional errors of NVSS sources.
}\label{fig:fig_discimg_bros}
\end{figure}

\begin{figure*}
 \begin{center}
	\includegraphics[width=165mm,bb=0 0 654 282]{timeline_ic170922a_02.pdf}
 \end{center}
\caption{
Timeline of the follow-up observations from the IceCube alert.
Kanata/HONIR, Kiso/KWFC, and Subaru/HSC observations
were shown in red, blue, and magenta, respectively.
The IceCube neutrino detection, its alert, and the first Fermi/LAT GCN
indicating TXS~0506+056\ is a possible counterpart
by \citep{tanaka2017_IC170922a} are indicated
in solid black at the top panel.
}\label{fig:fig_timeline}
\end{figure*}

\subsubsection{ASAS-SN monitoring data}
\label{sec:sec_asassndata}

ASAS-SN is covering a large fraction of the sky.
Long-term data nicely covering before and after the neutrino detection
are available and significant variablity and brightening of TXS~0506+056\
before the neutrino detection were reported in \citet{franckowiak2017}.
Original ASAS-SN $V$-band data
are available and taken from
ASAS-SN Sky Patrol\footnote{https://asas-sn.osu.edu}
(\cite{shappee2014}; \cite{kochanek2017}).
The data are
calibrated relative to
AAVSO Photometric All-Sky Survey (APASS; \cite{henden2014}).
As mentioned in \citet{icecube2018_science1},
the nearby west star\footnote{separated by 17.079~arcsec,
objID=114830773534352391, (RA, Dec) = (05:09:24.82, +05:41:35.7) in PS1 DR2.}
($g_{\rm{MeanPSF}}=14.7782\pm0.0024$, $r_{\rm{MeanPSF}}=14.4373\pm0.0010$,
in the PS1 DR2 catalog, corresponding to $V=14.587$, which is converted using an equation
in \cite{kostov2018}) flux
contaminates into the target flux
in photometry of
ASAS-SN Sky Patrol
because of the coarse spatial sampling (8~arcsec pixel$^{-1}$)
and large PSF ($\sim15$~arcsec FWHM)
of ASAS-SN.
We estimate the contribution to the TXS~0506+056\ to be 30\% of the nearby star flux
by calculating the contaminating flux of the nearby star
into the aperture of TXS~0506+056\
with the PSF size, shape (circular), and aperture adopted for the ASAS-SN photometry
and
subtract it
from the target flux
to extract only the flux of TXS~0506+056.
We here do not add additional errors to the original error of TXS~0506+056.
We also confirm this correction does not change our conclusions.
The brightness of the nearby west star is
confirmed to be almost constant over our observing period based on our data.

\subsection{Optical Spectroscopy}\label{sec:sec_optspec}

An obstacle to study the neutrino source in detail was uncertain
redshift determination of TXS~0506+056.
Redshift determination of BL~Lac objects is in general difficult
and require a large amount of observing time with a large aperture optical telescope
(\cite{landoni2015a}; \cite{landoni2015b})
even if their images can be easily taken with 1m-class telescopes.
This was the case for TXS~0506+056\ and the redshift had not been reliably determined \citep{halpern2003}.

TXS~0506+056\ was spectroscopically observed in the optical wavelength
and the redshift is
reported to be $z=0.336$ in \citet{ajello2014}.
However, the origin of this redshift determination is unclear from the figure and text in the paper.
MAGIC high-energy $\gamma$-ray detection (\cite{mirzoyan2017}; \cite{mirzoyan2018})
gives another constraint
on the redshift
through the measurement of $\gamma$-$\gamma$
attenuation effect.
In \citet{icecube2018_science1}, several conservative estimates are done with the MAGIC data and
$95$\%~confidence level upper limit on the source redshift is $z<0.98$
while the lowest upper-limit redshift is $z<0.41$.
Since there were
still debates on the redshift determination
before \citet{paiano2018} determined the redsfhit to be $z=0.3365$,
spectroscopic measurement was still important (\cite{steele2017}; \cite{morokuma2017_IC170922a}).

We took optical long-slit spectra with
Medium And Low-dispersion Long-slit Spectrograph (MALLS) on the 2-m Nayuta telescope,
HONIR \citep{akitaya2014} on the 1.5-m Kanata telescope,
Faint Object Camera and Spectrograph (FOCAS; \cite{kashikawa2002}) on the 8.2-m Subaru telescope,
and
Gemini Multi-Object Spectrograph (GMOS; \cite{hook2004}) on the 8.2-m Gemini-North telescope.
The spectral resolutions of the observations are mostly as low as $R\sim1000$ or even lower.
The exposure times are not so long, 1.5~hours at longest.
These are summarized in Table~\ref{tab:tab_obsspec}.
The FOCAS observations and spectra are described and shown in the supplement part of \citet{icecube2018_science1}.

All the spectra are reduced in a standard manner with IRAF,
including bias subtraction, flat-fielding, wavelength calibration,
sky subtraction, and source extraction,
while the Gemini-N/GMOS data is reduced with the Gemini IRAF package.
Wavelength calibration was done using lamp spectra or night sky lines.
Flux calibration is not applied and the obtained 1d spectra are normalized with continua
to see any weak emission and absorption lines.

\subsection{Optical and NIR Polarimetry}\label{sec:sec_optnirpol}

We also conducted optical and NIR polarimetry observations with HONIR \citep{akitaya2014}
at the Cassegrain focus on the 1.5-m Kanata telescope.
HONIR is equipped with a rotatable half-wave plate and a Wollaston prism
which enables us to conduct
simultaneous polarimetry observations
in an optical and an NIR channels.
We used
$R_C$-band and $J$-band in the optical and NIR channels, respectively.

The data were taken on 15~nights,
from $t=8$~days to
$t=213$~days
after the IceCube alert ($t=0$).
Each observation consisted of a set of 4~exposures at half-wave plate position angles
of
$0^{\circ}.0$,
$22^{\circ}.5$,
$45^{\circ}.0$, and
$67^{\circ}.5$.

The data were reduced
following a methodology of data analysis described in \citet{kawabata1999}
to derive polarization degrees and angles.
Instrumental polarization induced by the optical system
within the Kanata telescope and HONIR has been confirmed to be
as small as 0.1-0.2\% \citep{akitaya2014,itoh2017}.

Data taken with a fully-polarizing filter inserted in front of HONIR
are used for the correction of the wavelength-dependent origin points
of the position angles (originated from the multi-layered superachromatic
half-wave plate).
The obtained polarization degrees are
as stably high as $\gtrsim99$\%\
with this filter in both $R_C$ and $J$ bands
and we do not perform the depolarization correction.
A strongly polarized star, BD+64d106 \citep{schmidt1992},
was observed to correct
observed position angle in the celestial coordinate
and is used to calibrate the position angles of TXS~0506+056.
Galactic foreground polarization should be almost negligible
($P_{R}\lesssim 0.7$\%\ or $P_{J}\lesssim0.2$\%)
because the interstellar extinction toward TXS~0506+056\
is $A_{R}=0.235$ or $A_{R}=0.077$~mag \citep{schlafly2011,serkowski1975}.
These values are mostly smaller than the measured values
and then we do not adopt any correction for the observed polarization.

Observing epochs are separated into two as below.
The first epochs are $\sim1.5$~months after the IceCube alert.
The second epochs are around
$t\sim180$~days after the alert,
when larger polarization degrees of TXS~0506+056\
based on polarimetric data taken with Liverpool/RINGO
are reported \citep{steele2018}.
Motivated by this report, we also
took additional polarization data with Kanata/HONIR.

\section{Results and Discussions}\label{sec:sec_resultsdiscussion}

\subsection{Discovery of Rapid NIR Variability of TXS~0506+056}\label{sec:sec_discnirvari}

We examined the subtracted $J$-band HONIR images
(HONIR-HONIR and HONIR-2MASS) to see any rapid variability of
the BROS blazars and blazar candidates.
Figure~\ref{fig:fig_discimg_bros} shows
HONIR-HONIR subtraction images
for the 7 sources
in Table~\ref{tab:tab_brossources}, which were listed
in the preliminary BROS catalog at the observation time.
We found that TXS~0506+056\ showed a fading trend, by $\sim0.15$~mag,
from Sep. 23 to Sep. 24 in 2017
as shown in the top panels of Figure~\ref{fig:fig_discimg_bros}.
In $g$-band data taken with Kiso/KWFC, about $0.15$~mag decline was also detected.
For the other sources, we did not find any significant NIR variability
from the 2-night HONIR data
or did not detect the object, which was consistent with
non-detection or faint magnitudes recorded in the PS1 optical imaging data.

This rapid brightness change of TXS~0506+056\
may indicate a possible relation
with the neutrino detection,
motivating examining {\it Fermi} $\gamma$-ray all-sky monitoring data.
One of the co-authors of this paper led this effort,
found its $\gamma$-ray variability,
and reported it
via The Astronomer's Telegram (ATel) \citep{tanaka2017_IC170922a}.
This was further followed by
multi-wavelength follow-up observations
all over the wavelengths (Figure~\ref{fig:fig_timeline}; \cite{icecube2018_science1}).

Timeline from the IceCube alert \citep{kopper2017icecubegcn}
to the first Fermi ATel report \citep{tanaka2017_IC170922a} is summarized in
Figure~\ref{fig:fig_timeline}.
After the IceCube alert,
some observational reports
with monitoring and follow-up data
are distributed via GCN and ATel.
In summary,
no possible related objects to the IceCube neutrino were mentioned in any of the reports before \citet{tanaka2017_IC170922a}.
At the event time,
no significant $\gamma$-ray (INTEGRAL SPI-ACS in \cite{savchenko2017}; HAWC in \cite{martinez2017}; H.E.S.S. in \cite{denaurois2017atel}),
no significant neutrino (ANTARES in \cite{dornic2017,dornic2017atel}),
were detected.
In Swift/XRT follow-up observations,
nine sources were detected (including 8 known sources) \citep{keivani2017}
although no special notices were made for TXS~0506+056.

Optical long-term data in this field had been taken in the
ASAS-SN
project since October, 2012.
After the report by \citet{tanaka2017_IC170922a},
TXS~0506+056\ was reported by
\citep{franckowiak2017}
to be at the brightest phase in the ASAS-SN data
and to show a brightening of $\sim0.5$~mag in $V$-band over the last 50 days.

\subsection{Light Curves of TXS~0506+056}\label{sec:sec_lc}

\begin{figure*}
 \begin{center}
  \includegraphics[width=116mm,bb=0 0 566 1077]{test01f.pdf}
 \end{center}
\caption{
Daily variabilities of
optical ($g$; blue, $V$: cyan, $V$ by ASAS-SN: gray,
$R_C$ and $r$: green, $I_C$ and $i$: red, $z$: magenta)
and NIR ($J$: blue, $H$: green, $K_s$: red) fluxes,
optical-NIR colors
from the top to the 3rd panel.
Changes of
the structure functions in time scales of 1.0, 5.0, and 10.0~days
are shown in the 4th panel.
Daily changes of polarization degree ($R_C$: green, $J$: blue),
polarization angle ($R_C$: green, $J$: blue),
Fermi/LAT $\gamma$-ray fluxes in 200-800~MeV and 800~MeV-10~GeV energy bands
are shown from the 5th to the bottom panel.
In the polarization degree, Liverpool/RINGO3 data is plotted in open circles (MJD$\sim58192$).
The neutrino detection time is indicated as gray dashed vertical lines.
Vertical dash-dot gray lines indicate $-150, -100, -50, +50, +100, +150, +200$~days
with respect to the IceCube neutrino detection.
Galactic extinction is not corrected.
}\label{fig:fig_lc}
\end{figure*}

\begin{figure*}
 \begin{center}
  \includegraphics[height=207mm,bb=0 0 407 532]{test04bsraw.pdf}
 \end{center}
\caption{
Optical (top) and NIR (middle) structure functions of TXS~0506+056\
obtained in our observations in thick points and lines.
SFs of the neighboring stars are also shown in thin points and lines
to see how large the SF measurement errors are.
obtained in our observations in thick points and lines.
Power-law and exponential fitting results
for quasars \citep{vandenberk2004} are shown
in blue ($g$), green ($r$), and red ($i$) thick solid and dashed lines in $\Delta{t}>10^{1}$~days, respectively.
Power-law fitting results for FSRQs and
BL~Lacs \citep{bauer2009blazar} are also shown
in green dotted-dashed and dotted lines, respectively.
Flat parts in short time scales are dominated by measurement errors and
real SF amplitudes are thought to be lower than these lines.
Optical ASAS-SN structure function as a function of period (every 30~days)
is also shown in the bottom panel.
}\label{fig:fig_sf}
\end{figure*}

Time variabilities of the observables are shown in Figure~\ref{fig:fig_lc},
including those of
optical magnitudes,
NIR magnitudes,
optical and NIR colors,
optical and NIR polarization degrees and angles, and
$\gamma$-ray fluxes in the low (100-800~MeV) and high (800~MeV-10~GeV) energy bands.
In addition to our own data (\S\ref{sec:sec_optimg}),
optical $V$-band data from
ASAS-SN (see \S\ref{sec:sec_asassndata}),
RINGO3 optical polarimetry data (see \ref{sec:sec_optnirpol}),
and {\it Fermi} $\gamma$-ray data are also used.

The $\gamma$-ray fluxes are taken from the Fermi All-sky Variability
Analysis (FAVA; \cite{abdollahi2017})\footnote{https://fermi.gsfc.nasa.gov/ssc/data/access/lat/FAVA/} website.
Aperture photometry fluxes in the low (100-800~MeV) and high (800~MeV-30~GeV) energy bands
are used in this paper.
We understand that the FAVA light curves are preliminary as described on the FAVA website but
the temporal behavior is similar to those shown in \citet{icecube2018_science1} and
the FAVA data are
good enough for our purpose.

\subsection{Variability of TXS~0506+056}\label{sec:sec_vari}

We investigate
optical and NIR
variabilities in three time scales
including
daily (\S\ref{sec:sec_longtermvari}),
intranight (\S\ref{sec:sec_intranightvari}),
and second-scale (\S\ref{sec:sec_secondvari}) ones.
We construct
structure functions (SFs)
of optical and NIR variability of TXS~0506+056,
which is, in general, defined to be an ensemble variability of an object or a specific set of objects,
to quantify its time variability.
We compare SFs of TXS~0506+056\ with
those of other AGN in the literature.
A caveat and problem are pointed out in \citet{emmanoulopoulos2010} when an SF is used
for studying a time scale of blazar variability,
however, in this paper, we discuss only the variability amplitudes at given time scales,
not the time scale itself.

We adopt an usual concept on
a definition of the structure function (SF) $V(\Delta{t})$ as below,
\begin{eqnarray}
	V(\Delta{t})&=&\sqrt{<|\Delta{m}|^2>-<\sigma^2_{S/N}>}
\end{eqnarray}
where $\Delta{t}$ is a time interval between observations, and
$\Delta{m}$ is a magnitude difference between different observations
and $\sigma_{\rm{S/N}}$ is measurement errors in magnitude unit.
We calculate the SFs
for TXS~0506+056\ and neighboring stars with similar brightness
to TXS~0506+056.
These two SFs
are compared to see any significant flux variability of TXS~0506+056.
We estimate median and $1\sigma$ confidence level
with a bootstrap method.

The ASAS-SN data are useful for evaluating
variability behavior before and around the alert while
our own data enable us to investigate short time scale variabilities after the alert.

Obtained SFs
of TXS~0506+056\ are shown in thick points and lines
in Figure~\ref{fig:fig_sf}.
SFs of the neighboring stars are also shown
in thin faint points and lines in Figure~\ref{fig:fig_sf}.
The SFs of TXS~0506+056 are
compared with those of
blazars (BL~Lacs and FSRQs) studied in the Palomar-QUEST survey \citep{bauer2009blazar},
SDSS quasars \citep{vandenberk2004},
and CTA102 \citep{bachev2017}.
In \citet{bauer2009blazar},
94 BL~Lacs, 278 FSRQs, and 4 marginally classified (BL~Lacs or FSRQs) objects in total are observed
with 48-inch Samuel Oschin Schmidt Telescope and
QUEST2 Large Area Camera ($4.6^{\circ}\times3.6^{\circ}$ field-of-view; \cite{baltay2007})
over 3.5~years to construct the rest-frame SFs
of them as shown in Figure~\ref{fig:fig_sf}.
In \citet{vandenberk2004}, $>25,000$ quasars are observed in the 5 optical bands ($ugriz$)
and rest-frame (both in time scales and wavelengths) SFs
are derived.
These two SFs are derived with different formulations but quantitatively equivalent with each other.

\subsubsection{Long-Term Variability}\label{sec:sec_longtermvari}
Daily or longer time scale variability is summarized in this subsection.
First, optical and NIR fluxes are significantly variable as shown in Figure~\ref{fig:fig_lc}.
The peak-to-peak amplitude reaches up to 1~mag.
We also see some fluctuations in a time scale of a few days.
In addition, optical or NIR fluxes
are not at their peaks at the neutrino detection as well as in $\gamma$-ray \citep{icecube2018_science1}.
The overall trend around the neutrino detection indicates that TXS~0506+056\ is
in a declining phase.
About 180 days after the neutrino detection,
optical and NIR brightness increase again up to a brighter level
than that around the neutrino detection,
however, no neutrino detection is not reported in this period.

The optical SF amplitudes of TXS~0506+056\
are comparable with or slightly larger than
those of AGN in the previous studies
in time scales of $\Delta{t}>10^{1}$~days.
TXS~0506+056\ is more variable by a factor of $\sim2$ than the SDSS quasars
on average \citep{vandenberk2004}.
The SFs of TXS~0506+056\ in our measurements are larger than or comparable to
those of FSRQs and BL~Lacs \citep{bauer2009blazar} at their face values.
Note that the SFs of AGN in general have 0.1~dex or larger scatter (\cite{vandenberk2004}; \cite{bauer2009blazar}),
which make the envelope of the distribution is overlapped
between those of FSRQs and BL~Lacs
and our measurements.

For NIR variability,
TXS~0506+056\ in $J$-band in time scales from a few days to a few tens days,
and $H$ and $K_s$-band in a time scale of a few tens days
are significantly variable than neighboring stars.
The NIR variabilities of TXS~0506+056\
are
as large as optical ones in these time scales.
NIR variabilities in shorter time scales than a few days
are comparable to those of neighboring stars and are
almost at the limit of the measurement errors.
Typical
NIR variabilities of blazars in the literature are
$\sim0.1$~mag (\cite{sandrinelli2014}; \cite{li2018}),
which is comparable to our measurement for TXS~0506+056\
and this variability behavior of TXS~0506+056\ is not special.

We also calculated SFs
using the ASAS-SN $V$-band data
in each 30~days
from MJD=57793 (225~days before the neutrino detection)
to MJD=58243 (225~days after the neutrino detection),
covering the IceCube neutrino detection time
on MJD$\sim58018.87$.
Variability amplitudes in three time scales (1, 5, and 10~days) are
derived from the SFs and
shown in the 4th panel of Figure~\ref{fig:fig_lc}.
As a whole, there are no special periods
when significantly larger variability is detected than other periods.
It is not
clear but
the SFs marginally
indicates that
long-term (10~days) variability around $\sim60$~days before the neutrino detection is the largest
and larger than that in the neutrino detection period
with significance of $2.1\sigma$.
On the other hand, the short-term (1~day) variabilities are constant in time.
This might indicate the neutrino emission would be related to the 10-days-timescale variability in this epoch.
The hard $\gamma$-ray fluxes are also highly variable around this epoch (Figure~\ref{fig:fig_lc}).

Correlation between
the hard $\gamma$-ray fluxes (800~MeV-10~GeV)
and optical brightness is investigated
as seen in
Figure~\ref{fig:fig_optmag2bgamma}.
Optical magnitudes (brightness) are negatively (positively)
correlated with the $\gamma$-ray fluxes with
Spearman rank correlation coefficients of $-0.471$ and
$p$-values of $<0.03$.
This indicates that the correlation is significant and TXS~0506+056\ is brighter in optical
in brighter $\gamma$-ray phases, which is consistent with general trend of ISPs
or all types of blazars \citep{itoh2016,jermak2016}.

\begin{figure}[!htbp]
 \begin{center}
  \includegraphics[width=83.5mm,bb=0 0 396 432]{test02b4b.pdf}
 \end{center}
\caption{
$g$ or $V$-band magnitude
as a function of
$\gamma$-ray flux.
Colors of the points indicate
observing epochs (MJD) as shown in the color bar.
}\label{fig:fig_optmag2bgamma}
\end{figure}

\subsubsection{Intranight Variability}\label{sec:sec_intranightvari}
During 22 nights,
we contiguously took imaging data of TXS~0506+056\
for a few tens minutes or longer
with the 0.4-m Kyoto University telescope,
MITSuME,
Kanata/HONIR,
and IRSF/SIRIUS.
With these datasets, intranight variability can be investigated.

Magnified views of the light curves
on these densely observed nights are shown in
the 1st (optical) and 3rd (NIR) rows of Figure~\ref{fig:fig_lc_m}.
In this figure, photometry is done for each frame.
For a comparison, ensemble relative light curves of nearby stars
with similar brightness on average are also shown
in the 2nd (optical) and 4th (NIR) rows of the figure.
As seen in the daily light curves in Figure~\ref{fig:fig_lc},
brightness change from night to night is easily seen.
R.m.s. of the brightness is also shown in the figure.

\begin{figure*}
 \begin{flushleft}
  \includegraphics[width=221.9mm,bb=0 0 814 539]{test01m3.pdf}
 \end{flushleft}
\caption{
Light curves of TXS~0506+056\ based on photometry for each frame.
From top to bottom,
optical light curves of TXS~0506+056,
optical light curves of nearby similarly bright stars,
NIR light curves of TXS~0506+056,
NIR light curves of nearby similarly bright stars,
and optical and optical-NIR colors of TXS~0506+056\ are shown.
Dates in MJD for the data are shown at the top part in each column.
Colors used for the optical photometry are
blue for $g$,
cyan for $V$,
green for $R_C$ and $r$, and
red for $I_C$ and $i$.
Colors used for the NIR photometry are
blue for $J$,
green for $H$, and
red for $K_s$.
Standard deviations in magnitude unit within a time window are shown at the left-bottom part of the panels.
Optical PS1 photometry in $gri$-bands and NIR 2MASS photometry in $JHK_s$-bands are shown in dashed lines.
A thick gray line in the top second panel shows the length of 10~minutes.
Panels with no data are hatched in gray.
Large scatters ($\gtrsim0.05$~mag) are partly attributed to bad seeing.
}\label{fig:fig_lc_m}
\end{figure*}

Observed intranight variabilities are as small as $0.03-0.11$~mag, depending
on the filters and epochs (the data quality)
and no significant intranight variabilities are detected in our dataset.
The large scatters ($\gtrsim0.05$~mag) seen in some panels of the figure are partly attributed to bad seeing.
In \citet{sagar2004},
a few to 10 percent amplitude ($\lesssim0.1$~mag, $14.1$\% at maximum)
intranight $R_C$-band variability are seen for 11~blazars (6 BL~Lacs and
5 radio core-dominated quasars).
Typical observation duration in a night in \citet{sagar2004} is 6.5~hours
and the total number of the observed nights is 47.
They find that a duty cycle of intranight optical variability is $\sim60$\% for BL~Lacs.
Similarly, \citet{paliya2017} also monitored 17~blazars for $137$~hours
in $R_C$-band
and also obtained
a high duty cycle of $\sim40$\% for blazars.
Compared with these studies, the total duration (sum)
of our observations shown in Figure~\ref{fig:fig_lc_m}
is $\sim18$~hours
and much shorter by factor of $\sim17$ \citep{sagar2004} and $\sim8$ \citep{paliya2017}.
Then, it is difficult to say any consistency of our non-detection with
these two observations.

\subsubsection{Second-Scale Variability}\label{sec:sec_secondvari}

The fastest variability detected for blazars so far is a few minutes
in optical \citep{sasada2008}
and
$\gamma$-ray
\citep{albert2007,aharonian2007,vovk2013,vovk2015,ackermann2016},
which corresponds to a comparable size of a black hole,
which indicates a small emitting region in a relativistic jet.
A mass of a black hole hosting a blazar is expected
to be as massive as $\sim10^{9}$~M$_\odot$ in general \citep{castignani2013},
and the black hole mass of TXS~0506+056\ is estimated to be $\sim3\times10^{8}$~M$_\odot$ \citep{padovani2019}
by assuming the typical host galaxy luminosity of blazars \citep{paiano2018}
and black hole mass and bulge mass relation \citep{mclure2002}.
Detection of variability in a shorter time scale would put
a tight constraint
on the size of an emitting region
with a usually assumed Doppler boosting factor.
If second scale variability was detected, this might be attributed to
an apparent change in viewing angle to a bent relativistic jet \citep{raiteri2017}.

Second-scale variability is
investigated with the 2~fps
Tomo-e~Gozen\ data.
Four sets of 180-second (360~frames) exposure are obtained.
Detection of TXS~0506+056\ is sometimes marginal
and we use only photometric data of signal-to-noise ratios of $>10\sigma$.
As a result, the photometric data
especially
in the last (fourth) exposure are partly removed in our analysis and shown in Figure~\ref{fig:fig_lc_tomoe}.
The resultant SF is shown
in the leftmost part of the top panel of Figure~\ref{fig:fig_sf}.
In our dataset,
no significant rapid variability in a second time scale is detected.

\begin{figure}
 \begin{center}
  \includegraphics[width=95.4mm,bb=0 0 634 396]{test01tg.pdf}
 \end{center}
\caption{
Optical light curves
of TXS~0506+056\ (red) and
a nearby star (blue)
taken with Tomo-e~Gozen.
}\label{fig:fig_lc_tomoe}
\end{figure}

\subsection{Optical and NIR colors of TXS~0506+056}\label{sec:sec_sed}

Temporal changes of
optical, optical-NIR, NIR-NIR colors
are shown
in
Figure~\ref{fig:fig_lc} (entire light curves),
Figure~\ref{fig:fig_lc_m} (intranight light curves),
and Figure~\ref{fig:fig_optmag2a} (correlation between magnitudes and colors).

The range of the optical colors of
$g-r$ and $g-R$
is 0.2-0.5~mag.
A bluer-when-brighter
trend is seen in the color-magnitude diagrams
(Figure~\ref{fig:fig_optmag2a}).
This is consistent with an idea that TXS~0506+056\  is a BL~Lac-type blazar
with little contribution of its accretion disk
to optical-NIR emission \citep{bonning2012}.
In Figure~\ref{fig:fig_optmag2a},
bluish data points indicating the data obtained
around MJD=58200
are offset by $\sim0.4$~mag from the most crowded (redish) data region.
These bluish data also follow a tighter blue-when-brighter trend.
In summary, data points in different epochs follow different
color-magnitude relations.
This behavior is also observed for other blazars,
for example, OJ~287, and
indicates different activity states
between the different loci in the color-magnitude diagram
in different epochs \citep{bonning2012}.
These bluish points are data taken after the $\gamma$-ray flare
in March 2018 \citep{ojha2018} while the redish points are taken after
the IceCube neutrino detection \citep{icecube2018_science1}.
In these two epochs, increased $\gamma$-ray fluxes are detected with Fermi
but the optical and NIR color-magnitude relations are different from each other.
In general,
the brighter locations of the bluish points at a given color may be attributed
to
high-energy electron injection into jet-emitting region
or
an emergence of much brighter accretion disk than usual
possibly due to an accretion state change.
Bluer-when-brighter trends of BL~Lacs are sometimes attributed to a presumption that
the objects are in high state \citep{zhang2015}.
TXS~0506+056\ may
show this trend
in any state
considering
that almost featureless power-law continua are always observed and
that the equivalent widths (EWs) of the emission lines of TXS~0506+056\
(EW$_{\rm{[OIII]}}=0.17$\AA\ at most; \cite{paiano2018}) are small
compared to previously measured EWs of blazars
(although many of these measurements are upper limits) \citep{landt2004}.

The range of $r-J$ or $R_C-J$ colors of TXS~0506+056\ is
from 0.9 to 1.2 (Figure~\ref{fig:fig_optmag2a}),
roughly corresponding to 1.8 to 2.1 in the Vega system.
The $r-J$ and $R_C-J$
color changes are as small as $\sim0.3$~mag.
Although the previous studies examined different colors
(e.g., $V-J$ colors in \cite{ikejiri2011}),
these are typical for ISP blazars.
The bluer-when-brighter trend is also seen.
Note that no data points are shown in the panel
of Figure~\ref{fig:fig_optmag2a}
after the $\gamma$-ray flare in March 2018 reported by \citet{ojha2018}.

Intranight changes in optical colors are shown in the bottom panels of Figure~\ref{fig:fig_lc_m}.
The $g-r$ or $g-R_C$ colors are almost constant, about 0.4~mag all over the nights.
NIR colors of $J-H$ and $H-K_s$ are also almost constant in time, 0.3-0.4~mag.
Note that no significant intranight variability in any band is
detected (\S\ref{sec:sec_intranightvari}).

\begin{figure*}
 \begin{center}
  \includegraphics[width=164.5mm,bb=0 0 616 309]{test02a2.pdf}
 \end{center}
\caption{
Optical magnitudes
($r$ or $R_C$ in the left panel and
$g$ or $V$ in the right panel)
as functions of
optical-NIR color (left panel; $g-J$ or $V-J$)
and optical color (right panel; $g-r$ or $g-R_C$).
Different symbols indicate data taken with different telescopes/instruments
in the same way as in Figure~\ref{fig:fig_lc}.
}\label{fig:fig_optmag2a}
\end{figure*}

\subsection{Polarization of TXS~0506+056}\label{sec:sec_optnirpol_result}

\subsubsection{Temporal Changes of Polarization}\label{sec:sec_optnirpol_result_lc}

Our measurement results of polarization
in $R_C$ and $J$-bands
are summarized in Table~\ref{tab:tab_pol2},
including the optical polarimetric measurement taken with
RINGO3 on the 2-m Liverpool telescope \citep{steele2018}.
The temporal changes of the polarization degrees and angles are shown
in Figure~\ref{fig:fig_lc}.

\begin{table*}[!htbp]
  \tbl{Polarization measurements of TXS~0506+056\ with Kanata/HONIR.}{
  \begin{tabular}{ccrrrrl
      \hline
      Date (UT) & MJD & \multicolumn{2}{c}{Polarization Degree} & \multicolumn{2}{c}{Polarization Angle} & Note\\
                &     & $R_C$ & $J$ & $R_C$ & $J$ & \\
      \hline
      \input{honir_pol_table_20200508a}
	\hline
    \end{tabular}}
    \label{tab:tab_pol2}
\begin{tabnote}
(*) is calculated by assuming the observations is done at UT=0 hours.
\end{tabnote}
\end{table*}

For the first 5~data points taken within 1.5~months
after the alert (defined as ``first epoch''),
the polarization degrees are as small as 2--8\%\
as partly reported in \citet{yamanaka2017_IC170922a}.
About 6~months after the neutrino detection (defined as ``second epoch''),
\citet{ojha2018} reported that
a flare of the highest daily averaged $\gamma$-ray flux for TXS~0506+056\
was detected with {\it Fermi} on March 13, 2018.
Soon after that,
\citet{steele2018} carried out
polarimetric observation on the night of March~14, 2018
with Liverpool/RINGO3 and
found that optical polarization degree increases up to $\sim14$\%
at wavelengths roughly corresponding to
the $R_C$ and $I_C$-bands.
As shown in Figure~\ref{fig:fig_lc} and Table~\ref{tab:tab_pol2},
our subsequent observations with Kanata/HONIR
(12 polarization measurements for 23~days starting 8~days after the RINGO3 observation)
indicate that
optical polarization degree again decreases down to
1.5\%\ on March 22, 2018,
which is even lower than the observed level in the first epoch.
After this decrease, the polarization degrees gradually
increase up to $\sim20$\% during about two weeks
and decreases again down to 7.2\%, which is as low as those in the first epoch.
Throughout this period, the $J$-band polarization exhibits a
similar time variation behavior to that of
$R_C$-band.

In ISPs,
typical polarization degrees and their temporal change
are
$\sim30$~percent and
$\sim20$~percent or less, respectively, among the samples observed by
\citet{ikejiri2011}, \citet{itoh2016}, and \citet{jermak2016}.
The polarization degrees and its temporal change observed in TXS~0506+056\  is comparable
to or less than them, being consistent with those objects.

The polarization angles are roughly
constant at $\sim20^{\circ}\pm10^{\circ}$
in both $R_C$ and $J$-bands
in first epoch.
In the second epoch, our HONIR measurements indicate
that the polarization angles are 20--90~deg and significantly different from
those in the first epoch.
Following the polarization degree increase from 1.5\%\ to $\sim20$\%,
the polarization angles first
change from 84~deg to 22~deg
and back to 30-60~deg.
The position angles in $J$-band also show a similar behavior
to $R_C$-band.

\subsubsection{Correlations between Optical Brightness and Polarization}\label{sec:sec_optnirpol_result_corr}

The top panel of Figure~\ref{fig:fig_optmag2bpol} shows
a relation between the polarization degrees and
optical $r$ or $R_C$-band magnitudes.
Although
optical brightness and polarization degrees look roughly
coincident with each other as a whole (Figure~\ref{fig:fig_lc}),
almost no correlations are seen
during the rapid changes in the second epoch
as indicated by bluish points in Figure~\ref{fig:fig_optmag2bpol}.
This is discussed later in this subsection.

\begin{figure}[!htbp]
 \begin{center}
  \includegraphics[width=75.5mm,bb=0 0 396 734]{test02b4a3.pdf}
 \end{center}
\caption{
$r$ or $R_C$-band magnitude
and
$R_C$-band polarization angle
as a function of
$R_C$-band polarization degree.
All the data points are connected in time sequence.
Different symbols indicate data taken with different telescopes/instruments
in the same way as in Figure~\ref{fig:fig_lc}.
}\label{fig:fig_optmag2bpol}
\end{figure}

\begin{figure}[!htbp]
 \begin{center}
  \includegraphics[width=83.5mm,bb=0 0 382 508]{test01f2.pdf}
 \end{center}
\caption{
A magnified view of
temporal changes of optical magnitudes,
polarization degrees, and
polarization angles
around the $\gamma$-ray flare reported by \citet{ojha2018}.
The symbols and colors are the same as those in Figure~\ref{fig:fig_lc}.
$r$,$R_C$-band in green,
$i$,$I_C$-band in red,
and
$J$-band in blue.
}\label{fig:fig_lc2ndepoch}
\end{figure}

Polarization changes in optical wavelengths associated $\gamma$-ray flares
are reported in previous literature \citep{sasada2008,abdo2010_3c279}.
In the second epoch, the polarization change of TXS~0506+056\
seems associated to the $\gamma$-ray flare \citep{ojha2018}
although such a correlated behavior is not clear in the first epoch
because the polarization measurements are not done right after (or before) the neutrino detection,
the measurements are sparse,
and the number of the measurements is small.
In \citet{itoh2013},
no brightness flare was observed during their first polarization flare
for the famous blazar, CTA~102 (an FSRQ at $z=1.037$),
indicating that the polarization degree does not necessarily correlate with brightness.

The lack of strong correlations between optical brightness and polarization degree in our data for TXS~0506+056\
is similarly observed in previous literature for other blazars \citep{ikejiri2011,jermak2016},
although
a significant
positive correlation
between ``amplitudes'' of flux and polarization degree
is detected
for two blazars (AO~0235+164 and PKS~1510-089)
in 10-day time scale \citep{sasada2011}.
These observed weak correlations
could be partly due to ignorance of a possible time lag between
temporal changes of fluxes and polarization degrees \citep{uemura2017}.
Figure~\ref{fig:fig_lc2ndepoch} is a magnified view of Figure~\ref{fig:fig_lc} around the second epoch.
The peak of polarization degrees is around MJD$\sim$58211--58217~days
while that of optical brightness is around MJD$\sim$58210--58211~days.
This indicates that the optical brightness change proceeds the polarization degree change,
which is
the opposite sense to that observed for a BL~Lac, PKS~1749+096, reported by \citet{uemura2017}.
This lag partly makes the correlation worse in Figure~\ref{fig:fig_optmag2bpol}.
A positive correlation would be seen
if the data point in MJD=58216.5 (the right bottom point
in Figure~\ref{fig:fig_optmag2bpol}) is ignored.

The polarization angles also change in time by $\sim70$~deg in the second epoch
although the change is not so drastically large.
Figure~\ref{fig:fig_lc2ndepoch} indicates that
the polarization angles decrease as the polarization degrees increase
in the period of MJD$=58198-58207$.
The polarization degrees still increase after that, however,
the polarization angles do not show a systematic decrease.
These make the correlation between the polarization angles and degrees poor
as shown in the bottom panel of Figure~\ref{fig:fig_optmag2bpol}.
Changes of polarization angles are observed for many blazars \citep{itoh2016,hovatta2016}.
The observed change of the polarization angles for TXS~0506+056\ is not so large, which is
sometimes explained by
a curved structure of a relativistic jet
\citep{abdo2010_3c279,sorcia2014}.

\subsection{Spectra of TXS~0506+056}\label{sec:sec_spec}

All of our three new spectra of TXS~0506+056\ are shown in Figure~\ref{fig:fig_spec}.
The Subaru/FOCAS spectra in the two different setups shown in \citet{icecube2018_science1}
are also plotted as references.
All these spectra were taken roughly in similar epochs,
$\sim1$~week to $\sim1.5$~months
after the neutrino detection, and
before the GTC/OSIRIS spectrum
which conclusively determines the redshift of TXS~0506+056\ was taken \citep{paiano2018}.
In the epochs of our observations,
TXS~0506+056\ is slightly brighter ($g=14.8-15.1$)
than that in the GTC observation ($g=15.4$; \cite{paiano2018}),
which makes line detections more difficult.
All of our three spectra show basically featureless continua and
no significant emission or absorption lines are detected,
except for the weak emission line in the Subaru/FOCAS spectrum \citep{icecube2018_science1}.
These are quantitatively
consistent with the spectrum of \citet{paiano2018}.
Any significant changes between the different epochs are not detected.

\begin{figure}
 \begin{center}
  \includegraphics[width=90mm,bb=0 0 408 624]{spec_all_fig01.pdf}
 \end{center}
\caption{
Optical spectra of TXS~0506+056\ taken with
Nayuta/MALLS (red),
Kanata/HOWPol (blue),
Gemini-South/GMOS (green),
and
Subaru/FOCAS (black and gray).
The Subaru/FOCAS spectra are the same as those shown in \citet{icecube2018_science1}.
There are no changes in spectral features.
Noisy wavelength ranges are omitted.
All the spiky features seen in the spectra are all atmospheric.
}\label{fig:fig_spec}
\end{figure}

\section{Summary}\label{sec:sec_summary}

We first made optical and NIR imaging observations to search for a candidate of
the neutrino source of IceCube-170922A.
We found that TXS~0506+056\ was rapidly fading in NIR
in a day scale.
Motivated by this discovery of the rapid NIR variability,
the $\gamma$-ray flare was discovered in the Fermi/LAT monitoring data.
We conducted monitoring observations of TXS~0506+056\ with
Akeno/MITSuME,
Kiso/KWFC,
Kiso/Tomo-e~Gozen,
Kyoto 0.4-m,
Kanata/HONIR, and
Subaru/HSC.
Polarization imaging data were also taken with Kanata/HONIR at 17~epochs.
We also took optical spectra \sout{four times} five times with
Kanata/HOWPol, Nayuta/MALLS, Gemini-N/GMOS, and Subaru/FOCAS.

Combined with ASAS-SN optical monitoring data and Fermi/LAT $\gamma$-ray data in addition to our data described above, we find
\begin{itemize}
    \item daily variability is significant and the amplitude is as large as 1.0~mag,
    \item no significant intranight variability is detected,
    \item no special optical/NIR variability behavior over the observed period,
    but a very marginal sign of a larger variability in a time scale of $\sim10$ days about 2~months before the neutrino detection,
    \item a weak correlation between optical/NIR and $\gamma$-ray fluxes
    \item large changes in polarization degrees and angles about 180~days after the neutrino detection,
    \item a weak or no correlation between polarization degree and optical fluxes while some correlated behavior can be seen
    in a part of the March 2018 ($\sim180$~days after the neutrino detection) data,
    \item no significant optical spectral changes over the three months after the neutrino detection.
\end{itemize}

In summary, our data do not indicate that
TXS~0506+056\ is a special blazar (BL~Lac) among other blazars in terms of intranight, daily,
monthly optical/NIR variability, and optical polarization.
The neutrino detection is also not a special timing.

For future Icecube neutrino events, to further examine possible relations
between high-energy neutrinos and blazar
as well as blazar variability, more complete blazar catalogs like a recently
developed one, BROS \citep{itoh2018bros},
are required.
In addition, in optical and NIR wavelengths,
routine high-cadence imaging and polarization monitoring of blazars are desired.
Transient survey
observations are done
in several wide-field surveys such as Zwicky Transient Facility (ZTF; \cite{graham2019}) and
Tomo-e Gozen \citep{sako2018}.
Polarization monitoring is a more expensive observation program but is expected
to make unique science outputs
(e.g., \cite{sakimoto2012,itoh2016}).

\begin{ack}

We thank an anonymous referee for his/her comments to improve the manuscript.
Observations with the Kanata, Nayuta, MITSuME, IRSF, and Kiso Schmidt telescopes were supported by the Optical and Near-infrared Astronomy Inter-University Cooperation Program and the Grants-in-Aid of the Ministry of Education, Science, Culture, and Sport
JP23740143,
JP25800103,
JP16H02158,
JP18H04575,
JP15H02075,
JP16H06341,
JP18H05223,
JP17K14253,
JP17H04830,
JP26800103,
JP19H00693,
JP17H06363,
JP17H06362,
and
JP24103003.

This work was partly supported by the joint research program of the
Institute for Cosmic Ray Research (ICRR).
This work was also based in part on data collected at Subaru Telescope,
which is operated by the National Astronomical Observatory of Japan.
Based on observations obtained at the international Gemini Observatory, a program of NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. on behalf of the Gemini Observatory partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigaci\'{o}n y Desarrollo (Chile), Ministerio de Ciencia, Tecnolog\'{i}a e Innovaci\'{o}n (Argentina), Minist\'{e}rio da Ci\^{e}ncia, Tecnologia, Inova\c{c}\~{o}es e Comunica\c{c}\~{o}es (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea).

The IRSF project is a collaboration between Nagoya University and the
SAAO supported by the Grants-in-Aid for Scientific Research on
Priority Areas (A) (No.~10147207 and No.~10147214) and Optical \&
Near-Infrared Astronomy Inter-University Cooperation Program, from the
Ministry of Education, Culture, Sports, Science and Technology (MEXT)
of Japan and the National Research Foundation (NRF) of South Africa.

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

The data is partly processed using the Gemini IRAF package.
The authors acknowledge Dr.~N.~Matsunaga for a software to conduct comparisons between our image coordinates and the catalogue values.

\end{ack}