\section{Introduction}
Single impurity in a fermionic or bosonic environment has been investigated intensively over the past several decades in solid state physics,
and many important concepts such as Anderson's orthogonality catastrophe, Kondo effect, polarons have been introduced.
As well as bulk solid-state systems, ultracold atomic systems and nanostructured systems such as quantum dots coupled to some leads have also provided intriguing platforms to study the single impurity problems thanks to recent progress in experimental techniques.
In nanostructured systems, the Kondo effect has been observed in transport properties \cite{Goldhaber-Gordon1998,Goldhaber-Gordon1998-2,Cronenwett1998,Schmid1998}. More recently, polaronic systems have been experimentally realized in ultracold atomic systems for a strongly imbalanced mixture of two species of atoms \cite{Schirotzek2009,Nascimbene2009,Will2011,Koschorreck2012,Kohstall2012},
and the realization of Kondo physics with ultracold atomic gases has been proposed theoretically \cite{Bauer2013,Nishida2013,Nakagawa2015}.
Because ultracold atomic systems and nanostructured systems have high controllability,  more detailed studies of equilibrium properties and also nonequilibrium properties of impurity systems, which are difficult to study in bulk solid-state systems, have been performed both experimentally \cite{Catani2012,Fukuhara2013,Fukuhara2013a,Hild2014,Kogan2004,Latta2011} and theoretically
\cite{Nordlander1999,Lobaskin2005,Anders2005,Anders2006,Hettler1995,Ng1996,Schiller1996,Goldin1997,Nordlander1998,Lopez1998,Goldin1999,Kaminski1999,Kaminski2000,Lopez2001,Heyl2010,Knap2012,Knap2013}.

In particular, real-time dynamics of the standard impurity models, such as the Anderson model and the Kondo model, has been studied intensively in nonequilibrium conditions; dynamics under periodic driving \cite{Hettler1995,Ng1996,Schiller1996,Goldin1997,Nordlander1998,Lopez1998,Goldin1999,Kaminski1999,Kaminski2000,Lopez2001,Heyl2010} as well as post-quench dynamics (i.e. dynamics after a sudden change of the system parameters) \cite{Nordlander1999,Lobaskin2005,Anders2005,Anders2006} have been explored. These studies have focused on nanostructured systems in which some leads and quantum dots are coupled with each other, and therefore have concentrated on the transport quantities such as the current between two leads through a dot. We will investigate below the time-dependent local quantities such as the impurity spin polarization and the spin density in a fermionic bath. This is partly motivated by recent experimental advances in cold atomic systems which enable us to observe the time-dependent local densities, for example, by using the single-site-resolved imaging technique in ultracold atomic systems \cite{Fukuhara2013,Fukuhara2013a,Hild2014}. Also, we note that these quantities have not been investigated in detail in the previous studies focused on the nanostructured systems. We will elucidate nontrivial time-dependent phenomena, especially a novel resonance phenomenon, emerging in the long-time limit of a periodically driven Kondo system.

More importantly, the detailed analysis of such real-time dynamics of the Kondo system provides an opportunity to address a fundamental problem of current interest in condensed matter physics, that is, how the periodically driven quantum systems behave in the long time limit. This problem has attracted much attention recently because novel properties of matter emerge by applying time-periodic external field to quantum systems \cite{Oka2009,Lindner2011,Kitagawa2011,Aidelsburger2013,Miyake2013}. Actually, recent studies have revealed intriguing properties of periodically-driven {\it isolated} systems where the Magnus expansion \cite{Magnus1954,Blanes2009} (a series expansion in the inverse of the external driving frequency) is not guaranteed to converge \cite{DAlessio2013,Lazarides2014b,Ponte2015,Mori2015,Abanin2015a,Kuwahara2016,Mori2016}. On the other hand, for {\it open} quantum systems, there are many issues to be explored. One of the problems is whether periodically-driven open systems can be described by some equilibrium systems \cite{Kohler1997,Breuer2000,Kohn2001,Hone2009,Ketzmerick2010,Iadecola2013,Langemeyer2014,Iadecola2015b,Shirai2015,Liu2015,Iadecola2015a,Shirai2016}. Recently, it has been shown that the Gibbs distribution of the Floquet states can emerge only when certain special conditions are fulfilled \cite{Shirai2015,Liu2015,Iadecola2015a,Shirai2016}. In general, an asymptotic state of the periodically driven open systems is not described by the Floquet-Gibbs state, so that it is important to study such open quantum systems in more detail to figure out in what conditions equilibrium-state properties emerge

To address the above-mentioned problems, in this paper, we study the steady state of a periodically driven anisotropic Kondo impurity which is coupled to a fermionic bath with an infinite band width. The time dependence we consider is for the local magnetic field which couples to the impurity spin and also for the in-plane interaction strength.
We obtain the exact expression of time evolution for arbitrary time dependence of them by using the exact solution known as the Toulouse limit \cite{Toulouse1969}.
First, when the local magnetic field is periodically switched on and off and the in-plane interaction is time independent,
we find a non-monotonic behavior of the time average of the impurity spin polarization, which is regarded as a resonance phenomenon, as a function of the intensity of the local magnetic field.
Second, by taking the high frequency limit of the external driving, we show that the system recovers time translational invariance and the properties of the system are described by a mixture of the zero-temperature state and the infinite-temperature state in general.
In certain cases, the system has the properties either at zero or infinite temperature, and thereby can be described by the equilibrium state in these special cases.

This paper is organized as follows.
In Sec. \ref{Sec_ModelMethod}, we summarize the Kondo model at the Toulouse limit and the method to obtain the exact analytical expression of the time evolution in the Heisenberg picture when the time periodic local magnetic field is applied to the system.
In Sec. \ref{Sec_RectangularCase}, we investigate a specific case where the local magnetic field is periodically switched on and off.
 In Sec. \ref{Sec_TimeDepInt}, we extend the calculations in Sec. \ref{Sec_ModelMethod} to the case where the in-plane interaction strength is also time dependent. This model enables us to address the question posed in the introduction: in what conditions equilibrium-state properties emerge for an open system. For this purpose, we consider the case where the local magnetic field and the in-plane interaction are both time dependent, and analyze the properties of the high frequency limit of the external driving. A brief summary of our results is presented in Sec. \ref{Sec_Conclusion}.

\section{Model and Method \label{Sec_ModelMethod}}
We consider a periodically-driven impurity spin coupled to a fermionic bath via the anisotropic Kondo exchange interaction,
\begin{equation}
	\begin{aligned}
		H(t) = \sum _{\sigma} & \int dx : \psi ^{\dagger} _{\sigma} (x) \big( -iu \partial _x \big) \psi _{\sigma} (x) : \\
							  & +\sum _{i=x,y,z} J_i S^i s^i (0) -h(t)S^z.
	\end{aligned}
\end{equation}
The operator $  \psi _{\sigma} (x) \ (\sigma = \uparrow , \downarrow) $  annihilates a fermion in the bath and the colons $: \cdots :$ denote the normal ordering.
$s^i (x) = \sum _{ss^{\prime}} :\psi ^{\dagger} _s (x) \sigma ^i _{ss^{\prime}} \psi _{s^{\prime}}(x): \ (i=x,y,z)$ are the spin density of the fermionic bath, where $\sigma ^i _{ss^{\prime}}$ are the Pauli matrices.
$ S^i \ (i=x,y,z)$ are the impurity spin operators whose spin is 1/2 and the couplings $J_i \ (i=x,y,z)$ have anisotropy: $J_x =J_y = J_{\perp},\ J_z=J_{\|}$.

We consider a time-dependent local magnetic field with the periodicity $h(t + \tau)=h(t)$ and assume that the initial state is in the ground state of the Kondo model without a local magnetic field.
Because the scattering term is point like and only fermions whose angular momentum is zero (s-wave) are scattered, the effective model is reduced to a one-dimensional system \cite{Affleck1991a}.
To simplify the problem, we take the strength of $z$-component coupling as $ J_{\|}/\pi u = \sqrt{2} (\sqrt{2}-1) $ (Toulouse limit \cite{Toulouse1969,Zarand2000}).
It is known that the exact results can be obtained at this point while keeping most of essential properties of the Kondo model intact.
At this point, we can transform the Kondo model to a non-interacting resonant level model by using the bosonization/refermionization methods \cite{Giamarchi2003}:
\begin{equation}
	\begin{aligned}
		H(t) = & \int dx : \tilde{\psi}^{\dagger} (x) \big( -iu \partial _x \big) \tilde{\psi} (x): \\
				& + \frac{J_{\perp}}{\sqrt{2 \pi \alpha}} \big( \tilde{\psi} ^{\dagger} (0) \tilde{c}_d + \text{h.c.} \big) -h(t) \big( \tilde{c}_d ^{\dagger} \tilde{c}_d -\frac{1}{2} \big) ,
	\end{aligned}\label{Eq_ResonantLevel}
\end{equation}
where the newly introduced operators $ \tilde{\psi} (x) $ and $ \tilde{c}_d $ are fermionic.
The impurity spin polarization is represented by the occupation number of the resonant level, $S^z = \tilde{c}_d ^{\dagger} \tilde{c}_d -1/2 $.
The spin density of the fermionic bath is $s^z (x) = \sqrt{2}:\tilde{\psi}^{\dagger} (x) \tilde{\psi} (x): \ (x\neq 0)$.
The Kondo temperature, which provides a typical energy scale characterizing this quantum impurity system, can be determined by the impurity spin contribution to the specific heat as $T_K = \pi w J_{\perp} ^2/4 \pi \alpha u$, where $ w=0.4128 $ is known as the Wilson number \cite{Lobaskin2005}.

We calculate the time evolution of the system in the Heisenberg picture.
The time evolution operator of this system is $ U(t) = \mathcal{T} \exp \big[ -i \int _0 ^t du H(u) \big] $, where $ \mathcal{T} $ is the time-ordering operator. Because Hamiltonians of different time do not commute with each other, it is difficult to obtain the explicit expression for $ U(t) $.
Thus, we use the following strategy to obtain the time evolution of operators.
First, we divide the time-periodic local magnetic field $ h(t) $ into discrete $M$ time-steps:
\begin{equation}
	h(t) = \sum _{N=0} ^{\infty} \sum _{n=1} ^M \theta \big( t-(N-(n-1)/M) \tau \big) \theta \big( (N+n/M)\tau -t \big) h^{(n)}.
\end{equation}
When $ h(t) $ is divided into $M$ time-steps, the time evolution of operators in each period reads
\begin{equation}
	\tilde{c}_l (\tau) = e^{i H^{(1)} \tau /M} \cdots e^{i H^{(M)} \tau /M} \tilde{c}_l e^{-i H^{(M)} \tau /M} \cdots e^{-i H^{(1)} \tau /M},
\end{equation}
where $ l=k,d $ and $ \tilde{c}_k $ is a Fourier coefficient of $ \tilde{\psi}(x) $ $( \tilde{\psi} (x) = 1/\sqrt{L} \sum _k e^{ikx} \tilde{c}_k )$.
$ H^{(n)} $ is a time-independent Kondo model with a local static magnetic field $ h^{(n)} $.
As the Hamiltonian $ H^{(n)} $ is quadratic, the time evolution by $ H^{(n)} $ is denoted as,
\begin{equation}
	\begin{aligned}
		c_l (t) &= e^{iH^{(n)}t}c_l e^{-iH^{(n)}t} \\
				&= \sum _{l^{\prime}} G^{(n)} _{ll^{\prime}} (t) c_{l^{\prime}} , \ (l=k,d) ,
	\end{aligned}
\end{equation}
where $ G^{(n)} _{ll^{\prime}} (t) $ is the transition matrix from $ l^{\prime} $ to $ l $.
Thus, the time evolution in a period is given by
\begin{equation}
	\begin{aligned}
		& \tilde{c}_l (\tau) = \sum _{l^{\prime}} M_{ll^{\prime}} \tilde{c}_{l^{\prime}} \\
		& M_{ll^{\prime}} = \sum _{ \{ l_i \} } G^{(M)} _{ll_1} \big(\tau /M \big) G^{(M-1)} _{l_1 l_2} \big(\tau /M \big) \cdots G^{(1)} _{l_{M-1} l^{\prime}} \big(\tau /M \big) .
	\end{aligned}
\end{equation}
The explicit expression of $ M_{ll^{\prime}} $ is shown in the Appendix \ref{appendix_a}.
By taking the continuum limit, the time evolution of the operators $ \tilde{\psi}(x) $ and $ \tilde{c}_d $ is obtained in the long-time limit,
\begin{equation}
	\begin{aligned}
		& \begin{aligned}
			\tilde{\psi} (x,t) = & \tilde{\psi}(x-ut) \\
								 & -i \sum _k \frac{J_{\perp}}{\sqrt{2 \pi \alpha u^2}} \theta (x) T_{dk} (t-x/u) \tilde{c} _k \\
		  \end{aligned} \\
		& \tilde{c}_{d}(t) = \sum _k T_{dk} (t) \tilde{c} _k \\
		& \begin{aligned}
			T_{dk} (t) & = \frac{J_{\perp}}{\sqrt{2 \pi \alpha L}} \frac{-ie^{-i\epsilon _k t}}{1-e^{i\epsilon _k \tau} M_{dd}} \\
						& \times \int _0 ^{\tau} ds \exp \Big[ i \int _{\tau -s} ^{\tau} h(s^{\prime} +t) ds^{\prime} - \Delta \tau +i \epsilon _k s \Big] \\
		  \end{aligned} \\
		& M_{dd} = \exp \Big[ i \int _0 ^{\tau} h(t) dt - \Delta \tau \Big] ,\ \Delta = \frac{J_{\perp} ^2}{4\pi \alpha u}.
	\end{aligned}
	\label{Eq_TransitionMatrix}
\end{equation}
These expressions are applicable for an arbitrary time-periodic function $h(t)$ and contains the whole information about the dynamics of the system including the fermionic bath. We note that similar exact analytical expressions for time-dependent current between the two leads through a nanostructure are obtained by using the Keldysh Green function method \cite{Wingreen1993,Schiller1996}.

The relation between the spin density of the fermionic bath and the time derivative of the impurity spin is also obtained:
\begin{equation}
	\braket{s^z(x,t)}=\begin{cases}
						-\frac{\sqrt{2}}{u}\braket{\frac{d}{dt} S^z (t-x/u)} & (x>0) \\[5pt]
						0													 & (x<0)
					  \end{cases}
					  \label{Eq_SpinRelation}
\end{equation}
This relation means that a fermion whose spin is antiparallel to the impurity spin is scattered and propagates with the Fermi velocity $u$
because the fermionic bath and the impurity spin are coupled by antiferromagnetic interaction and the fermionic bath has a linear dispersion.

\section{Time-Dependent Properties and emergent Resonance Phenomenon  \label{Sec_RectangularCase}}
\subsection{Time evolution of observables \label{SubSec_RectangularTimeEvolution}}
\begin{figure*}
	\centering
	\includegraphics[width=150mm]{RectangularTimeEvolution.eps}
	\caption{Time evolution of the impurity spin polarization (a), (c) and the spin density of the fermionic bath (b), (d) in a period for several choices of local magnetic field intensity $h$
			($h=10\Delta$, blue solid line; $h=2.5\Delta$, red dashed line; $h=0.625\Delta$, yellow dashed dotted line).
			(a) and (b) are the figures for $ \tau = 10/\Delta $, while (c) and (d) are for $ \tau = 0.3125/\Delta $.\label{Fig_RectangularTimeEvolution}}
\end{figure*}
Here, we discuss a specific case: the local magnetic field is periodically switched on and off. Note that the following analysis and discussion can be straightforwardly applied to more generic periodically driven systems.
The time dependence of $ h(t) $ we consider is
\begin{equation}
	h(t) = \sum _{N=0} ^{\infty} \theta \big( t-N\tau \big) \theta \big( N\tau + \tau /2 -t \big) h .
\end{equation}
Decompose time into periodic intervals:	$t=N\tau +s , \ N \in \mathbb{N} , \ 0<s<\tau$.
Then, the time evolution of the impurity spin polarization in the steady state
in the long-time limit
is
\begin{equation}
	\braket{S^z (t)} \xrightarrow{N\rightarrow \infty}
		\begin{cases}
			\begin{aligned}
				A_1 + A_2 & e^{-(T_K /\pi w)s} +\text{Re} \big[ A_3 (s) e^{ihs}] \\
						& (0<s<\tau /2)
			\end{aligned} \\
			\begin{aligned}
				A_1 ^{\prime} +A_2 ^{\prime} & e^{-(T_K /\pi w)(s-\tau /2)} \\
															& (\tau /2 <s<\tau).
			\end{aligned}
		\end{cases}
\end{equation}
The prefactors $ A_i ^{(\prime)} (s) $ can be obtained by the momentum integral in the fermionic-bath sector.
The second term is the relaxation term due to the switching of the field intensity.
The third term in $ 0<s<\tau /2 $ represents the oscillation caused by the local magnetic field, which comes from the interference between the impurity spin and the spins in the fermionic bath. Recall that
the spin density of the fermionic bath is represented by the time derivative of the impurity spin polarization (see Eq. (\ref{Eq_SpinRelation})).

The time evolution of the impurity spin polarization $ \braket{S^z (t)} $ and the spin density of the fermionic bath $ \braket{s^z(x,t)} $ in a period are shown in Fig. \ref{Fig_RectangularTimeEvolution}.
Because the impurity spin polarization $ \braket{S^z (t)} $ is time periodic in the steady state, the spin density of fermionic bath $ \braket{s^z(x,t)} $ is spatially periodic as well as time periodic, as seen from Eq. (\ref{Eq_SpinRelation}).
We thus show only the data at the point $x=Mu\tau \ (M\in \mathbb{N})$.
In the first half of the period the local magnetic field is switched on (i.e. $ h(t) = h $), while in the second half it is off (i.e. $ h(t) = 0 $).
In the long driving period (see the figures of $ \tau = 10 /\Delta $), the impurity spin polarization oscillates with a frequency proportional to the field intensity $h$, and eventually relaxes to its equilibrium value in the time scale proportional to $1/T_K$.
Wave packets whose width is of the order of the inverse of the Kondo temperature are observed in the spin density of the fermionic bath since the spin density of the fermionic bath is given by the time derivative of the impurity spin polarization.
These wave packets propagate with the Fermi velocity $u$ without decaying due to the linear dispersion.
In the short driving period (see the figures of $ \tau = 0.3125 /\Delta$), the system cannot keep up with the temporal change of the local magnetic field and thus approaches constant values.
Note that a large cusp is observed at a half of the period in the spin density of the fermionic bath, because it is given by the time derivative of the impurity spin polarization and thus amplified by the factor $\Omega=2\pi /\tau$.

\subsection{Emergent resonance in the time averaged observables}
In the limit of fast driving, the impurity spin polarization $ \braket{S^z (t)} $ approaches a temporally constant value (see Fig. \ref{Fig_RectangularTimeEvolution}). Under this condition, we investigate how the time average of the impurity spin polarization, $ \overline{S^z} $, depends on the local magnetic field intensity $h$, and elucidate their intriguing resonance phenomenon. Here we define the time average $ \overline{S^z} $ and its variance $ \Delta S^z $:
\begin{equation}
	\begin{aligned}
		& \overline{S^z} = \lim _{N \rightarrow \infty} \int _0 ^{\tau} \frac{ds}{\tau} \braket{S^z (N\tau +s)} \\
		& \Delta S^z = \lim _{N \rightarrow \infty} \int _0 ^{\tau} \frac{ds}{\tau} \Big| \braket{S^z (N\tau +s)} -\overline{S^z} \Big| ,
	\end{aligned}
\end{equation}
where $ \lim _{N\rightarrow \infty} $ is taken to extract the steady state.

The obtained results are shown in Fig. \ref{Fig_RectangularAverage} as a function of the field intensity $h$ for a choice of the driving period, $\tau=0.3125\Delta$. Naively, we expect that $ \overline{S^z} $ increases monotonically and saturates to $1/2$ with increasing $h$. Interestingly, however, an emergent oscillation is observed as a function of $h$ and its period is determined by the driving period of the local magnetic field. When $ h \gg T_K $, the time average of the impurity spin polarization $ \overline{S^z}(h) $ becomes periodic as a function of $h$
\begin{equation}
	\overline{S^z} (h+4\pi /\tau) = \overline{S^z}(h) \ (h \gg T_K) ,
\end{equation}
where $ \tau $ is the driving period of the local magnetic field.

To see the essence of this behavior, let us define the ``power" $P$ and the ``work done in a period" $ W $ for this system, based on an analogy to the power in classical mechanics,
\begin{equation}
	\begin{aligned}
		&P(t)=\frac{dS^z}{dt} \frac{dh(t)}{dt} \\
		&W= \lim _{N \rightarrow \infty} \int _0 ^{\tau} P(N \tau +s) ds .
	\end{aligned}
\end{equation}
It is expected that $ dS^z /dt $ corresponds to the velocity of a particle $v=dx(t)/dt$, and $ dh(t)/dt $ corresponds to the force $F$ in classical mechanics, that is, the following correspondence to classical mechanics is expected:
\begin{equation}
	\frac{dS^z (t)}{dt} \frac{dh(t)}{dt} \leftrightarrow \bm{v} \cdot \bm{F} = \frac{d}{dt} \text{(Kinetic energy)}
\end{equation}
The quantity $W$ is shown as a function of $h$  in Fig. \ref{Fig_RectangularAverage}. The emergent oscillation in the time average of the impurity spin polarization as a function of $h$ can be attributed to a resonance phenomenon between the oscillation of the impurity spin polarization (see Fig. \ref{Fig_RectangularTimeEvolution}) and switching of the local magnetic field.
The impurity spin polarization oscillates with the period of $ 2 \pi / h $ in the first half of the period, and the local magnetic field is switched off just at the half of the period.
Analogy to an oscillator in classical mechanics would suggest that if the time derivative of the impurity spin polarization is negative at the half of the period where the local magnetic field is switched off, the amplitude toward the negative direction would become larger, and therefore the time average would be smaller.
Note that the ``work done in a period" $W$ oscillates with the same period but an opposite phase, whereas the variance of the impurity spin polarization  $ \overline{S^z} $ is small but oscillates with the same phase as the ``work". These behaviors are indeed consistent with our interpretation. The resonance condition deduced by the analogy mentioned above is
\begin{equation}
	\frac{2 \pi n}{h} = \frac{\tau}{2} \Leftrightarrow h=\frac{4 \pi n}{\tau}  \ (n \in \mathbb{N}) ,
\end{equation}
which reproduces the period observed in Fig. \ref{Fig_RectangularAverage}.
\begin{figure}[t]
	\centering
	\includegraphics[width=80mm]{RectangularAverage.eps}
	\caption{Time-averaged quantities as a function of the local magnetic field intensity $h$: Time-averaged impurity spin polarization $\overline{S^z}$ (blue solid line), Variance $\Delta S^z$ (red dashed line), Work $W$ (yellow dashed dotted line).\label{Fig_RectangularAverage}}
\end{figure}

\section{Some novel properties driven by time-dependent interactions\label{Sec_TimeDepInt}}
When the in-plane interaction $ J_{\perp} $ is also time-periodic, the system shows quite different behaviors from those without the in-plane interaction driving. The method in the Sec. \ref{Sec_ModelMethod} can be easily extended to the case where both the local magnetic field $h(t)$ and the in-plane interaction $J_{\perp}(t)$ are time-periodic functions.
The time evolution of the operators in this condition is shown in the Appendix \ref{appendix_a}, where the relation between the spin density of the fermionic bath and the time derivative of the impurity spin polarization (Eq. (\ref{Eq_SpinRelation})) also holds in this case.

In particular, we investigate below our driven system in the high frequency limit of the external driving and elucidate the characteristic long-time behavior of the periodically driven system.
Note that the Magnus expansion cannot be applied to this system due to the infinite band width and the Kondo model can be regarded as an open quantum system where a localized spin couples to a reservoir.
Therefore, it is instructive to study the properties of this system from the viewpoints of periodically driven open quantum systems as well as the divergence of the Magnus expansion. This is indeed related to the question whether the driven system can be described by some equilibrium states.

\subsection{Dynamics of the observables}
Here, we consider a specific case where the time dependence of the local magnetic field $h(t)$ and the in-plane interaction parameter $J_{\perp}(t)$ is given b
\begin{equation}
	h(t)=h\sin (2\Omega t) , \ J_{\perp}(t)=J_{\perp}\sin \Big( \Omega t +\frac{\pi}{4} \Big) .
\end{equation}
Note that the time dependence of the local magnetic field $h(t)$ equals that of the square of the in-plane interaction parameter $ J_{\perp} ^2 (t) $, which has the essence of dynamics when both the local magnetic field and the in-plane interaction parameter are time dependent.
The dynamics of the impurity spin polarization in the steady state is
\begin{equation}
	\begin{aligned}
		& \braket{S^z(N\tau +s)} \xrightarrow{N \rightarrow \infty} e^{x \cos (2\Omega s)} A(s) -\frac{1}{2} \\
		& A(t) = \int _{-\infty} ^0 \frac{d \omega}{2 \pi} \Big| \sum _{p \in \mathbb{Z}} \frac{J_p (z)-J_{p+1} (z)}{\omega -(2p+1)\Omega /\overline{\Delta} +i} e^{2ip(\Omega t+\pi /4)} \Big| ^2 \\
		& \overline{\Delta} = \frac{1}{4\pi \alpha u} \int _0 ^{\tau} \frac{dt}{\tau} J_{\perp} ^2 (t) , \ x=\frac{\overline{\Delta}}{\Omega} , \ z=\frac{h+i\overline{\Delta}}{2\Omega} ,
	\end{aligned}
\end{equation}
where $J_n (z)$ is the integer Bessel function.
Because $A(t+\tau /2) = A(t)$, the impurity spin polarization is time periodic whose period is half of the driving period $\tau$: $ \braket{S^z(t+\tau /2)}=\braket{S^z(t)} $.
This is because the in-plane interaction parameter acts as an effective hopping parameter, and thereby the square of the in-plane interaction parameter appears in the observables.

\begin{figure*}[t]
	\centering
	\includegraphics[width=150mm]{InteractionTimeEvolution.eps}
	\caption{Time evolution of the impurity spin polarization (a), (c) and the spin density of the fermionic bath (b), (d) in a half of the driving period for several choices of local magnetic field $h$
			($h=10\overline{\Delta}$, blue solid line; $h=2.5\overline{\Delta}$, red dashed line; $h=0.625\overline{\Delta}$, yellow dashed dotted line).
			(a) and (b) are the figures for $ \tau = 20/\overline{\Delta} $, while (c) and (d) are for $ \tau = 0.625/\overline{\Delta} $.\label{Fig_InteractionSinusoidalDynamics}}
\end{figure*}

In Fig. \ref{Fig_InteractionSinusoidalDynamics}, we show the time evolution of the impurity spin polarization $ \braket{S^z (t)} $ and the spin density of a fermionic bath $ \braket{s^z(x,t)} \ (x=Mu\tau ,\ M\in \mathbb{N})$ in a half of the driving period.
Let us define the effective Kondo temperature by $ T_K ^{\eff} = \pi w \overline{\Delta} $ (see Appendix \ref{appendix_b} for the meaning of this definition).
When the driving period is much longer than the time scale of the inverse of the effective Kondo temperature (see the figures of $ \tau = 20 / \overline{\Delta} $), the system follows a temporal change in the external driving in the time scale of $1/T_K ^{\eff}$. Therefore, the behavior of time evolution is similar to  the case without in-plane interaction driving (see Sec. \ref{Sec_RectangularCase}).
When the driving period is much smaller than the inverse of the effective Kondo temperature (see the figures of $ \tau = 0.625 / \overline{\Delta} $), these values also approach temporally constant values as we have seen for the case without in-plane interaction driving.

It should be noted, however, that the time average of the impurity spin polarization does not vanish even though the local magnetic field does not have the static component (i.e. zeroth mode of Fourier coefficients).
This is because the square of the in-plane interaction parameter becomes small while the local magnetic field is negative.
Namely, the impurity spin is effectively decoupled from the fermionic bath when the local magnetic field is negative.
In the high frequency limit $(\overline{\Delta}/\Omega \rightarrow 0)$, the time average of impurity spin polarization is $\overline{S^z} = \sum _{n \in \mathbb{N}} J_n(h/2\Omega)J_{n+1}(h/2\Omega)$.
This behavior is consistent with the results by Heyl and Kehrein \cite{Heyl2010} in the linear response regime when the coupling term is periodically switched on and off.

\subsection{High frequency limit: analysis of asymptotic steady states}

We now demonstrate that intriguing steady states appear in the high frequency limit. The following analysis is valid for any type of time dependence of $h(t)$ and $J_{\perp}(t)$. Note that the time evolution of the system is determined solely by the function $T_{dk}(t)$ in the steady state (see Eq. (\ref{Eq_TimeEvolution})), so it is sufficient to study the properties of $T_{dk}(t)$.
In the language of the resonant level model, $T_{dk}(t)$ denotes a transition matrix from a state with wave-number $k$ in the fermionic bath to the resonant level. For arbitrary $h(t)$ and $J_{\perp} (t)$, each element of the transition matrix $T_{dk}(t)$ is represented as,
\begin{equation}
	\begin{aligned}
		&T_{dk}(t) = F(t) \sum _{n \in \mathbb{Z}} T_{dk} ^{(n)} (t) \\
		&T_{dk} ^{(n)} (t) = \frac{\tilde{J}^{(n)}}{\sqrt{2\pi \alpha L}} \frac{e^{-i(\epsilon _k -n\Omega)t}}{\epsilon _k +\overline{h}-n\Omega +i\overline{\Delta}},
	\end{aligned}
\end{equation}
where $\overline{f}=\int _0 ^{\tau} f(t) dt /\tau$ is a time average in a period and the coefficients $\tilde{J}^{(n)}$ and $F(t)$ are determined from the Fourier coefficients of the external driving $h(t)$ and $J_{\perp}(t)$.
The transition matrix without driving $T_{dk} ^{\text{bare}}(t)$ in the long time limit is
\begin{equation}
	T_{dk} ^{\text{bare}} (t) \xrightarrow{t \rightarrow \infty} \frac{J_{\perp}}{\sqrt{2\pi \alpha L}} \frac{e^{-i\epsilon _k t}}{\epsilon _k +h+i\Delta}.
\end{equation}
Thus, $T_{dk} ^{(n)} (t)$ can be interpreted as a "$n$th photon assisted transition matrix".

By taking the high frequency limit $(\overline{\Delta}/\Omega, \ \overline{h}/\Omega \rightarrow 0)$, the impurity-spin polarization $\braket{S^z (t)}$ and the impurity-spin time correlation function, $\braket{S^z (t)S^z(t^{\prime})}$, recover time translational invariance.
We thus end up with the analysic formuae in this limit,
\begin{widetext}
	\begin{equation}
		\begin{aligned}
			&\braket{S_z (t)} \xrightarrow{t \rightarrow \infty, \ \overline{\Delta} /\Omega,\ \overline{h}/\Omega \rightarrow 0} \frac{1}{\pi} \Bigg[ \frac{\Delta ^{(0)}}{\overline{\Delta}} \int _{-\infty} ^{0} + \sum _{n=-1} ^{-\infty} \frac{\Delta ^{(n)}}{\overline{\Delta}} \int _{-\infty} ^{\infty} \Bigg] \frac{d \omega}{(\omega +\overline{h}/{\overline{\Delta}})^2 +1} -\frac{1}{2} \\
			&\begin{aligned}
				&\braket{S_z (t) S_z (t^{\prime})} - \braket{S_z (t)} \braket{S_z (t^{\prime})} \\
				& \begin{aligned}
						\xrightarrow{t,\ t^{\prime} \rightarrow \infty, \ \overline{\Delta} /\Omega,\ \overline{h}/\Omega \rightarrow 0} \Bigg( \Bigg[ \frac{\Delta ^{(0)}}{\overline{\Delta}}& \int _{-\infty} ^{0} +\sum _{n=-1} ^{-\infty} \frac{\Delta ^{(n)}}{\overline{\Delta}} \int _{-\infty} ^{\infty} \Bigg] \frac{d \omega}{\pi} \frac{e^{i\omega \overline{\Delta} (t-t^{\prime})}}{(\omega + \overline{h} / \overline{\Delta})^2+1} \Bigg) \\
																				 & \times \Bigg( \Bigg[ \frac{\Delta ^{(0)}}{\overline{\Delta}} \int _{-\infty} ^{0} +\sum _{n=1} ^{\infty} \frac{\Delta ^{(n)}}{\overline{\Delta}} \int _{-\infty} ^{\infty} \Bigg] \frac{d \omega}{\pi} \frac{e^{i\omega \overline{\Delta} (t-t^{\prime})}}{(\omega - \overline{h} / \overline{\Delta})^2+1} \Bigg) ,
				  \end{aligned}
			 \end{aligned}
		\end{aligned}
		\label{Eq_FastDrivingLimit}
	\end{equation}
\end{widetext}
where $\Delta ^{(n)} = |\tilde{J}^{(n)}|^2 / 4\pi \alpha u$.
In the equilibrium Kondo model with a local magnetic field $h$ at Toulouse limit and temperature $T$,
the impurity-spin polarization $ \braket{S^z}_{eq} $ and the impurity-spin time correlation function \cite{Guinea1985,Leggett1978} $ \braket{S^z (t) S^z }_{eq} $ are expressed as follows:
\begin{equation}
	\begin{aligned}
		&\begin{aligned}
			\braket{S_z}_{eq} &= \int _{-\infty} ^{\infty} \frac{d\omega}{\pi} \frac{f(\omega \Delta /T)}{(\omega +h/\Delta)^2 +1} -\frac{1}{2}
		 \end{aligned} \\
		& \begin{aligned}
			&\braket{S_z(t)S_z}_{eq} -\big( \braket{S^z}_{eq} \big) ^2 \\
								& = \Big( \int _{-\infty} ^{\infty} \frac{d\omega}{\pi} \frac{f(\omega \Delta /T)e^{i\omega \Delta t}}{(\omega +h/\Delta)^2 +1} \Big) \Big( \int _{-\infty} ^{\infty} \frac{d\omega}{\pi} \frac{f(\omega \Delta /T)e^{i\omega \Delta t}}{(\omega -h/\Delta)^2 +1} \Big) ,
		\end{aligned}
	\end{aligned}
	\label{Eq_EquilibriumSystem}
\end{equation}
where $f(\epsilon)=1/(e^{\epsilon}+1)$ is the Fermi distribution function.

Comparison between Eq. (\ref{Eq_FastDrivingLimit}) and Eq. (\ref{Eq_EquilibriumSystem}) gives instructive implications for the asymptotic states.
Because $f(\omega \Delta /T) \xrightarrow{T\rightarrow 0} \theta (-\omega)$ and $f(\omega \Delta /T) \xrightarrow{T \rightarrow \infty} \frac{1}{2}$,
the integrals from $-\infty$ to $0$ in Eq. (\ref{Eq_FastDrivingLimit}) have the properties of the zero temperature state (the time correlation function shows a power law decay $\braket{S^z (t) S^z (t^{\prime})} \propto (t-t^{\prime})^{-2} \ (t\gg t^{\prime})$), while the integrals from $-\infty$ to $\infty$ in Eq. (\ref{Eq_FastDrivingLimit}) have the properties of the infinite temperature state (the time correlation function shows an exponential decay $\braket{S^z(t)S^z(t^{\prime})} \propto e^{-2\overline{\Delta}|t-t^{\prime}|}$).
Thus, the properties of the system are described by {\it a mixture of the zero-temperature state and those of the infinite-temperature state} in general.
This means that the reservoir we consider here cannot prevent the localized spin from heating up completely.
The properties of the zero temperature state come from the bare transitions (i.e. transition processes without photons) while the properties of the infinite temperature come from  the photon-assisted transition processes.
If $\Delta ^{(0)} = 0$ and $\Delta ^{(n)}=\Delta ^{(-n)}$ $(\Delta ^{(n\neq 0)} = 0)$ the properties of the system are described by those at zero (infinite) temperature, respectively. In these special cases, the asymptotic state can be described by the equilibrium state.

Two concrete examples exhibiting either zero- or infinite-temperature properties for asymptotic states are in order here. First, consider the case where the local magnetic field is sinusoidal $h(t)=h\sin (\Omega t)$ and the in-plane interaction is time independent, $\tilde{J}^{(n)} = J_n (h/\Omega)$. Then, the impurity-spin polarization and the impurity-spin time correlation function are exactly reduced to those of the infinite temperature state for the intensities $h$ which satisfy the zeros of the zeroth Bessel function.
This is easily confirmed from Eq. (\ref{Eq_FastDrivingLimit}), where we can put $\Delta^{(0)}=0$. The second example is the $h(t)/\Omega \rightarrow 0$ case: if we take the limit of $h(t)/\Omega \rightarrow 0$, then $\tilde{J}^{(n)} = J_{\perp} ^{(n)}$, where $J_{\perp} ^{(n)}$ is the $n$th Fourier coefficient of $J_{\perp}(t)$. In this case, the steady state cannot be described by the time averaged Hamiltonian even though the frequency of the external driving is set to infinity.
While the impurity-spin polarization and the impurity-spin time correlation function are exactly equivalent to those of the infinite-temperature state when $J_{\perp} ^{(0)} = 0 \ (\Leftrightarrow \overline{J_{\perp}}=0)$,
they are exactly equivalent to those of the zero-temperature state when $J_{\perp} ^{(n)} = 0 \ (\forall n \neq 0) \ (\Leftrightarrow J_{\perp} \text{ is time independent})$.
The other cases cannot be described effectively by any equilibrium systems due to the fluctuation dissipation theorem \cite{Heyl2010}.

To close this section, we make a brief comment on the nonequilibrium Fermi distribution function.
Note that the above asymptotic behaviors originate from the duplication of the Fermi edge (or equivalently duplication of the impurity level) \cite{Tien1963,Kohler2005}.
In calculating observables, we evaluate $\braket{\tilde{c}^{\dagger} _d (t) \tilde{c}_d (t^{\prime})}$, which can be written down explicitly as,
\begin{equation}
	\begin{aligned}
		&\braket{\tilde{c}^{\dagger} _d (t) \tilde{c}_d (t^{\prime})}=\sum _k f(\beta \epsilon _k) \sum_n \Big( T_{dk} ^{(n)} (t) \Big) ^{\ast} T_{dk} ^{(n)} (t^{\prime})\\
		&\propto \sum _k \sum _n |\tilde{J}^{(n)}|^2 f(\beta (\epsilon _k +n\Omega)) \Big( T_{dk} ^{\text{bare}} (t) \Big) ^{\ast} T_{dk} ^{\text{bare}} (t^{\prime})
	\end{aligned}
\end{equation}
The above observation enables us to define the nonequilibrium distribution function $f_{\text{noneq}} (\beta \epsilon) \propto \sum _n |\tilde{J}^{(n)}|^2 f(\beta (\epsilon +n\Omega))$.
The nonequilibrium distribution function has Fermi edges at $\epsilon = n\Omega \ (n\in \mathbb{N})$ and the weight of each edge is $|\tilde{J}^{(n)}|^2$.
Because the Fermi edge at $\epsilon =0$ can pick up the contribution of the density of states of the resonant level and the others cannot, a transition without photons results in the zero-temperature behavior while the others result in the infinite temperature behavior.

\section{Conclusions\label{Sec_Conclusion}}
We have studied the dynamics of a periodically driven anisotropic Kondo model. By carrying out calculations at the Toulouse limit, the exact analytical expressions of the time evolution for an arbitrary periodic time dependence of the local magnetic field and the in-plane interaction have been obtained.

Focusing on the specific case where the local magnetic field is periodically switched on and off, we have specified the time scale of the dynamics and found the intriguing resonance phenomenon between the switching of the local magnetic field and the oscillation of the impurity spin polarization. This resonance results in a characteristic non-monotonic behavior in the time average of the impurity spin polarization as a function of the local magnetic field.

We have also elucidated the intriguing properties emerging in the steady state in the high frequency limit. In that limit, the impurity-spin polarization and the impurity-spin time correlation function recover the time translational invariance. In particular, we have found that the properties of the system are in general described by a mixture of the zero-temperature and infinite-temperate properties, and for special cases, the system can be described either by zero-temperature or infinite-temperature properties. In the latter special cases, the system has the equilibrium-state properties.
If we regard the Kondo model as an open quantum system, the above ``mixture" of states means that the fermionic bath cannot prevent the localized spin from heating up completely. We expect that this behavior may be common to the open quantum systems where the Hamitonian is given by a bilenear form.  Our simple system gives an example which captures the essential structure of such steady states. Extending the present analysis to more generic cases beyond the bilinear Hamitonian is an interesting issue to be explored in the futrure study.

\begin{acknowledgments}
This work was partly supported by a Grand-in-Aid for Scientific Research on Innovative Areas (JSPS KAKENHI Grant No. 15H05855) and also JSPS KAKENHI (No. 25400366).
\end{acknowledgments}


\section{Introduction}

The spacetimes that are {\em Locally Rotationally Symmetric} (LRS) have been studied in detail and discussed many times in the literature in the cosmological context, i.e. with a fluid matter source  (see for example \cite{Ellis_1967,Ellis_1968,Elst_Ellis_1996} and the references therein). For these spacetimes there exists a continuous isotropy group at each point and hence there is a multiply-transitive isometry group acting on the spacetime manifold. As we know, the isotropies around a point in a spacetime with a fluid can occur as a one-dimensional or three dimensional subgroup of the full group of isometries that leaves the normalised 4-velocity of the matter flow invariant. A three dimensional group of isotropies at each point implies  that the spacetime is isotropic at every point and gives rise to the homogeneous and isotropic Friedmann-Lema\^{i}tre-Robertson-Walker (FLRW) models. While a one dimensional group of isotropies at each point  corresponds to anisotropic and in general spatially inhomogeneous models\footnote{generically with one or two centres where the isotropy group is 3-dimensional}, but includes also some spatially homogeneous (Bianchi and Kantowski-Sachs) models \cite{Mac_Ellis_1969, King_Ellis_1973}. The important property of LRS spacetimes is that they exhibit locally (at each point) a unique preferred spatial direction, covariantly defined, (for example, by a vorticity vector field, a non-vanishing non-gravitational acceleration of the matter, or a density gradient).

LRS spacetimes with a perfect fluid matter source have been completely analysed and classified by Stewart and Ellis using tetrad methods \cite{Ellis_1968}.
 Using a semi-tetrad covariant formalism it was shown that the Einstein field equations can be written as a set of first order equations of geometrical scalars \cite{Elst_Ellis_1996, Clarkson:2007yp}. By analysing the consistency conditions of the field equations, it was rigorously proved that a perfect fluid LRS spacetime cannot have a simultaneous fluid rotation and spatial twist of the preferred spatial direction. Based on this observation, the perfect fluid LRS spacetimes can be divided into three distinct classes. {\em Class I} spacetimes are those where the rotation is non zero but the twist vanishes. This class was shown to be non-expanding, non-distorting and stationary and the solutions generalise the well known G{\"o}del solution. In {\em Class II} spacetimes both the rotation and the twist vanish and these consist of the spherical, hyper-spherical, and plane symmetric (cylindrical) solutions. {\em Class III} spacetimes have no rotation or acceleration but non-zero twist of the preferred spatial direction. These spacetimes are spatially homogeneous.

Though all these classes are of
 interest, and LRS-II solutions have been used extensively to study spherically symmetric astrophysical objects, none of them are suitable for modelling a dynamical rotating star (gravitational collapse of a rotating star, for example). For LRS-I, the rotation is non-zero but the spacetime is stationary, while the other two classes allow dynamical solutions with vanishing rotation. In this study the three key questions are: {\em By relaxing the condition of a perfect fluid, that is by introducing pressure anisotropy and heat flux, is it possible to have dynamical solutions with non-zero rotation and non-zero twist? If yes, can these solutions be physical? What are the local geometrical properties  of such solutions?}

In this paper we investigate in detail the above questions by using the semi-tetrad 1+1+2 covariant formalism \cite{Clarkson:2002jz,Betschart:2004uu,Clarkson:2007yp}. We first establish the existence of such solutions and then find the constraints on the thermodynamic quantities of matter that generate such solutions. We also demonstrate that there exist physically realistic solutions where the matter satisfies physically reasonable energy conditions.

The paper is organised as follows: In the next two sections we describe briefly the basic concepts of local semi-tetrad 1+3 and 1+1+2 covariant formalisms. In the subsequent sections we discuss the properties of LRS spacetimes and the field equations written in terms of the 1+1+2 geometrical variables. In section 5, we then proceed to show the existence of dynamic solutions for imperfect fluids (with pressure anisotropy and heat flux) with non-zero rotation and spatial twist. We also investigate transparently the constraints that the thermodynamic quantities of the matter must satisfy for such solutions to exist. We provide a brief algorithm indicating how to solve the system of field equations with the given initial data. Finally we briefly discuss about how these solutions can be used as a first approximation to spherical symmetry in order to study rotating, inhomogeneous and dynamic astrophysical objects.

Unless otherwise specified, we use natural units ($c=8\pi G=1$)  and $(-,+,+,+)$ signature throughout this paper.
The symbol $\nabla$ represents the usual covariant derivative. The Riemann tensor is defined by
\begin{equation}
R^{a}{}_{bcd}=\Gamma^a{}_{bd,c}-\Gamma^a{}_{bc,d}+ \Gamma^e{}_{bd}\Gamma^a{}_{ce}-\Gamma^e{}_{bc}\Gamma^a{}_{de}\;,
\end{equation}
and the Ricci tensor is obtained by contracting the {\em first} and {\em third} indices
\begin{equation}\label{Ricci}
R_{ab}=g^{cd}R_{cadb}\;.
\end{equation}
The Hilbert--Einstein action in the presence of matter is given by
\begin{equation}
{\cal S}=\frac12\int d^4x \sqrt{-g}\left[R-2\Lambda-2{\cal L}_m \right]\;,
\end{equation}
variation of which gives the Einstein field equations as
\begin{equation}
G_{ab}+\Lambda g_{ab}=T_{ab}\;
\end{equation}
where $G_{ab} := R_{ab} - \frac{1}{2}R g_{ab}$, $R :=R^a_{\,\,\,a}$, and $\Lambda$ is the cosmological constant.

\section{1+3 decomposition of spacetime}

With respect to a timelike congruence, the spacetime can be locally decomposed into time and space parts. One natural way to define such a timelike congruence would be along the matter flow, with the \textit{four-velocity} defined as
\begin{equation}
u^a = \frac{dx^a}{d\tau}, \quad \mbox{with} \quad u^a u_a = -1 ,
\end{equation}
where $\tau$ is the proper time. Given the \textit{four-velocity} $u^a$, we have the unique \textit{projection tensors}
\begin{eqnarray}
U^a{}_b &=& -u^a u_b ,
\\
h^a{}_b &=& g^a{}_b+u^au_b,
\end{eqnarray}
where $h^a{}_b$ is the projection tensor that projects any 4D vector or tensor onto the local 3-space orthogonal to $u^a$. It follows that
\begin{align*}
U^a{}_c U^c{}_b &= -U^a{}_b &  U^a{}_b u^b &= u^a, & U^a{}_a &=1, \\
h_{ab}u^b &= 0 & h^a{}_c h^c{}_b &= h^a{}_b, & h^a{}_a &=3.
\end{align*}

With the choice of this timelike vector, we have two well defined directional derivatives. We have the vector $u^{a}$ which is used to define the \textit{covariant time derivative} (denoted by a dot) for any tensor $ S^{a..b}{}_{c..d}$, given by
\begin{equation}
\dot{S}^{a..b}{}_{c..d}{} = u^{e} \nab_{e} {S}^{a..b}{}_{c..d}
\end{equation}
and we have the tensor $h_{ab}$ which is used to define the fully orthogonally \textit{projected covariant derivative} $D$ for any tensor $
S^{a..b}{}_{c..d} $:
\begin{equation} D_{e}S^{a..b}{}_{c..d}{} = h^a{}_f
h^p{}_c...h^b{}_g h^q{}_d h^r{}_e \nab_{r} {S}^{f..g}{}_{p..q}\;,
\end{equation}
with total projection on all the free indices.
The splitting of the spacetime gives a 3-volume element
\begin{equation}
\ep_{abc}= \eta_{abcd}u^d, \; \mbox{where} \; \ep_{abc} = \ep_{[abc]} \; \mbox{and} \; \ep_{abc}u^c=0.
\end{equation}
Since $\eta_{abcd}$ is the four-dimensional volume element, i.e., $\eta_{abcd}=\sqrt{|\mbox{det}\,g|}\delta^0_{\left[ a
\right. }\delta^1_b\delta^2_c\delta^3_{\left. d \right] }$, we have
\begin{equation}
\eta_{abcd} = 2u_{\left[ a \right.}\ep_{\left. bcd\right]}.
\end{equation}
Since $\eta_{abcd}$ is skew-symmetric, the following contractions hold
\begin{eqnarray}
\ep_{abc} \ep^{def} &=& 3!h^d{}_{\left[ a \right. }h^e{}_bh^f{}_{\left. c \right] }, \\
\ep_{abc} \ep^{dec} &=& 2h^d{}_{\left[ a \right.}h^e{}_{\left. b \right]}, \\
\ep_{abc} \ep^{dbc} &=& 2h^d{}_a, \\
\ep_{abc} \ep^{abc} &=& 3.
\end{eqnarray}

The covariant derivative of $u^a$ can be decomposed as
\begin{equation}
\nabla_a u_b = -u_a A_b + D_a u_b,
\end{equation}
where $D_a$ totally projects derivatives onto the 3-space. $D_a u_b$ can be decomposed into the trace part, the trace-free symmetric part and the trace-free anti-symmetric part, i.e.,
\begin{equation}
\nabla_a u_b = -u_a A_b +\frac13\Theta h_{ab}+\sigma_{ab}+\ep_{abc}\om^{c},
\end{equation}
where $A_b=\dot u_b$ is the acceleration, $\Theta=D_au^a$ represents the expansion of $u_a$,
$\sigma_{ab}=\bra{h^c{}_{\left( a \right.}h^d{}_{\left. b \right)}-\sfr13 h_{ab} h^{cd}}D_cu_d$ is the shear tensor that denotes the distortion and $\om^{c}$ is the vorticity vector denoting the rotation.

The \textit{Weyl curvature tensor} $C_{abcd}$, which gives the locally free gravitational field, is defined by the equation
\begin{equation}
C^{ab}{}_{cd} := R^{ab}{}_{cd}-2g^{\left[ a \right.}{}_{\left[c \right.}R^{\left. b\right]}{}_{\left. d\right]}+\sfr13Rg^{\left[ a \right.}{}_{\left[c \right.}g^{\left. b\right]}{}_{\left. d\right]}.
\end{equation}
Since the Weyl tensor is trace-free on all its indices $\bra{C^c{}_{acb}=0}$, the Ricci tensor $R_{ab}$ is the trace of $R_{abcd}$, and $C_{abcd}$ is the trace-free part. The Weyl tensor can be split relative to $u^a$ into the \textit{electric} and \textit{magnetic Weyl curvature} parts as
\begin{eqnarray}
E_{ab} &=& C_{abcd}u^bu^d  \\
&\Rightarrow&  E^a{}_{a} = 0,~~E_{ab} = E_{\langle ab \rangle},~~E_{ab}u^b = 0,
\end{eqnarray}
and
\begin{eqnarray}
H_{ab} &=& \sfr12\ep_{ade}C^{de}{}_{bc}u^c
\\
&\Rightarrow&  H^a{}_{a} = 0,~~ H_{ab} = H_{\langle ab \rangle},~~ H_{ab}u^b = 0.
\end{eqnarray}
The energy momentum tensor of matter can be decomposed similarly as
\begin{equation}
T_{ab}=\mu u_au_b+q_au_b+q_bu_a+ph_{ab}+\pi_{ab}\;,
\end{equation}
where $p=(1/3 )h^{ab}T_{ab}$ is the isotropic pressure, $\mu=T_{ab}u^au^b$ is the energy density, $q_a=q_{\langle a\rangle}=-h^{c}{}_aT_{cd}u^d$ is the 3-vector that defines the heat flux, and $\pi_{ab}=\pi_{\langle ab\rangle}$ is the anisotropic stress.

\section{1+1+2 decomposition of spacetime}

The 1+1+2  decomposition is a natural extension of 1+3 decomposition, where with respect to a given spatial direction the 3-space is further decomposed, that is we now have another split along a preferred spatial direction. We choose a spacelike vector field $e^a$  such that
\begin{equation}
u^a e_a =0 \quad \mbox{and} \quad e^a e_a=1.
\end{equation}

Then the new projection tensor is given by
\begin{equation} \label{N1}N_a^{~b}\equiv
h_a^{~b}-e_ae^b=g_{a}^{~b}+u_au^b-e_ae^b\, . \end{equation}
This tensor projects vectors onto local 2-spaces, defined as \textit{sheets} (note that these are not subspaces of the 3-space if the twist of $e^a$ is nonzero). Thus
\begin{equation}\label{N2}
e^a N_{ab}=0=u^a N_{ab}, \quad N^a{}_a=2.
\end{equation}
The volume element of this sheet is
\begin{equation} \varepsilon_{ab}\equiv\ep_{abc}e^c = u^d\eta_{dabc}e^c\;\end{equation}
Using the definitions of  $\varepsilon_{ab}$ and $N_{ab}$, we have the following conditions
\begin{eqnarray}
 \varepsilon_{ab}e^b &=&0=\varepsilon_{(ab)}\,,
\\
 \varepsilon_{abc} &=&  e_a \varepsilon_{bc} - e_b  \varepsilon_{ac} + e_c  \varepsilon_{ab}\,,
\\
 \varepsilon_{ab} \varepsilon^{cd} &=& N_a{}^c N_b{}^d - N_s{}^d N_b{}^c ,
\\
 \varepsilon_a{}^c \varepsilon_{bc} &=& N_{ab} ,
\\
 \varepsilon^{ab} \varepsilon_{ab} &=& 2.
\end{eqnarray}
Any 3-vector $\psi^a$ can now be irreducibly split into a scalar, $\Psi$, which
is the vector component parallel to $e^a$, and a vector,
$\Psi^a$ that lies in the sheet as follows:
\begin{eqnarray}
\psi^a&=&\Psie^a+\Psi^{a},~~~\mbox{where}~~~\Psi\equiv \psi_ae^a\
,\nonumber\\&&~~~\mbox{and}~~~\Psi^{a}\equiv N^{ab}\psi_b\equiv \psi^{\bar a}\label{psia},
\end{eqnarray} where the bar over the index denotes projection with $N_{ab}$.
Similarly, the same can be done for any 3-tensor, $\psi_{ab}$,
\begin{equation}
\psi_{ab}=\psi_{\langle
ab\rangle}=\Psi\bra{e_ae_b-\sfr12N_{ab}}+2\Psi_{(a}e_{b)}+\Psi_{{ab}}\
, \label{tensor-decomp}
\end{equation} where
 \begin{eqnarray}
\Psi&\equiv &e^ae^b\psi_{ab}=-N^{ab}\psi_{ab}\ ,\nonumber\\
\Psi_a&\equiv &N_a^{~b}e^c\psi_{bc}=\Psi_{\bar a}\ ,\nonumber\\
\Psi_{ab}&\equiv &
\bra{N_{(a}^{~~c}N_{b)}^{~~d}-\sfr{1}{2}N_{ab}N^{cd}}\psi_{cd}
\equiv \Psi_{\{}%{\lceil ab\}}%{\rfloor}\ . \label{PSTF-TT}
\end{eqnarray}
Apart from the `{\it time}' (dot) derivative, we introduce two new derivatives, which for any tensor $ \psi_{a...b}{}^{c...d}  $:
\begin{eqnarray}
\hat{\psi}_{a..b}{}^{c..d} &\equiv & e^{f}D_{f}\psi_{a..b}{}^{c..d}~,
\\
\delta_f\psi_{a..b}{}^{c..d} &\equiv & N_{a}{}^{p}...N_{b}{}^gN_{h}{}^{c}..
N_{i}{}^{d}N_f{}^jD_j\psi_{p..g}{}^{i..j}\;.
\end{eqnarray}
The derivative along the $e^a$ vector-field in the surfaces orthogonal to $ u^{a}$ is called the hat-derivative, while the derivative projected onto the sheet is called the $\delta$ -derivative. This projection is on every free index.

Using (\ref{psia}) and (\ref{tensor-decomp}), the usual 1+3 kinematical and Weyl quantities can now be split into the irreducible set \[{\cal D}_1 :=  \{\Theta,{\cal A},\Omega,\Sigma,{\cal E},{\cal H},{\cal A}^a,\Sigma^a,{\cal
E}^a,{\cal H}^a,\Sigma_{ab},{\cal E}_{ab},{\cal H}_{ab}\}. \] The 4-acceleration,vorticity and shear split as
\begin{eqnarray}
\dot{u}^a&=&{\cal A} e^a+{\cal A}^a,\\
\omega^a&=&\Omega e^a+\Omega^a,\\
\sigma_{ab}&=&\Sigma\bra{e_ae_b-\sfr{1}{2}N_{ab}}+2\Sigma_{(a}e_{b)}+\Sigma_{ab}.
\end{eqnarray}
For the electric and magnetic Weyl tensors we get
\begin{eqnarray}
E_{ab}&=&{\cal E}\bra{e_ae_b-\sfr{1}{2}N_{ab}}+2{\cal E}_{(a}e_{b)}+{\cal E}_{ab},\\
H_{ab}&=&{\cal H}\bra{e_ae_b-\sfr{1}{2}N_{ab}}+2{\cal H}_{(a}e_{b)}+{\cal
H}_{ab}.
\end{eqnarray}
Similarly the fluid variables, $q^a$ and $\pi_{ab}$, may be split as follows
\begin{eqnarray}
q^a&=&Q e^a+Q^a,\\
\pi_{ab}&=&\Pi\bra{e_ae_b-\sfr{1}{2}N_{ab}}+2\Pi_{(a}e_{b)}+\Pi_{ab}.
\end{eqnarray}

By decomposing the covariant derivative of $e^a$ in the direction orthogonal to $u^a$ into it's irreducible parts, i.e., the spatial derivative of $e^a$, we get
\begin{equation}
 {D}_{a}e_{b} = e_{a}a_{b} + \frac{1}{2}\phi N_{ab} +
\xi\epsilon_{ab} + \zeta_{ab}~,
\end{equation}
where
\begin{eqnarray}
a_{a} &\equiv & e^{c}{\rm D}_{c}e_{a} = \hat{e}_{a}~, \\
\phi &\equiv & \delta_ae^a~, \\  \xi &\equiv & \frac{1}{2}
\epsilon^{ab}\delta_{a}e_{b}~, \\
\zeta_{ab} &\equiv & \delta_{\{}%{\lceil a}e_{b \}}%{\rfloor }~.
\end{eqnarray}
Here, $\phi$ represents the \textit{expansion of the sheet},  $\zeta_{ab}$ is the \textit{shear}, i.e., the distortion of the sheet, $a^{a}$ its \textit{acceleration} and $\xi$ is its (spatial) \textit{vorticity}, i.e., the ``twisting'' or rotation of the sheet.

\section{LRS spacetime field equations}

As discussed in the Introduction, in LRS spacetime there exists a unique, preferred spatial direction at each point and this preferred direction is covariantly defined.  This direction creates a local axis of symmetry, i.e., all observations are identical under rotations about it. In particular, they are the same in all spatial directions that are perpendicular to that direction.
Hence the 1+1+2 decomposition described in the previous section is ideally suited for the study of LRS spacetimes. We can immediately see that if we choose the spacelike unit vector $e^a$ along the preferred spatial direction of the spacetime, then by symmetry all the sheet vectors and tensors vanish identically.
 Thus, all the non-zero 1+1+2 variables are covariantly defined scalars. The geometrical scalar variables that fully describe LRS spacetimes are \[ {\cal D}_2 := \brac{{\cal A}, \Theta,\phi, \xi, \Sigma,\Omega, {\mathcal E}, \H, \mu, p, \Pi, Q }. \] Decomposing the Ricci identities for $u^a$ and $e^a$ and the doubly contracted Bianchi identities, we now get the following field equations for LRS spacetimes:
\smallskip

\textit{Evolution}:
\begin{eqnarray}
   \dot\phi &=& \bra{\sfr23\Theta-\Sigma}\bra{{\cal A}-\sfr12\phi}
+2\xi\Omega+Q\ , \label{phidot}
\\
\dot\xi &=& \bra{\sfr12\Sigma-\sfr13\Theta}\xi+\bra{{\cal A}-\sfr12\phi}\Omega
\nonumber \\ && +\sfr12 {\mathcal H},  \label{xidot}
\\
\dot\Omega &=& {\cal A}\xi+\Omega\bra{\Sigma-\sfr23\Theta}, \label{dotomega}
\\
\dot {\mathcal H} &=& -3\xi{\mathcal E}+\bra{\sfr32\Sigma-\Theta}{\mathcal H}+\Omega Q
\nonumber\\ && +\sfr32\xi\Pi.
\end{eqnarray}
\smallskip

\textit{Propagation}:
\begin{eqnarray}
\hat\phi  &=&-\sfr12\phi^2+\bra{\sfr13\Theta+\Sigma}\bra{\sfr23\Theta-\Sigma}
    \nonumber\\&&+2\xi^2-\sfr23\bra{\mu+\Lambda}
    -{\mathcal E} -\sfr12\Pi,\,\label{hatphinl}
\\
\hat\xi &=&-\phi\xi+\bra{\sfr13\Theta+\Sigma}\Omega , \label{xihat}
\\
\hat\Sigma-\sfr23\hat\Theta&=&-\sfr32\phi\Sigma-2\xi\Omega-Q\
,\label{Sigthetahat}
 \\
  \hat\Omega &=& \bra{{\cal A}-\phi}\Omega, \label{Omegahat}
\\
\hat{\mathcal E}-\sfr13\hat\mu+\sfr12\hat\Pi&=&
    -\sfr32\phi\bra{{\mathcal E}+\sfr12\Pi}+3\Omega{\mathcal H}
 \nonumber\\&&   +\bra{\sfr12\Sigma-\sfr13\Theta}Q , \label{Ehatmupi}
\\
\hat {\mathcal H} &=& -\bra{3{\mathcal E}+\mu+p-\sfr12\Pi}\Omega
\nonumber\\&&-\sfr32\phi {\mathcal H}-Q\xi,
\end{eqnarray}
\smallskip

\textit{Propagation/evolution}:
\begin{eqnarray}
   \hat{\cal A}-\dot\Theta&=&-\bra{{\cal A}+\phi}{\cal A}+\sfr13\Theta^2
    +\sfr32\Sigma^2 \nonumber\\
    &&-2\Omega^2+\sfr12\bra{\mu+3p-2\Lambda}\ ,\label{Raychaudhuri}
\\
    \dot\mu+\hat Q&=&-\Theta\bra{\mu+p}-\bra{\phi+2{\cal A}}Q \nonumber \\
&&- \sfr32\Sigma\Pi,\,
\\    \label{Qhat}
\dot Q+\hat
p+\hat\Pi&=&-\bra{\sfr32\phi+{\cal A}}\Pi-\bra{\sfr43\Theta+\Sigma} Q\nonumber\\
    &&-\bra{\mu+p}{\cal A}\ ,
\end{eqnarray}\\
\begin{eqnarray}
\dot\Sigma-\sfr23\hat{\cal A}
&=&
\sfr13\bra{2{\cal A}-\phi}{\cal A}-\bra{\sfr23\Theta+\sfr12\Sigma}\Sigma\nonumber\\
        &&-\sfr23\Omega^2-{\mathcal E}+\sfr12\Pi\, ,\label{Sigthetadot}
\\
\dot{\mathcal E} +\sfr12\dot\Pi +\sfr13\hat Q&=&
    +\bra{\sfr32\Sigma-\Theta}{\mathcal E}-\sfr12\bra{\mu+p}\Sigma \nonumber \\
  && -\sfr12\bra{\sfr13\Theta+\sfr12\Sigma}\Pi+3\xi{\mathcal H} \nonumber\\
    &&+\sfr13\bra{\sfr12\phi-2{\cal A}}Q
\label{edot}.
\end{eqnarray}

\textit{Constraint:}
\begin{equation}
{\mathcal H} = 3\xi\Sigma-\bra{2{\cal A}-\phi}\Omega. \label{H}
\end{equation}

Also we give the commutation relation for the {\it dot} amd {\it hat} derivatives, for LRS spacetimes:
\begin{equation}
\hat {\dot \Psi}-\dot {\hat \Psi} = -{\cal A} \dot\Psi+\bra{\sfr13\Theta+\Sigma} \hat\Psi, \label{psihatdot}
\end{equation}
which holds true for any scalar $\Psi$.
\\

\textbf{Perfect fluids:} Now, if we consider a perfect fluid with $Q=\Pi=0$, the consistency conditions of LRS spacetimes (which are derived from the fact that the propagation equations should be identically evolved in time) demand that \begin{equation}
\Omega\xi=0.
\end{equation} Due to this condition, the spacetime is divided into 3 distinct subclasses \cite{Ellis_1968,Elst_Ellis_1996}:
\begin{enumerate}
\item\textbf{LRS class I (Rotating solutions):} $\Omega\neq 0$ \\
$e^b$ is hypersurface orthogonal and $u^b$ is twisting. When $\Omega\neq 0$, we see that $\Rightarrow \xi=0=\Theta=\Sigma, \, \dot f=0$. Therefore models with LRS class I solutions can neither expand nor distort. These models are stationary as the dot of all the scalar quantities vanish.

\item\textbf{LRS class II (The Inhomogeneous orthogonal family):} $\xi=0=\Omega$ \\
Here, \textit{both} $e^b$ and $u^b$ are hypersurface orthogonal. When $\xi=0=\Omega$, there exist 3-surfaces orthogonal to the fluid flow. All models in this dynamic and spatially inhomogeneous LRS class have vanishing `magnetic part' of the {Weyl curvature tensor}.

\item\textbf{LRS class III (Homogeneous orthogonal models with twist):} $\xi\neq 0$ \\
$e^b$ is twisting and $u^b$ is hypersurface orthogonal. When $\xi\neq 0 \Rightarrow \Omega=\phi={\cal A}=0$, \textit{all} spatial derivatives vanish and all scalars \textit{f} are spatially homogeneous. We see that the spacetimes themselves are orthogonally spatially homogeneous (OSH).
\end{enumerate}
\section{A new class of LRS with $\Omega, \xi \neq 0$}

As discussed earlier we would now like to relax the perfect fluid condition, that is we introduce pressure anisotropy and heat flux in the matter, and look for existence of solutions that have both rotation and twist of the preferred direction. To do this, first let us derive an important result for LRS spacetimes.
We can write the full covariant derivatives of the vectors $u^a$ and $e^a$ in terms of the LRS scalars in the following way:
\begin{eqnarray}
\nabla_au_b&=&-{\cal A} u_ae_b+\bra{\sfr13\Theta+\Sigma}e_ae_b \nonumber
\\&&+\bra{\sfr13\Theta-\sfr12\Sigma}N_{ab}+\Omega\varepsilon_{ab}, \label{delaub}
\end{eqnarray}
\begin{eqnarray}
\nabla_ae_b &=& -{\cal A} u_au_b +\bra{\sfr13\Theta+\Sigma}e_au_b \nonumber
\\ && +\sfr12\phi N_{ab}+\xi\varepsilon_{ab}. \label{delaeb}
\end{eqnarray}
Contracting the above with $\varepsilon_{ab}$ and using
\begin{equation}
 \varepsilon^{ab}u_a =\varepsilon^{ab}e_a=0; \; \varepsilon^{ab}\varepsilon_{ab}=2, \label{epabua},
\end{equation}
we get
\begin{equation}
\Omega = \sfr12\varepsilon^{ab}\nabla_au_b, \quad \xi = \sfr12\varepsilon^{ab}\nabla_ae_b . \label{omegaxi}
\end{equation}
Now, for any scalar function `$\Psi$', we have
\begin{equation}
\nabla_b \Psi= -\dot\Psi u_b+ \hat\Psi e_b.
\end{equation}
Differentiating again we have
\begin{eqnarray}
\nabla_a \nabla_b \Psi&=& -\bra{\nabla_a \dot\Psi}u_b - \dot \Psi\bra{\nabla_a u_b} \nonumber
\\&&+\bra{\nabla_a \hat \Psi}e_b+\hat \Psi \bra{\nabla_a e_b} . \label{deladelbf}
\end{eqnarray}
Contracting with $\varepsilon^{ab}$, and noting that $\nabla_a\nabla_b \Psi$ is symmetric in \textit{a} and \textit{b}, we see that the LHS of (\ref{deladelbf}) vanishes. Using equations (\ref{epabua}) and (\ref{omegaxi}) we get an important result:
\begin{equation}\label{scalarcons}
\forall \Psi, \,\,\, \dot\Psi\Omega = \hat\Psi \xi.
\end{equation}
This equation implies self-similarity, for it applies to all scalars, and is unchanged under the transformation $t \rightarrow at,$ $r \rightarrow ar,$ where $t$ and $r$ are the curve parameters of the integral curves of $u$ and $e$.\\

 From the above equation it is clear that if $\Omega\ne0, \xi=0$, the dot derivatives of all the scalars vanish, making the spacetime stationary. On the other hand if $\Omega=0, \xi\ne0$, the hat derivatives of all scalars vanish, making the spacetime spatially homogeneous. Thus we arrive to an important result:
 \begin{thm}
 For LRS spacetimes with non-zero rotation and spatial twist, there always exists a conformal Killing vector in the $[u,e]$ plane. When one of these quantities vanish then the conformal Killing vector becomes a Killing vector. This Killing vector is timelike for vanishing spatial twist and it is spacelike for vanishing rotation. However when both the rotation and spatial twist vanish no such symmetry is guaranteed.
 \end{thm}

 Another important point to be noted here is that $\Omega$ and $\xi$ do not evolve independently. Supposing that both are not equal to zero, first considering $\Psi=\Omega$ in equation
(\ref{scalarcons}) and then $\Psi=\xi$ and using the field equations (\ref{dotomega}), (\ref{dotomega}), (\ref{dotomega}), (\ref{xihat}) and (\ref{H}) we get the constraint
\begin{equation}
\frac{\Omega}{\xi}=-\frac{\phi}{\Sigma-\sfr23\Theta} \,.
\label{constraint1}
\end{equation}

Now to establish the existence of solutions with non-zero rotation and spatial twist, we state and prove the following theorem:
\begin{thm}
Evolution of all the independent geometrical scalars of LRS spacetimes that have non-zero rotation and spatial twist, obey a common second order linear hyperbolic partial differential equation and the existence of a initial spacelike Cauchy surface is guaranteed. Subject to the initial Cauchy data on this surface these geometrical scalars can be uniquely determined, and hence unique solutions of the field equations exist.
\end{thm}
\begin{proof}
Taking the {\it dot} derivative and then the {\it hat} derivative of the equation (\ref{scalarcons}), subtracting them and then using the commutation relation (\ref{psihatdot}) and the  field equations (\ref{dotomega}), (\ref{dotomega}), (\ref{dotomega}), (\ref{xihat}) and (\ref{H}), we obtain the following equation $\forall \Psi,$
\begin{eqnarray}\label{waveeqn}
 -\Omega^2\ddot\Psi+\xi^2\hat{\hat{\Psi}}&-&\dot\Psi\Omega\left[\xi({\cal A}-\phi)+\Omega(\Sigma-\sfr23\Theta)\right]\nonumber\\
&+&\hat{\Psi}\xi\left[2\Omega\Sigma-\sfr13\Omega\Theta-\phi\xi\right]=0.
\end{eqnarray}
We can easily see that the above equation is a hyperbolic (wave like) second order linear partial differential equation  for $\Omega,\xi\ne0$, that governs the evolution of all independent geometrical scalars that describe a LRS spacetime. By the properties of hyperbolic partial differential equations, there exists a unique solution subject to Cauchy initial data on a spacelike Cauchy surface. To check whether such a 3 dimensional spacelike surface exist, let us consider the Lie derivative of the tensor $N^{ab}$ with respect to the spacelike vector $e^a$. We know that
\begin{equation}
(\mathcal{L}_{\bm{e}}N)^{ab}=e^c\nab_cN^{ab}-N^{cb}\nab_ce^a-N^{ac}\nab_ce^b\;.
\end{equation}
Using (\ref{N1}),(\ref{N2}),(\ref{delaub}),(\ref{delaeb}) and (\ref{epabua}) we see that
\begin{equation}
(\mathcal{L}_{\bm{e}}N)^{ab}=\phi N^{ab}\;,
\end{equation}
which implies
\begin{equation}
(\mathcal{L}_{\bm{e}}N)^{ab}u_a=(\mathcal{L}_{\bm{e}}N)^{ab}u_b=0\;,
\end{equation}
that is neither the vector $e^a$ and the tensor $N^{ab}$, nor the Lie derivative of $N^{ab}$ with respect to $e^a$ has any component along the timelike vector $u^a$.
This clearly shows that the tensor product of $e^a$ and $N^{ab}$ indeed spans
a spacelike 3-surface  where we can specify the Cauchy initial data to obtain a unique solution of (\ref{waveeqn}) for all the independent geometrical and thermodynamic scalars of the LRS spacetime.
\end{proof}
Furthermore, the hyperbolic nature of the above equation dictates the existence of two families of characteristics. In analogy with the incoming and outgoing waves, these characteristics describe the expanding/collapsing branches of the solutions.

\subsection{Constraints on thermodynamic variables}

We will now describe the constraints on the thermodynamic variables for the energy momentum tensor of the matter field, that generates LRS solutions with non-zero rotation and spatial twist. We first observe that the common wave like equation (\ref{waveeqn}) was obtained by the Ricci identities of the timelike vector $u^a$ and spacelike vector $e^a$. To obtain the constraints on the matter variables, we need to look at the consistencies of doubly contracted Bianchi identities carefully.
We state and prove the following theorem here:
\begin{thm}
The necessary condition for a LRS spacetime to have non-zero rotation and spatial twist simultaneously is non-zero heat flux which is bounded from both sides.
\end{thm}
\begin{proof}
Taking the \textit{time-like} derivative for the equation (\ref{Ehatmupi}) and using (\ref{psihatdot}) and the field equations, we get
\begin{equation}
\Omega\xi\bra{\mu+p+\Pi}+Q\bra{\Omega^2+\xi^2}=0. \nonumber
\end{equation}
Simplifying the above equation we get
\begin{equation}
\sfr{\sfr\Omega\xi}{1+\bra{\sfr\Omega\xi}^2} = \frac{-Q}{\mu+p+\Pi}.\label{omegaxi2}
\end{equation}
From the above equation it is clear that if we demand both $\Omega$ and $\xi$ are well defined and non-zero, and all the energy conditions to be satisfied we must have $Q \neq 0$. Also it is interesting to note the the ratio of the rotation and spatial twist can be described in terms of the thermodynamic quantities only.
 Now using (\ref{omegaxi2}) to solve for $\sfr\Omega\xi$ gives:
\begin{equation}\label{constraint2}
\frac{\Omega}{\xi} = \frac{-\bra{\mu+p+\Pi}\mp \sqrt{\bra{\mu+p+\Pi}^2-4Q^2}}{2Q}.
\end{equation}
From the above equation it is clear that for the rotation and spatial twist to be well defined, real and non-zero, we must have $(\mu+p+\Pi)^2  > 4Q^2$. Thus the thermodynamic quantities must satisfy the following constraint
\begin{equation}
-\frac12(\mu+p+\Pi) <Q<\frac12(\mu+p+\Pi)\;;\; Q \neq 0 .\label{ConsCon}
\end{equation}
\end{proof}
From \cite{santos,chan} we can see that the above constraints are consistent with the Dominant energy conditions (DEC) for the matter field. Hence we do have matter that obeys the physically reasonable energy conditions that can generate a LRS spacetime with non-zero rotation and spatial twist. Also the rest of the propagation equations evolve identically in time and give no new constraints.

\subsection{Other constraints and solution finding algorithm}

Let us now try to reduce the number of independent geometrical scalars of an LRS spacetime by using equation (\ref{scalarcons}).
Inserting the scalar variables $\phi$, $\bra{\Sigma-\sfr23\Theta}$ and ${\mathcal H}$, and using the field equations, we get the following set of equations:
\begin{eqnarray}
\Omega Q-2\xi^3+2\Omega^2\xi-\sfr13\xi\Theta\Sigma+\xi\Sigma^2+\sfr23\xi\mu && \nonumber \\
-\sfr29\xi\Theta^2 +\xi {\mathcal E}+\sfr12\xi\Pi+A\phi\xi &=& 0, \\
-\Omega A\phi-\sfr13\Omega\Theta\Sigma+\sfr29\Omega\Theta^2-2\Omega^3-\Omega E &&  \nonumber \\
 +\sfr12\Omega\Pi+\sfr13\Omega\mu+\Omega p + \xi\phi\Sigma+2\Omega\xi^2+Q\xi &=& 0, \\
\sfr92\Omega\xi\Sigma^2-3\Omega\xi\Theta\Sigma+\Omega^2Q+\Omega\xi\Pi+\sfr92\xi^2\phi\Sigma && \nonumber \\
\xi\Omega\mu+\xi\Omega p+Q\xi^2 &=& 0.
\end{eqnarray}
Solving the above system of equations for ${\mathcal E}, p$ and $\phi$, we get
\begin{eqnarray}
p &=& -\frac{\Omega^2Q+\xi\Omega\mu+\Omega\xi\Pi+Q\xi^2}{\Omega\xi}, \\
\phi &=& -\frac{\bra{3\Sigma-2\Theta}\Omega}{3\xi},
\end{eqnarray}
which are same as (\ref{omegaxi2}) and (\ref{constraint1}). Also we get a new algebraic relation for ${\mathcal E}$:
\begin{eqnarray}\label{E}
{\mathcal E} &=&  \frac\Omega\xi {\cal A} (\Sigma-\frac23\Theta)-\Sigma^2+\frac13\Theta\Sigma+\frac29\Theta^2\nonumber\\
&&+2(\xi^2-\Omega^2)- \frac\Omega\xi Q -\frac12\Pi-\frac23\mu\;.
\end{eqnarray}

The above equation along with equation (\ref{H}) completely describes the Weyl tensor in LRS spacetimes.

Now taking into account the results for perfect fluid LRS spacetimes \cite{Elst_Ellis_1996}, we see that the above phenomenon is true for any LRS spacetime and we can state this interesting theorem:

\begin{thm}
The symmetry of LRS spacetimes makes the Weyl tensor obey an algebraic constraint with other 1+1+2 geometrical variables.  Hence the doubly contracted Bianchi identities that describe the propagation and evolution of the Weyl tensor become redundant.
\end{thm}

Now we can see that the number of independent geometrical 1+1+2 scalars that describe a LRS spacetime is reduced considerably. For example, specifying \[ {\cal D}_3 :=\brac{{\cal A}, \Theta, \xi, \Sigma, \mu, p(\mu,\Pi,Q), \Pi, Q }, \] will automatically specify \[ {\cal D}_4 :=\brac{\Omega, \phi, {\mathcal E}, \H} ,\] via the constraint equations (\ref{constraint1},\ref{constraint2},\ref{E},\ref{H}). Hence, we can give the initial Cauchy data on any spacelike Cauchy surface for the independent variables using any suitable chosen equation of state $p(\mu, \Pi, Q)$, and then determine their evolution via equation (\ref{waveeqn}), which applies equally to all the variables in ${\cal D}_3$. This will then provide us with a unique self-similar dynamical solution for the LRS spacetimes with non-zero rotation and spatial twist. The nature of the matter required for such solutions to exist follows from Theorem 2, where there is a non trivial condition on the presence of heat flux. There are no other conditions on the density, pressure or pressure anisotropy, but for a physically realistic solutions these must obey the Dominant energy conditions (DEC).

As described in detail in \cite{Ellis_1968}, if a spacetime exhibit local rotational symmetry in an open neighbourhood of a point $P$, then the coordinate freedoms can be used to describe the local metric in the neighbourhood in $(t,r,x,y)$ coordinates in the following way:
\begin{eqnarray}\label{metric}
ds^2&=&-F^2(t,r)dt^2+X^2(t,r)dr^2\nonumber\\
&&+Y^2(t,r)[dx^2+D(x)dy^2]\nonumber\\
&&+g(x)F^2(t,r)[2dt-g(x)dy]dy\nonumber\\
&&-h(x)X^2(t,r)[2dr-h(x)dy]dy
\end{eqnarray}
We can immediately see that $g(x)=h(x)=0$ and $D(x)=\sin^2x$ gives a general spherically symmetric metric which is of LRS class II. However we have already established that LRS spacetimes with non-vanishing rotation and spatial twist must be self similar. Hence the functions $F$, $X$ and $Y$ can be written in terms of a single variable $z\equiv t/r$. Hence only self similar spherically symmetric solutions can be obtained in the limit $g(x)\rightarrow 0, h(x)\rightarrow 0$ or equivalently $\Omega\rightarrow 0,\xi\rightarrow 0$.
Therefore to study the interior of a rotating, radiating and inhomogeneous star as a first approximation from the spherical symmetry, we can start with a self similar spherically symmetric spacetime and add sufficiently small $g(x)$ and $h(x)$, with respect to some covariant scale in the problem (the Misner Sharp mass of the spherical star for example) and solve the field equations with the matter source that obeys all the restrictions as imposed by Theorem 2 and the energy conditions.

\section{Discussion}

In this paper we transparently showed that it is possible to have a Locally Rotationally Symmetric spacetime with non-zero rotation and spatial twist simultaneously if we allow for non-zero and bounded heat flux. We investigated in detail all the covariant geometrical properties of such spacetimes and proved an interesting result that evolution of all the covariant scalars obey a single common hyperbolic linear second order partial differential equation. The existence of spacelike Cauchy surface, where initial Cauchy data can be provided is guaranteed. It was also shown that these solutions are self similar as they possess a conformal Killing vector in the $[u,e]$ plane.\\

As these solutions are neither stationary nor spatially homogeneous in general, with suitable equations of state, perhaps with the temperature $T$ as an internal variable in the equations of state for $P$, $\Pi$, and $Q$, they have the potential to give exact general relativistic models for rotating and dynamic and radiating stellar structures as they definitely have non zero heat flux in the interior. These solutions will then provide a relativistic description of a rotating stellar interior
with quadrupole and other higher multipole moments and this may account for physical features of stars that cannot be explained by Newtonian dynamics.

 \begin{acknowledgments}
SS. GFRE and RG are supported by National Research Foundation (NRF), South Africa. SDM
acknowledges that this work is based on research supported by the South African Research Chair Initiative of the Department of
Science and Technology and the National Research Foundation.
\end{acknowledgments}