\section{\label{sec1}Introduction}
Temperature is an important parameter of physical systems. For many problems, the temperature scale of interest is far below the optical gap (or the excitation energies). For example, typical electronic excitation energies in molecules are of the order of a few electron volts (or $\mathrm{eV}$) , which is much larger than room temperature ($\sim 25 \mathrm{meV}$). In such problems where we are interested only in electronic degrees of freedom, it suffices to know the ground electronic state and perhaps a few low-lying excited states. To access these states, we solve the time-independent Schr\"{o}dinger equation. As is well known,\cite{dirac_quantum_1929} this is a very complicated problem and one generally relies on a series of approximate methods such as Hartree-Fock, perturbation theory, configuration interaction (CI), coupled cluster theory\cite{crawford_introduction_2000,bartlett_coupled-cluster_2007} (CC), Monte Carlo methods,\cite{hammond_monte_1994, foulkes_quantum_2001, zhang_quantum_2003, zhang_quantum_2004, al-saidi_auxiliary-field_2006} among many more.\cite{white_density-matrix_1993,ostlund_thermodynamic_1995,georges_dynamical_1996}

There are, however, many interesting problems and applications where one may be interested in temperature scales that are comparable or even larger than the excitation gap. Examples include metallic compounds with small gap that host an unconventional superconductivity persisting at relatively high temperatures,\cite{lee_doping_2006} ultra-cold chemistry,\cite{balakrishnan_perspective_2016,bohn_cold_2017} geochemical processes which generally involve very high temperatures and pressure,\cite{guillot_interiors_1999} etc. In these problems, we can no longer make do with a few electronic states and must evaluate properties as thermal averages weighted over an appropriate ensemble of states. For a system in thermal equilibrium at inverse temperature $\beta$, the ensemble of choice is generally canonical or grand-canonical and the expectation value of an observable $A$ is defined by
\begin{equation}
    \langle A \rangle_\beta = \frac{1}{\mathcal{Z}} \mathrm{Tr} ( A \rho ),
    \label{trace-avg}
\end{equation}
where $\rho$ is the thermal density matrix, and $\mathcal{Z}$ is the partition function, given by
\begin{equation}
    \rho = e^{-\beta H'},  \quad \mathcal{Z} = \mathrm{Tr} ( \rho ),
\end{equation}
with $H' = H$ for the canonical ensemble and $H' = H-\mu N$ for the grand canonical ensemble, where $\mu$ is the chemical potential. In this paper, we shall work explicitly with the grand canonical ensemble.

Exactly computing $\rho$ (or equivalently $\mathcal{Z}$) requires information about the entire spectrum of the Hamiltonian, which is far from feasible. Accordingly, just as for zero temperature, a series of approximate methods is needed. Several methods have been proposed over the years to evaluate thermal averages of operators.\cite{mermin_stability_1963,sokoloff_consequences_1967,lichtenstein_finite-temperature_2001,verstraete_matrix_2004,feiguin_finite-temperature_2005, white_minimally_2009, stoudenmire_minimally_2010,pittalis_exact_2011, pizorn_real_2014, hermes_finite-temperature_2015, czarnik_variational_2016, santra_finite-temperature_2017, zgid_finite_2017, claes_finite-temperature_2017, white_time-dependent_2018, hummel_finite_2018} Most of these methods can be broadly categorized into deterministic methods such as diagrammatic perturbation theory based on the Matsubara formalism,\cite{matsubara_new_1955} and stochastic quantum Monte Carlo methods.
Wavefunction methods are particularly convenient in the study of zero-temperature ground-state properties of finite sized systems and clearly, their thermal equivalents are highly desirable. However, the development of such thermal wavefunctions methods has been rather challenging, primarily because the thermal density matrix cannot be expressed in terms of a single wavefunction in the original Hilbert space. Nevertheless, several such methods have been introduced over the years such as the Ancilla density matrix renormalization group\cite{verstraete_matrix_2004,feiguin_finite-temperature_2005} and finite-temperature perturbation theories,\cite{he_finite-temperature_2014,santra_finite-temperature_2017} among many more.

Given its features, especially size-extensivity and success with weakly correlated systems, the coupled cluster ansatz is an ideal candidate to study finite-temperature properties.
A thermal analogue of the CC method\cite{sanyal_thermal_1992,sanyal_systematic_1993,mandal_thermal_1998,mandal_finite-temperature_2003} was proposed by Mukherjee \emph{et. al.} and has been further elaborated recently in independent works by White \emph{et al.}\cite{white_time-dependent_2018} and Hummel.\cite{hummel_finite_2018} This formulation uses a thermal Wick's theorem to compactly represent the imaginary time evolution operator as a \emph{thermal normal ordered} exponential of some cluster operator and a number.

In this paper, we present an alternative approach to thermal coupled cluster based on the thermofield dynamics (TFD),\cite{matsumoto_thermo_1983,semenoff_functional_1983,umezawa_methods_1984,evans_heisenberg_1992} using a framework we recently explored in Ref.~\citenum{harsha_thermofield_2019}. Thermofield dynamics provides a convenient way to represent the thermal density matrix via a wavefunction which evolves in temperature according to the imaginary-time evolution Schr\"odinger equation. Undoubtedly, TFD has the potential to study many-electron systems in quantum chemistry and condensed-matter physics.\cite{suzuki_thermo_1985, hatsuda_mean_1989, walet_thermal_1990, de_vega_thermofield-based_2015,borrelli_quantum_2016, nocera_symmetry-conserving_2016, chen_finite_2017, borrelli_simulation_2017, wu_variational_2018} A connection between thermofield formalism and thermal cluster-cumulant theory was also explored by Mukherjee \emph{et. al.} in Ref.~\citenum{sanyal_systematic_1993}. Here, we parametrize this so-called thermal wavefunction as an exponential ansatz to integrate the evolution equation.

\section{\label{sec2}Coupled Cluster Theory}

The coupled cluster method is one of the most widely used methods in quantum chemistry. Introduced first in nuclear physics by Coester and K\"{u}mmel\cite{coester_bound_1958, coester_short-range_1960} and later reformulated for electronic structure theory by \v{C}i\v{z}ek and Paldus,\cite{cizek_correlation_1971} it uses an exponential wavefunction ansatz
\begin{equation}
    |\Psi \rangle = e^{T} |\Phi \rangle,
\end{equation}
to solve the time-independent Schr\"odinger equation. Here $|\Phi\rangle$ is some Hartree-Fock (HF) Slater determinant reference and $T$ contains particle-hole excitations defined on this reference state,
\begin{subequations}
    \begin{align}
        T &= T_1 + T_2 + \ldots, \\
        T_1 &= \sum_{i,a} t_i^a c^\dagger_a c_i,\\
        T_2 &= \frac{1}{4} \sum_{i,j,a,b} t_{ij}^{ab} c^\dagger_a c^\dagger_b c_j c_i, \ldots
    \end{align}
\end{subequations}
where we have followed the standard notation for labeling orbital indices, i.e., occupied orbitals are denoted by indices $i,j,k,\ldots $, while unoccupied orbitals are denoted by $a, b, c, \ldots$. Since the wavefunction $\vert \Psi \rangle$ is constructed as excitations on a single Slater determinant $\vert \Phi \rangle$, the traditional CC is also known as single-reference CC. In what follows, unless explicitly mentioned, we consider only the single-reference CC. The equations for the ground state energy and the unknown $t$-amplitudes can be obtained by left-projecting the Schr\"odinger equation
\begin{equation}
    \bar{H} \left| \Phi \right \rangle = E \left | \Phi \right \rangle, \quad \bar{H} = e^{-T} H e^{T}
    \label{schro-eq}
\end{equation}
with various Slater determinants. For example, if the cluster operator $T$ is truncated to single and double excitations only (CCSD), the energy and amplitude equations are obtained by solving the following equations
\begin{align}
    E &= \left \langle \Phi \right \rvert \bar{H} \left \lvert \Phi \right \rangle, \\
    0 &= \left \langle \Phi_{i}^{a} \right \rvert \bar{H} \left \lvert \Phi \right \rangle, \\
    0 &= \left \langle \Phi_{ij}^{ab} \right \rvert \bar{H} \left \lvert \Phi \right \rangle,
\end{align}
where $\langle \Phi_{i}^{a} \vert$ and $\langle \Phi_{ij}^{ab} \vert $ are singly- and doubly-excited slater determinants.

As introduced so far, CC describes a correlated ansatz for the \textit{ket} wavefunction. In order to compute expectation values other than that of the Hamiltonian, one also needs a correlated \textit{bra} state. A linear response wavefunction is generally employed for this purpose,\cite{arponen_variational_1983,helgaker_analytical_1988} i.e. one makes the energy functional
\begin{equation}
    E = \langle \Phi \vert (1 + Z) e^{-T} H e^{T} \vert \Phi \rangle
    \label{lin-resp-cc-func}
\end{equation}
stationary with respect to $T$ and $Z$, where
\begin{subequations}
    \label{Z-operator}
    \begin{align}
        Z &= Z_1 + Z_2 + \ldots, \\
        Z_1 &= \sum_{i,a} z_i^{a} c_i^\dagger c_a, \\
        Z_2 &= \frac{1}{4}\sum_{i,a} z_{ij}^{ab} c_i^\dagger c_j^\dagger c_b c_a.
    \end{align}
\end{subequations}
By realizing that the cluster operator $T$ and CI operator $Z$ are composed of particle-hole excitation and de-excitation operators respectively, the \textit{bra} state can be re-written as an explicit CI wavefunction
\begin{equation}
    \langle \Psi'\vert = \langle \Phi  \vert (1 + Z) e^{-T} = \langle \Phi \vert (1 + W) e^{w_0},
    \label{ci2cc}
\end{equation}
where $w_0$ is a constant and $W$ has the same operator-form as $Z$.

A similar formulation known as equation of motion CC\cite{rowe_equations--motion_1968} (EOM-CC) can be used for excited states. Coupled cluster gives highly accurate results for weakly correlated systems and its success can be attributed to the polynomial computational scaling ($\mathcal{O}(N^{6})$ for CCSD) as well as to the fact that the computed properties are size-extensive (scale linearly with the particle number) in the thermodynamic limit.

\section{\label{sec3}Thermofield Dynamics}
Thermofield dynamics is a real-time thermal field theory that treats both time and temperature on equal footing, and was proposed as an alternative to the Keldysh approach in the Matsubara imaginary time formalism. It provides a prescription for purification of the thermal density matrix, allowing us to construct a wavefunction $\vert \Psi(\alpha, \beta) \rangle$, with $\alpha=\beta \mu$, generally known as the \textit{thermal vacuum}, \textit{thermofield double state} or simply as the \textit{thermal state}, so that the trace over an ensemble of states in Eq.~\ref{trace-avg} can be replaced by an expectation value over this wavefunction, i.e.
\begin{equation}
    \langle A \rangle = \frac{ \langle \Psi (\alpha, \beta) \vert A \vert \Psi (\alpha, \beta) \rangle }{ \langle \Psi (\alpha, \beta) \vert \Psi (\alpha, \beta) \rangle }.
    \label{tfd-exp}
\end{equation}
That the thermal state $\vert \Psi (\alpha, \beta) \rangle$ cannot be a pure state in the physical Hilbert space $\mathcal{H}$ is easily established since the density matrix $\rho$ represents mixed states to begin with. In TFD, one therefore introduces a fictitious, conjugate copy of the original Hilbert space, known as the \textit{tilde}-conjugate space or $\tilde{\mathcal{H}}$, and the thermal state is then defined in the doubled space $\mathcal{H} \otimes \tilde{\mathcal{H}}$ as
\begin{equation}
    \vert \Psi (\alpha, \beta) \rangle = e^{(\alpha N-\beta H)/2} \vert \mathbb{I} \rangle, \quad \vert \mathbb{I} \rangle = \sum_m \vert m, \tilde{m} \rangle,
    \label{thermal-vacuum}
\end{equation}
where $H$ is the Hamiltonian, $\{\vert m \rangle \}$ is some orthonormal basis in the Hilbert space $\mathcal{H}$ and $\vert \tilde{m} \rangle $ is the tilde-state corresponding to $\vert m \rangle$, while we have used the shorthand notation $\vert m, \tilde{m}\rangle = \vert m \rangle \otimes \vert \tilde{m} \rangle$. The state $\vert \mathbb{I} \rangle$ is the state with maximal entanglement between $\mathcal{H}$ and $\tilde{\mathcal{H}}$ and is the exact thermal state at infinite temperature (i.e. $\beta=0$) and $\alpha=0$. Moreover, it is invariant under any transformation of the basis. Further details about the TFD formalism and the structure of the \textit{tilde}-conjugate space $\tilde{\mathcal{H}}$ can be found in Ref.~\onlinecite{harsha_thermofield_2019} and the references therein.

By construction, the thermal state satisfies the following imaginary time evolution equations,
\begin{subequations}
    \label{beta-mu-schro}
    \begin{align}
        \frac{ \partial }{\partial \beta} \vert \Psi (\alpha, \beta) \rangle &= -\frac{1}{2} H \vert \Psi (\alpha, \beta) \rangle.
        \label{beta-schro-eq}
        \\
        \frac{ \partial }{\partial \alpha} \vert \Psi (\alpha, \beta) \rangle &= \frac{1}{2} N \vert \Psi (\alpha, \beta) \rangle.
        \label{alpha-schro-eq}
    \end{align}
\end{subequations}
One can solve for $\vert \Psi (\alpha, \beta)\rangle $ by integrating Eq.~\ref{beta-mu-schro}, generally starting from $(\alpha,\beta)=(0,0)$, where the initial thermal state is known exactly.
Exactly evolving the thermal state $\vert \Psi (\alpha, \beta)\rangle$ is equivalent to the exact diagonalization of the density matrix $\rho$. Clearly, approximations need to be introduced in the process.

The simplest approximation that we can invoke is to use the mean-field Hamiltonian $H_0$ instead of $H$. For a many-electron system, the state $\vert \mathbb{I}\rangle$ can be expressed in terms of single-particle Fock states,
\begin{equation}
    \vert \mathbb{I} \rangle = \prod_{p\in \mathrm{levels}} \left ( \vert 0,\tilde{0}\rangle_p + \vert 1,\tilde{1}\rangle_p \right ),
\end{equation}
where $\vert 0 \rangle_p$ and $\vert 1 \rangle_p$ mean that the orbital $p$  is empty or occupied, respectively.
Moreover, if we chose to work with the eigen-basis of $H_0$, i.e.
\[
    H_0 = \sum_p \epsilon_p c^\dagger_p c_p
\]
the normalized mean-field thermal state can be written as
\begin{equation}
    \vert 0 (\alpha, \beta) \rangle = \prod_{p \in \mathrm{levels}} \left ( x_p \vert 0, \tilde{0} \rangle_p + y_p \vert 1, \tilde{1} \rangle \right ),
    \label{mf-thermal-state}
\end{equation}
where $x_p$ and $y_p$ are related to the Fermi-Dirac statistics,
\begin{subequations}
    \begin{align}
        x_p &= \frac{1}{\sqrt{1 + e^{(\alpha - \beta \epsilon_p)}}}, \\
        y_p &= \frac{ e^{(\alpha - \beta \epsilon_p)/2} }{ \sqrt{1 + e^{(\alpha - \beta \epsilon_p)}} },
    \end{align}
\end{subequations}
with $x_p^2 + y_p^2=1$. The  mean-field thermal state in Eq.~\ref{mf-thermal-state} allows us to introduce a thermal Bogoliubov transformation,
\begin{equation}
    \begin{bmatrix}
        a_p  \\
        \tilde{a}_p^\dagger
    \end{bmatrix} =
    \begin{bmatrix}
        x_p & -y_p\\
        y_p & x_p
    \end{bmatrix}
    \begin{bmatrix}
        c_p\\
        \tilde{c}_p^\dagger
    \end{bmatrix},
    \label{hfb-transformation}
\end{equation}
such that
\[
    a_p \vert 0 (\alpha, \beta) \rangle = 0 = \tilde{a}_p \vert 0(\alpha, \beta) \rangle.
\]
Correlated methods can be built with either $\vert \mathbb{I} \rangle$ or $\vert 0(\alpha, \beta) \rangle$ as the reference, while integrating Eq.~\ref{beta-mu-schro}. As we have discussed in Ref~\citenum{harsha_thermofield_2019}, the former choice, which we call the fixed-reference formalism, performs well only in the vicinity of $\beta=0$, while the latter, called the covariant formalism, yields accurate results for the entire range of $\beta$. In the next section, we present details to integrate Eq.~\ref{beta-mu-schro} through the CC ansatz in the covariant formalism. We have also included a short discussion on the fixed-reference approach in Appendix~\ref{app2:fixed}. For brevity of notation, we will henceforth refer the  mean-field thermal state $\vert 0(\alpha, \beta)\rangle$ by $\vert \Psi_0 \rangle$.

\section{\label{sec4}Thermal Coupled Cluster}

As explained in Eq.~\ref{lin-resp-cc-func}, the CC expectation value of any operator $A$ can be evaluated as an asymmetric expectation value,
\begin{equation}
    \langle A \rangle_{\mathrm{CC}} = \frac{ \langle \Psi' \vert A \vert \Psi \rangle }{
        \langle \Psi' \vert \Psi \rangle
    },
\end{equation}
where both the \textit{ket} $\vert \Psi \rangle$ and the \textit{bra} $\langle \Psi' \vert$ states are approximations to the same thermal state, and consequently evolve according to Eq.~\ref{beta-mu-schro} and its adjoint respectively.

Given that the Bogoliubov transformation in Eq.~\ref{hfb-transformation} is BCS-like, i.e. orbitals $p$ (in $\mathcal{H}$) and $\tilde{p}$ (in $\mathcal{\tilde{H}}$) in the mean-field thermal state (Eq.~\ref{mf-thermal-state}) are coupled in just the same way as two opposite momentum single-particle levels in a BCS-wavefunction, we parametrize the \textit{ket} state as an exponential of quasiparticle creation operators\cite{henderson_quasiparticle_2014} acting on an $\alpha$- and $\beta$-dependent mean-field thermal reference, $\vert \Psi_0 \rangle$, defined in Eq.~\ref{mf-thermal-state},
\begin{subequations}
    \begin{align}
        \vert \Psi \rangle &= e^{S(\alpha, \beta)} \vert \Psi_0 \rangle,
        \label{covar-cc-ansatz}
        \\
        S &= s_0
        + \sum_{p,q} s_{pq} a^\dagger_p \tilde{a}^\dagger_p
        + \frac{1}{(2!)^2} \sum_{p,q,r,s} s_{pqrs} a^\dagger_p a^\dagger_q \tilde{a}^\dagger_s \tilde{a}^\dagger_r
        + \ldots
        \label{covar-cc-cluster}
    \end{align}
\end{subequations}
On the other hand, the \textit{bra} state, as in the traditional CC formalism, is approximated as a linear CI-like wavefunction, i.e.
\begin{equation}
    \langle \Psi'\vert = \langle \Psi_0 \vert \: (1 + Z)e^{z_0}e^{-S},
\end{equation}
and as explained in Eq.~\ref{ci2cc}, it can be expressed as an effective CI wavefunction
\begin{subequations}
    \begin{align}
        \langle \Psi'\vert &= \langle \Psi_0 \vert (1 + W) e^{w_0}, \\
        W &= \sum_{p,q} w_{pq} \tilde{a}_q a_p + \frac{1}{4} \sum_{p,q,r,s} w_{pqrs} \tilde{a}_r \tilde{a}_s a_q a_p + \ldots
    \end{align}
\end{subequations}
For both the \textit{bra} and the \textit{ket} states, the reference $\vert \Psi_0 \rangle$ evolves continuously as we evolve the Schr\"odinger Eq.~\ref{beta-mu-schro} (hence the name `covariant'). Accordingly, both the amplitudes ($s_{pq}$, $w_{pq}$, etc.) and the quasiparticle operators $a^\dagger, \tilde{a}^\dagger$ carry $\alpha$- and $\beta$-dependence.
The $\alpha$- and $\beta$-evolution of $\langle \Psi' \vert$, a thermal CI wavefunction, is governed by
\begin{subequations}
    \label{cov-ciAll-beta}
    \begin{align}
        \langle \Psi_0 \vert \left ( \frac{\partial W}{\partial \alpha} + (1 + W) \frac{\partial w_0}{\partial \alpha} \right )
        &= \frac{1}{2} \langle \Psi_0 \vert  N_{\textrm{CI}},
        \\
        \langle \Psi_0 \vert \left ( \frac{\partial W}{\partial \beta} + (1 + W) \frac{\partial w_0}{\partial \beta} \right )
        &= -\frac{1}{2} \langle \Psi_0 \vert H_{\textrm{CI}},
    \end{align}
\end{subequations}
where $H_{\textrm{CI}}$ and $N_{\mathrm{CI}}$ are effective CI Hamiltonian and Number operators respectively, and are given by
\begin{align*}
    H_{\textrm{CI}} &= (1 + W)\,H - H_0\,(1 + W) ,
    \\
    N_{\textrm{CI}} &= W\, N - N\, W.
\end{align*}
A detailed discussion on thermal CI and the derivation of these equations can be found in Ref.~\onlinecite{harsha_thermofield_2019}.

For the evolution of $\vert \Psi \rangle$, substituting the CC ansatz from Eq.~\ref{covar-cc-ansatz} into the Schr\"odinger Eq.~\ref{beta-mu-schro} gives
\begin{subequations}
    \label{covar-cc-gov}
    \begin{align}
        e^{-S} \left ( \frac{\partial}{\partial \alpha} e^{S} \right )\: \vert \Psi_0 \rangle &=
        \frac{1}{2} \left ( e^{-S} N e^{S} - N \right ) \: \vert \Psi_0 \rangle,
        \\
        e^{-S} \left ( \frac{\partial}{\partial \beta} e^{S} \right )\: \vert \Psi_0 \rangle &=
        - \frac{1}{2} \left ( e^{-S} H e^{S} - H_0 \right ) \: \vert \Psi_0 \rangle.
    \end{align}
\end{subequations}
The evolution equations for the amplitudes can be obtained by left projecting Eq.~\ref{covar-cc-gov} with the respective determinants.

\subsection{Wilcox identity}
The process of reducing Eq.~\ref{covar-cc-gov} to evolution equations for the amplitudes is complicated by the fact that the derivative of the cluster operator does not commute with the operator itself, i.e.
\[
    \left [ \frac{\partial S}{\partial x} , S \right ] \neq 0,
\]
where $x = \alpha, \beta$. The derivative of the exponential cluster operator is appropriately performed by making use of the Wilcox identity,\cite{wilcox_exponential_1967} which states that the derivative of the exponential of an operator $M$ with respect to some parameter $\lambda$ can be evaluated as
\begin{equation}
    \frac{\partial }{\partial \lambda}e^{M(\lambda)} = \int_0^1 dy \: e^{(1-y)M} \frac{\partial M}{\partial \lambda} e^{y M}.
    \label{wilcox}
\end{equation}
With this, the left-hand side of Eq.~\ref{covar-cc-gov} becomes
\begin{subequations}
    \begin{align}
        e^{-S} \left ( \partial_x \: e^{S} \right ) &= \int_0^1 dy \: e^{-y S} \: (\partial_x S) \:  e^{y S}, \label{wilcox-demo-1}
        \\
        &= (\partial_x S) + \frac{1}{2!} [ (\partial_x S), S ] \nonumber
        \\
        & \quad + \frac{1}{3!} \left[ [(\partial_x S), S], S \right] + \ldots, \label{wilcox-demo-2}
    \end{align}
\end{subequations}
where we have used the shorthand $\partial_x$ for $\partial/\partial x$, and made use of the Baker-Campbell-Hausdorff expansion in going from Eq.~\ref{wilcox-demo-1} to \ref{wilcox-demo-2}. Finally, breaking the derivative $\partial_x S$ into the amplitude ($\partial_{\mathrm{amp}}S$) and operator ($\partial_{\mathrm{op}} S$) derivatives,
\[
    \partial_x S = \partial_{\mathrm{amp}}S + \partial_{\mathrm{op}} S,
\]
and realizing that the former commutes with $S$, we can compactly write the left-hand side of Eq.~\ref{covar-cc-gov} as
\begin{equation}
    e^{-S} \left ( \partial_x \: e^{S} \right ) = \frac{\partial_{\mathrm{amp}} S}{\partial x}+ S_x, \quad
    S_x = \int_0^1 dy \: e^{-yS}  \frac{ \partial_{\mathrm{op}} S}{ \partial x}   e^{yS},
\end{equation}
where the integration over $y$ in the second equation can be carried out analytically, as explained in Eq.~\ref{wilcox-demo-2}. With these details, Eq.~\ref{covar-cc-gov} can be further simplified as
\begin{subequations}
    \label{covar-cc-gov2}
    \begin{align}
        \frac{\partial_{\mathrm{amp}} S}{\partial \alpha}\: \vert \Psi_0 \rangle &= \left [ \frac{1}{2} \left ( e^{-S} N e^{S} - N \right ) - S_{\alpha} \right ] \: \vert \Psi_0 \rangle
        \\
        \frac{\partial_{\mathrm{amp}} S}{\partial \beta}\:  \vert \Psi_0 \rangle &= -\left [ \frac{1}{2} \left ( e^{-S} H e^{S} - H_0 \right ) + S_{\beta} \right ] \: \vert \Psi_0 \rangle,
    \end{align}
\end{subequations}
which can then be left-projected with various thermal quasiparticle states to yield a set of differential equations governing the evolution of the $s$-amplitudes in the chemical potential - temperature or $\alpha$-$\beta$ space. Complete expressions for the CCSD evolution equations are included in Appendix \ref{app1:eqns}.

\begin{figure*}[!ht]
    \subfloat[2-site, $U/t=1$\label{fig:2SitesHalfFilling}]{
        \centering
        \includegraphics[width=0.48\textwidth]{fig1a_hub_2s_u1.pdf}
    }
    \subfloat[6-site, $U/t=2$\label{fig:6SitesHalfFilling}]{
        \centering
        \includegraphics[width=0.48\textwidth]{fig1b_hub_6s_u2.pdf}
    }
    \caption{Error in internal energy for thermal HF, covariant CISD, fixed-reference and covariant thermal CCSD for (a) two-site, and (b) six-site Hubbard models with $U/t = 1,2$ respectively at half filling on average.}
    \label{fig:HubbardHalfFilling}
\end{figure*}

\begin{figure*}[t]
    \subfloat[6-levels, $G=0.2$\label{fig:6LevelHalfFilling_G2}]{
        \centering
        \includegraphics[width=0.48\textwidth]{fig2a_pairing_gpoint2.pdf}
    }
    \subfloat[6-levels, $G=0.5$\label{fig:6LevelHalfFilling_G5}]{
        \centering
        \includegraphics[width=0.48\textwidth]{fig2b_pairing_gpoint5.pdf}
    }
    \caption{Error in internal energy for thermal HF, covariant CISD and CCSD for the six-level pairing model with (a) $G = 0.2$, and (b) $G = 0.5$ respectively at half filling on average.}
    \label{fig:pairingHalfFilling}
\end{figure*}

\begin{figure}[t]
    \includegraphics[width=0.96\columnwidth]{fig3_be_sto3g.pdf}
    \caption{Error in internal energy for Be atom as function of temperature ($\theta$) in thermal HF, covariant CISD and CCSD as calculated in STO-3G basis.}
    \label{fig:Beryllium}
\end{figure}

\section{Implementation}
\subsection{Integration}
An explicit Runge-Kutta method of order (4)5 with step-size control\cite{dormand_family_1980} has been employed for integrating the resulting set of differential equations, where step-size is adaptively modified so as to keep the relative error in evolution within $10^{-8}$. Furthermore, for applications that require a fixed average filling fraction, a tolerance of $10^{-5}$ in the number of electrons is used to find the target chemical potential.

\subsection{Choice of $H_0$}
As mentioned earlier, for the covariant thermal CC, a mean-field Hamiltonian $H_0 = \sum_p \epsilon_p a_p^\dagger a_p$ is used to construct the mean-field thermal reference state $\vert \Psi_0 \rangle$. There are several alternatives to pick this mean-field Hamiltonian -- in thermal HF theory,\cite{mermin_stability_1963, sokoloff_consequences_1967} one uses an $H_0$ that optimizes the applicable free energy, and it therefore depends on $\alpha$ and $\beta$. Such a choice, while it may yield an excellent mean-field thermal reference state, will result in convoluted expressions for the mean-field evolution as well as the operator derivatives, and is beyond the scope of current work.

In our implementation, we use zero-temperature ground-state HF, which may or may not preserve the ground-state symmetries, to construct the energy eigenvalues in $H_0$. Using an $\alpha$- and $\beta$-independent $H_0$ is analogous to the imaginary-time interaction picture formalism and simplifies our implementation since, (\emph{i}) an $\alpha$- and $\beta$-independent $H_0$ leads to clean analytical forms for the mean-field evolution operator as well as the $\alpha$- and $\beta$-derivatives of the thermal quasiparticle operators, and (\emph{ii}) $H_0$ being diagonal yields a straightforward thermal Bogoliubov transformation in Eq~\ref{hfb-transformation}. We note, however, that this choice is different from the one used in Ref.~\onlinecite{harsha_thermofield_2019}, where we use simply the one-electron Hamiltonian to construct $H_0$.

\section{\label{sec5}Results}

Armed with the working equations, we now proceed to present results for the application of the thermal CC in the covariant formalism, truncated at singles and doubles (CCSD), to various many-electron systems, \emph{viz.} the one-dimensional Hubbard model,\cite{hubbard_electron_1963} the pairing or the reduced BCS model, as well as chemical systems (atomic Beryllium and molecular $\mathrm{H}_2$).
In order to make correspondence with the canonical ground state limit, we present results for all of these systems with a fixed number of particles on average - at each $\beta$ grid-point, we evolve the thermal states in $\alpha$ to fix the average number of particles before evolving again in $\beta$.
We compare our results with full configuration interaction (FCI) results.

We first apply the CCSD methods to the one-dimensional Hubbard model with periodic boundary conditions. Having already presented results for thermal CI truncated to singles and doubles (CISD) in Ref.~\onlinecite{harsha_thermofield_2019}, this model system seems to be the right place to start comparing thermal CC with thermal HF and CI. The Hamiltonian is given by
\begin{equation}
    H = -t \sum_{\langle p,q \rangle,  \sigma} \left(c_{p,\sigma}^\dagger \, c_{q,\sigma} + \textrm{h.c.}\right) + U \, \sum_{p} n_{p,\uparrow} \, n_{p,\downarrow},
\end{equation}
where $\langle,\rangle$ denotes that the sum is carried over sites connected in the lattice, $t$ denotes the strength of the kinetic energy term, $U$ denotes the strength of the on-site Coulomb repulsion, and $n_{p,\sigma} = c_{p,\sigma}^\dagger \, c_{p,\sigma}$ is the number operator for lattice site $p$ and spin $\sigma$.  The ratio $U/t$ characterizes the correlation strength.

\begin{figure}[t]
    \includegraphics[width=0.96\columnwidth]{fig4_h2_074ang_sto3g.pdf}
    \caption{Error in internal energy for H$_2$ molecule at bond length of $0.74 \textup{\AA}$ as function of temperature $\theta$ in thermal HF, covariant CISD and CCSD as calculated in STO-3G basis.}
    \label{fig:H2_sto3g}
\end{figure}

Figure~\ref{fig:2SitesHalfFilling} shows the temperature dependence of the error in internal energy for a two-site Hubbard model at half-filling on average with $U/t=1$ as computed by thermal HF, covariant CISD, thermal CCSD as well its fixed-reference formulation. Figure~\ref{fig:6SitesHalfFilling} presents the same for a six-site Hubbard model with $U/t=2$. For both cases, restricted HF (or RHF) eigenvalues and integrals have been used to construct $H_0$ and $H$ respectively. Thermal CCSD clearly outperforms CISD, especially for larger systems where CI is a less accurate wavefunction ansatz. Moreover, the covariant thermal CISD and CCSD go to the appropriate ground-state restricted CCSD in the limit $\beta \rightarrow \infty$ (or $\theta \rightarrow 0$, where $\theta = 1/\beta$ is the temperature). We note that while CCSD and CISD are exact in describing the ground state of the two-site Hubbard model, they are not exact at finite temperatures since we are working in the grand canonical ensemble.
The fixed-reference CCSD, on the other hand, performs poorly except for a small window around $\beta = 0$. Such a behaviour can be anticipated since this method uses the thermal reference corresponding to $\beta = 0$ as its starting point. Accordingly, for all other results that follow, we present only the covariant methods.

Next, we consider the reduced BCS or the pairing model, the Hamiltonian for which is given by
\begin{equation}
    H = \sum_p \epsilon_p N_p - G \sum_{p,q} P^\dagger_p P_q
\end{equation}
where $N_p$ counts the number of electrons, $\epsilon_p$ denotes the energy, and $P^\dagger_p$/$P_p$ respectively creates/annihilates a pair of electrons in the $p^{\mathrm{th}}$-level, while $G$ quantifies the attractive pair-hopping interaction. Here we choose the energy levels with a uniform spacing of 1 unit, i.e., $\Delta \epsilon = \epsilon_{p+1}-\epsilon_p = 1$.
Figure~\ref{fig:6LevelHalfFilling_G2} describes the temperature dependence of the error in internal energy for a six-level pairing model with $G = 0.2$ (weakly correlated) at average half-filling. Figure~\ref{fig:6LevelHalfFilling_G5} shows the same for $G = 0.5$ (near critical regime). Again, we use RHF eigenvalues and integrals to construct $H_0$ and $H$, and we see that the covariant thermal CCSD improves significantly over both the HF and CISD.
We also recover the zero-temperature ground-state limit for thermal CISD and CCSD.
Figures~\ref{fig:Beryllium} and \ref{fig:H2_sto3g} show similar trends for atomic Beryllium and molecular H$_2$ at bond length $0.74 \textup{\AA}$ in STO-3G basis sets.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.96\columnwidth]{fig5_uhf_vs_rhf_6x1_u5.pdf}
    \caption{Error in internal energy for thermal HF and thermal CCSD, both based on RHF as well as UHF, for six-site Hubbard model with $U/t=5$ at half filling on average.}
    \label{fig:HubbardStrongCorr}
\end{figure}

\begin{figure}[tbh!]
    \centering
    \includegraphics[width=0.96\columnwidth]{fig6_sz_corr_10x1_u2_fixedMu.pdf}
    \caption{Trends in $z$-component spin-spin correlation function for the ten-site Hubbard model with $t=1,\, U=2$ at various temperatures $\theta$ and a fixed chemical potential, $\mu = 0.9$ (which ensures half-filling in the large $\beta$ limit), computed using thermal CCSD. Exact (FCI) and RCCSD results for the ground-state correlation functions are also included for reference.}
    \label{fig:spin-corr}
\end{figure}

In all of our applications so far, $H_0$ is constructed using RHF energy eigenvalues. Accordingly, in the zero-temperature limit, thermal HF approaches ground-state RHF and thermal CCSD approaches ground-state RCCSD. In strongly-correlated regime, where zero-temperature RCCSD fails to converge, the evolution of thermal CCSD also eventually diverges. Such an issue can be avoided by using unrestricted HF (or UHF) energy eigenvalues and integrals to construct $H_0$ and $H$ respectively. Figure~\ref{fig:HubbardStrongCorr} shows the performance of RHF and UHF based thermal CCSD internal energies against FCI results for six-site Hubbard model at half-filling on average with $U/t=5$.

In addition to the internal energy, we can also compute other physical properties and correlation functions at any temperature / chemical potential. In Figure~\ref{fig:spin-corr}, we show the $z$-component spin-spin correlation function
\begin{equation}
    \chi(i,j) = \langle S^z(i)\:S^z(j) \rangle
\end{equation}
for the ten-site Hubbard model with $U/t = 2$. Here, the expectation values are computed using linear-response density matrices and orbital relaxation effects have not been considered. As one would expect, at very high temperature $\theta$, there is no correlation between adjacent spins. As $\theta$ is reduced, the correlation appears and becomes maximal in the zero-temperature limit. Again, as with the internal energies, we see that the correlation function approaches the ground state CC in the limit $\beta \rightarrow \infty$.

\section{Conclusions}

We have demonstrated that the framework of thermofield dynamics can be exploited to formulate a finite-temperature coupled cluster theory.
We use the CCSD approximation to benchmark our method on various many-electron Hamiltonians and find that it performs substantially better than the thermal Hartree-Fock and thermal CI, just as one would expect for their ground-state counterparts.
This improved accuracy comes with the same asymptotic $\mathcal{O}(N^6)$ scaling ($N$ being the number of spin-orbitals or basis functions) as does standard quasiparticle CCSD, with a modestly larger pre-factor, though note that we must solve these equations at each grid point in the imaginary-time evolution.
We also observe that in the zero-temperature limit, thermal HF and CCSD approaches the appropriate ground-state HF and CCSD.

The thermofield based formulation of CC provides the basis for generalization of more sophisticated CC ansatze, e.g. unitary and multi-reference CC, to finite-temperatures, something that will be explored in future work.
Moreover, while CCSD with a perturbative treatment of triple excitations, or CCSD(T), has claimed the status of ``gold standard'' among ground-state methods for achieving chemical accuracy in weakly correlated systems, it is not clear how a similar notion can be defined for finite-temperature methods. More benchmark studies need to be performed to understand and improve accuracy of thermal CC.

\begin{acknowledgments}
        This work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Computational and Theoretical Chemistry Program under Award No. DE-FG02-09ER16053. G.E.S. acknowledges support as a Welch Foundation Chair (No. C-0036).
\end{acknowledgments}


\section{Introduction}
In numerous circumstances, more than just the content of a message has to be hidden from the adversary.
Unlike covertness which aims to deny that any communication is taking place \cite{DeniableComm:KadheJB14, HidingInformationInNoise:BashGT15, CovertComm:Bloch16}, we consider the case where it is known that a group of peers communicate but it is desired to hide who is communicating with whom \cite{TorlikeAnonymityInWirelessNets:SafakaCA15}.
It is well known that identities of peers communicating over a network can be identified via rather simple network traffic analysis techniques \cite{NetsWithUserObservability:PfitzmannW85}.
Anonymity mixes were introduced by David Chaum in 1980's as a general framework for implementing anonymous message exchange \cite{Mixes:Chaum81}.
They are sophisticated network routers that pass messages such that no one (except the mix itself) can link an ingoing message to an outgoing message. Today, some form of a mix is often a part of anonymity preserving solutions (e.g., PetMail, Mixminion, Panoramix) or data transfer services (e.g., Onion routing, Freenet).

A mix typically collects messages and forwards them in batches according to a fixed deterministic rule or a randomized strategy (see e.g., \cite{BatchMix:KesdoganAP02, DummyMessageMix:BertholdL02, Mix:diazS03}). This allows hiding the origin of the outgoing messages, but also introduces delay in message transfer.
The incurred delay of the mixes is the most concerning cost of anonymity they provide.
For instance anonymous web browsing platform ToR, which currently has more than 2 million users, does not implement sophisticated mixing to keep a low latency platform, even though it is shown to be vulnerable to deanonymization attacks based on network traffic analysis \cite{ChallengesDeployingTor:DingledineMS05, Tor:dingledineMS07,TorTraffAnalysis:ChakravartyBP14}.

Appropriate modeling of the mixes is crucial to study their delay vs.\ anonymity tradeoff.
Stability conditions and delay characteristics of a mix naturally depend on its system parameters which also determine its anonymity preserving capabilities.
In this paper, we propose and study two queueing models for batch mixes that are designed and used against passive adversarial attacks.
Note that, we do not consider active attacks that involve traffic injection into the network, which have also been shown to successfully deanonymize users on popular anonymity platforms \cite{SampledTrafficAnalysis:MurdochZ07, TorTraffAnalysis:ChakravartyBP14}.

We propose a mix model that implements the well known deterministic \textit{batch} mixing algorithm \cite{BatchMix:KesdoganAP02}.
We observe the close connection of the model to assembly queues, which was used to model and study the operational process of assembling multiple items into a product \cite{AssemblyQs:Harrison73}.
Using the proposed model, we find that batch mix provides a well defined anonymity guarantee that gets better with the batch size, on the other hand, its incurred delay grows quickly as the batch size gets close to the number of senders connected to the mix.

Our study of the batch mix led us to consider a new randomized mixing algorithm.
We show that the randomized model achieves better delay scaling in terms of the batch size compared to its deterministic counterpart. However, it can provably preserve anonymity \emph{only if} the adversary can not infer the state of the mix, and is in general vulnerable to anonymity attacks under low traffic.

There are many measures of anonymity and privacy (see \cite{Unlinkability:PfitzmannH05, InfoLeakage:IssaKW16}). We are here concerned with preserving \emph{unlinkability}, which is ensuring that no sender/receiver pair is exposed. Our study shows that delay of the mix can be reduced by sacrificing some anonymity, which would eventually lead to complete deanonymization of all the sender-receiver pairs.
However,
message transfer sessions are of finite duration in practice, and minimum amount of time required for an attack to destroy anonymity is a concern regardless of the anonymity measure.
Previous papers that are concerned with the delay of anonymity schemes ignore the queueing dynamics within the mix (see e.g., \cite{DelayAnonymityTradeoffInMixNetworks:JavidbakhtV17, AnonymousNetworkingWithMinimumLatency:VenkitasubramaniamT08} and references therein).
We believe that queueing models are necessary for studying the time-to-deanonymize metric, and this paper is a first step towards understanding this metric.

This paper is organized as follows.
Sec.~\ref{sec:sec_batchmix_model} describes the batch mix model.
Sec.~\ref{sec:sec_batchmix_anonymity} presents the anonymity guarantee implemented by a batch mix.
Sec.~\ref{sec:sec_batchmix_stability_delay} shows a stability criterion for the batch mix and presents an approximate method for analyzing its incurred delay.
Sec.~\ref{sec:sec_samplingmix} introduces a randomized batch mixing strategy, and discusses its anonymity and delay properties. Sec.~\ref{sec:sec_conclusion} gives a summary and conclusions.

\section{A Batch Mix  and its Anonymity}
\subsection{Mix Model}
\label{sec:sec_batchmix_model}
A batch mix has $n$ senders connected, and buffers the messages incoming from each sender in a separate first-in first-out queue with an infinite buffer space. As soon as any $k \geq 2$ queues become non-empty, one message from each is dispatched (see Fig.~\ref{fig:fig_43_mix}).
The recipient sets of each sender are assumed to be disjoint and of at most size $m$.

\begin{figure}[t]
  \centering
  \includegraphics[width=0.4\textwidth, keepaspectratio=true]{fig_43_mix.pdf}
  \caption{Illustration of a $(4, 3)$ batch mix. As long as there are less than three non-empty queues, messages are blocked (Left). As soon as a message arrival forms a group of three non-empty queues, one message from each is dispatched (Right).}
  \label{fig:fig_43_mix}
\end{figure}

Each sender is assumed to generate an independent Poisson
message traffic at rate $\lambda$.
Delay added by the mix is assumed to come only from the message queueing time. We ignore any message reception or transmission delay.

\begin{definition}
  An \ul{$(n, k)$ batch mix} is a system of $n$ first-in first-out queues, each receiving messages from an independent Poisson process of rate $\lambda$.
  Messages are blocked as long as the mix has less than $k$ non-empty queues. As soon as $k$ queues become non-empty, one message from each is dispatched.
\end{definition}

\subsection{Attack Model and Anonymity}

\label{sec:sec_batchmix_anonymity}
 We assume that the adversary monitoring the traffic going in and out of the mix can observe 1) who the sender of each incoming message is and 2) who the recipient of each outgoing message is. Thus, if a message arrival triggers a message departure, the adversary can identify the sender-receiver pair. His goal is to identify the receivers of a particular sender, which we refer to as the \textit{target sender}.

Forwarding messages in batches of size $k$ prevents the adversary from immediately finding out the exact destination of an incoming message, as it can be any of the $k$ message recipients. However, the adversary can, over time, collect multiple size-$k$ receiver sets, each containing a potential recipient of the target sender. Intersecting such sets would eventually reveal the receivers linked to the target sender. We refer to such attacks as \textit{intersection attacks} \cite{IntersectionAttack:WrightAL02}.

We say that the mix preserves anonymity, when it ensures that no sender/receiver pair is exposed, that is, no sender and receiver can be linked.

\begin{theorem}[Anonymity under intersection attack]
  Consider a target sender connected to an $(n, k)$ batch mix that is under intersection attacks.
  When $k < n$, all $m$ receivers of the target can be identified if $m \leq \floor*{n/k}$.
  All $m$ receivers cannot be identified surely otherwise.
\label{thm_batchmix_anonymity}
\end{theorem}
\begin{proof}
  This theorem is a reformulation of \cite[Claim~1]{BatchMix:KesdoganAP02}.
  Let adversary wait and observe $m$ \emph{mutually disjoint} sets $R_1, \ldots, R_m$ of size $k$ that include the possible receivers of Alice. These $m$ sets can be disjoint only if $k m \leq n$. Adversary is thus sure that there is exactly one receiver of Alice in each observed recipient set $R_i$. Afterwards, adversary refines each of these sets by observing new recipient sets that intersect with only one of the prior sets. This means, a new recipient set $R$ is useful if $R \cap R_i \neq \emptyset$ and $R \cap R_j = \emptyset$ for all $j \neq i$, then $R_i$ can be refined to $R \cap R_i$.
  Note that if $R$ intersects with multiple prior recipient sets, then refining all intersecting sets  may remove the actual receivers of Alice.
  The correct refinement process can be continued until each of the sets $R_1, \ldots, R_m$ contains only one receiver. Remaining $m$ receivers in the refined recipient sets are clearly the communication partners of Alice.

  As described above, intersection attacks will surely identify all receivers of a target only if adversary can observe $m$ \emph{disjoint} sets of size $k$. This is the only way for adversary to isolate each receiver of the target in a separate set so that a newly observed set can be intersected with one of these sets \emph{correctly}, that is, intersection will not surely end up removing the true receiver from the set.
  When $k m>n$, adversary can never observe $m$ disjoint sets of size $k$, hence can never surely identify all $m$ receivers of the target.
\end{proof}

\noindent
\section{Stability and Delay}
\label{sec:sec_batchmix_stability_delay}
A batch mix consists of $n$ FIFO queues, each buffering messages arriving from an i.i.d.\ Poisson process. A message arrival triggers a batch departure if it finds $k-1$ other non-empty queues in the mix, and the arriving message departs immediately with the batch. Therefore, there can be at most $k-1$ non-empty queues in the mix at any time. Since all the queues and the associated arrival processes are identical, system state can be represented as the Markov process $L(t) = (l_1(t), \dots, l_{k-1}(t))$ where $l_i(t)$ denotes the length of the $i$th longest queue in the system at time $t$.

An $(n, n)$ batch mix behaves as an assembly queue, found  to be unstable in \cite{AssemblyQs:Harrison73}.
Stability here refers to the existence of an invariant probability measure for the system state process.

\begin{theorem}
  An $(n, k)$ batch mix is stable if $k < n$.
\label{thm_nkmix_isstable}
\end{theorem}
\begin{proof}
  A Markov process is stable if and only if it is positive recurrent.
  Given that transition rates of $L(t)$ are neither too ``slow'' nor too ``fast'', its positive recurrence is implied by the positive recurrence of its embedded discrete chain $S_t$.
  We here use the Foster-Lyapunov criterion to show the positive recurrence of $S_t$ as interpreted from \cite[Thm.~2]{MarkovChainStability:FossK04}.

  For system state $s = (s_1, \dots, s_{k-1})$, let
  \[ W(s) \coloneqq s_{k-1}^{\log_2(n/n-1)}. \]
  Recall that $s_{k-1} = \min\{s_i, ~i=1, \dots, k-1\}$.

  Note that $\sup_s W(s) = \infty$ as required.
  One step drift for any state $s \in \{s, ~W(s) > 0\}$ is
  \[ \mathbb{E}\left[W(S_1) - W(S_0) \mid S_0=s\right] < 0. \]
  and we have
  \[ \sup_{\{s,\, W(s) \leq 0\}} \mathbb{E}[S_1 \mid S_0=s] < 1 < \infty. \]
  Thus $S_t$, hence $L(t)$ is positive recurrent.

\end{proof}

\begin{figure*}[t]
  \centering
  \begin{subfigure}[]{.32\textwidth}
    \centering
    \includegraphics[width=1\textwidth, keepaspectratio=true]{plot_nkmix_n40_varyingk.pdf}
  \end{subfigure}
  \begin{subfigure}[]{.32\textwidth}
    \centering
    \includegraphics[width=1\textwidth, keepaspectratio=true]{plot_nkmix_n40_varyingar.pdf}
  \end{subfigure}
  \begin{subfigure}[]{.32\textwidth}
    \centering
    \includegraphics[width=1\textwidth, keepaspectratio=true]{plot_nkmix_vs_randmix_n40.pdf}
  \end{subfigure}
  \caption{(Left, Middle) Average delay in a batch mix with $n=40$; Left: fixed $\lambda$ and varying $k$, Middle: fixed $k$ and varying $\lambda$. (Right) Comparison between the average delay in a batch mix and in a sampling mix with $p_a = 1/n$.}
  \label{fig:plot_ED_sim_vs_approx}
\end{figure*}

\noindent
There are three scenarios that a message can experience upon arrival to the mix:
1) If a message arrives to an empty queue and finds $k-1$ other non-empty queues in the mix, then it will immediately depart with no queueing.
2) If a message arrives to an empty queue and finds fewer than $k-1$ other non-empty queues in the mix, then it has to wait for a formation of $k$ non-empty queues (i.e., a batch).
3) If a message arrives to a non-empty queue, it has to first wait to the HoL (head of the line) in its queue, and then wait for the next batch formation.

In a tagged queue, batch formation delay experienced by a message is completely characterized by the number of non-empty queues $R$ (excluding the tagged queue) seen by the message once it moves to HoL. If $R < k-1$, message will be blocked until any $k-1-R$ of the $n-1-R$ empty queues receive at least a message.
Using the memoryless property of message inter-arrival times, batch formation delay is distributed as the $(k-1-R)$th order statistic of $n-1-R$ i.i.d.\ exponentials, which we denote as $X_{n-1-R:k-1-R}$\footnote{$X_{i:j} \coloneqq 0$ if $i < j$ or $j = 0$.}.
Overall, a message moving to HoL may find from $0$ up to $k-1$ other non-empty queues, hence there are $k$ possible different distributions for the batch formation delay.

When $k=2$, system state is just the longest queue length and defines a birth-death process. Exact analysis is formidable when $k > 2$ because of the infamous state space explosion problem. We first present the exact analysis for $k=2$, then present an approximate method for $k > 2$, which is similar to an approximation presented for assembly queues in \cite{AssemblyQs:LipperS86}.

\subsection{Exact analysis of $(n, 2)$-mix}
\label{subsec:subsec_n_2}
In $(n, 2)$-mix for $n > 2$, there can be at most one non-empty queue at any time, hence the system state is captured by the length of the longest queue $\mathcal{L}(t)$. It defines a single dimensional birth-death Markov process as shown in Fig.~\ref{fig:fig_n_2_markovprocess}.
\begin{figure}[t]
  \centering
  \includegraphics[width=0.35\textwidth, keepaspectratio=true]{fig_n_2_markovprocess.pdf}
  \caption{Markov process for $(n, 2)$-mix. State here denotes the length of the longest queue in the mix; $0$ length means the system is empty.}
  \label{fig:fig_n_2_markovprocess}
\end{figure}

Exact analysis in this case is straightforward. Let $p_l$ be the stationary probability for state $l$. From global balance equations we find
\[ p_0 = \frac{n-2}{2(n-1)} \quad\quad p_l = \frac{n(n-2)}{2(n-1)^{i+1}};~l=1,2,\ldots \]
Ergodicity implies that fraction of the time an arbitrary queue is non-empty (i.e., average load on the queue) is $\rho = (1-p_0)/n = 1/2(n-2)$. Larger $n$ gives higher frequency of emptiness at the servers, which is natural since queues empty out faster when the mix receives messages at a higher rate.

Using the stationary state probabilities, first two moments of the length of an arbitrary queue are given as
\[
\begin{split}
  \mathbb{E}[L] &= \frac{1}{n}\sum_{l=1}^{\infty} l\; p_l = \frac{1}{2(n-2)} \\
  \mathbb{E}[L^2] &= \frac{1}{n}\sum_{l=1}^{\infty} l^2\; p_l = \frac{n}{2(n-2)^2}
\end{split}
\]

We next derive some simple conclusions for the steady state delay experienced by an arriving message. Using PASTA \cite{Pasta:Wolff82}, an arbitrary message finds the system empty with probability $p_0$ and will have to wait for the first arrival to one of the other $n-1$ queues. Since arrivals are Poisson, waiting time distribution for the message is minimum of $n-1$ $\mathrm{Exp}(\lambda)$'s, that is $\mathrm{Exp}\left((n-1)\lambda\right)$. An arriving message may find its corresponding queue with $l$ messages with probability $p_l/n$. In this case, waiting time distribution for the message is sum of $l+1$ independent $\mathrm{Exp}\left((n-1)\lambda\right)$'s, which is $\mathrm{Erlang}\left(l+1, (n-1)\lambda\right)$. Finally, it may also find its corresponding queue empty with probability $p_l(n-1)/n$ if there is another queue with $l$ messages. Then the message will not be queued and will depart immediately upon arrival together with the first message in the busy queue. Using the law of total probability, distribution of waiting time $D$ for an arbitrary message is then given as
\begin{equation*}
\begin{split}
  \Pr\{D > w\} &= p_0\; \Pr\left\{\mathrm{Exp}\left((n-1)\lambda\right) > w\right\} \\
  &+ \frac{1}{n}\sum_{l=1}^{\infty} p_l\; \Pr\left\{\mathrm{Erlang}\left(l+1, (n-1)\lambda\right) > w\right\}.
\end{split}
\end{equation*}

\subsection{Approximate analysis of $(n, k>2)$-mix}
\label{subsec:subsec_n>2}
We here adopt the following approximating assumption;
a message \emph{upon moving to head of the line (HoL)} in its queue finds each other queue independently non-empty with probability $p$. Given that, and the fact that there can be at most $k-1$ non-empty queues at any time, the number of non-empty queues seen by a message moving to HoL is distributed as $R \sim B | \{B \leq k-1\}$ where $B \sim \mathrm{Binom}(n-1, p)$.
Given $R = r$, message will have to wait before getting dispatched for the first $k-1-r$ among all the $n-1-r$ empty queues to receive at least one arrival, that is, the message will experience a batch formation delay of $V|\{R=r\} \sim X_{n-1-r:k-1-r}$.
Then $V$ for an arbitrary message, which \emph{arrives to a non-empty queue in the first place}, is approximately distributed as
\begin{equation}
  \Pr\{V \leq v\} = \mathbb{E}_R\left[\Pr\{X_{n-1-R:k-1-R} \leq v\}\right],
  \label{eq:eq_MG1efs_Vcdf}
\end{equation}
where $\mathbb{E}_R$ denotes expectation with respect to $R$.

Batch formation delay for messages that arrive to an empty queue is differently distributed (than $V$ above) because they find each other queue non-empty with a different probability than messages that arrive to a non-empty queue.
Let a tagged queue be left empty by a departing message $m$. If the queue was left non-empty, the next message in line would have immediately moved to HoL. Then according to our earlier assumption, the number of non-empty queues left behind non-empty by $m$ is distributed as $R$.
Including the next arrival to the tagged queue, say message $m^+$, the next batch formation requires $k-R$ of the $n-R$ empty queues to receive at least an arrival.
Given that $m^+$ is among these first $k-R$ arrivals and messages are generated from i.i.d. streams, probability that $m^+$ is the $j$th among the $k-R$ arrivals is $1/(k-R)$.
Thus, batch formation delay experienced by a message arriving to an empty queue is approximately distributed as
\begin{equation}
  \Pr\{V_e \leq v\} = \mathbb{E}_R\left[\sum_{j=1}^{k-R} \frac{\Pr\{X_{n-R-j:k-R-j} \leq v\}}{k-R} \right].
\label{eq:eq_MG1efs_Vecdf}
\end{equation}

\begin{proposition}[Approximation by decoupling the queues]
  An $(n, k)$ batch mix approximately behaves to each sender as an $M/G/1/\text{efs}$ queue with regular service times distributed as \eqref{eq:eq_MG1efs_Vcdf} and exceptional first service times distributed as \eqref{eq:eq_MG1efs_Vecdf}.
\end{proposition}

Approximation requires estimating $p$, for which a natural estimate would be the average load $\rho$ for a queue, which is known for an $M/G/1/\text{efs}$ queue to be \cite{MG1ExceptionalFirstBusy:Welch64}
\[ \rho = \frac{\lambda \mathbb{E}[V_e]}{1 - \lambda\left(\mathbb{E}[V] - \mathbb{E}[V_e]\right)}. \]
Moments of $V$ and $V_e$ depend on $p$, hence on its estimate $\rho$. The equality above can be solved numerically to find a value for $\rho$.
Simulated and approximated values of delay are compared in Fig.~\ref{fig:plot_ED_sim_vs_approx} for a $(40, k)$ batch mix. Despite the strong independence assumptions employed in deriving the approximation, it compares well with the simulations for low values of $k$, which is the practically relevant case since the incurred delay must be kept below a threshold in practice.

\section{Sampling mix}
\label{sec:sec_samplingmix}
A sampling mix also implements an $(n, k)$ model; mix buffers the messages from each of the $n$ connected senders in a separate FIFO queue and each sender communicates with a disjoint set of at most $m$ receivers. However, buffered messages are forwarded differently compared to batch mix; as soon as a message arrives to the mix, $k$ queues are randomly selected and released. Releasing a queue allows it to forward a message if it is non-empty.
Specifically with probability $p_a$, the queue that receives the arrival is selected together with $k-1$ queues chosen uniformly at random from the remaining $n-1$ queues, or with probability $1-p_a$, the queue that receives the arrival is skipped and $k$ queues are selected uniformly at random from the rest of the queues.

\begin{theorem}
  \underline{Average load} of a queue in an $(n, k)$ sampling mix is given as
  \[
  \rho = (1 - p_a)/(k - p_a)
  \]
  and the \underline{average delay} experienced by a message is given as
  \[
  (1 - p_a)/\lambda(k-1).
  \]
\end{theorem}
\begin{proof}
  The length of a particular queue in the mix defines a birth-death process with a state space of non-negative integers and transition rates given for $i \geq 0$ as
  \[ \Pr\{i \to i+1\} = \lambda(1 - p_a), ~~ \Pr\{i+1 \to i\} = (n-1)\lambda p_o, \]
  where $p_o = p_a(k-1)/(n-1) + (1-p_a)k/(n-1)$.
  Stationary state probabilities are easily derived, using which average length of a queue is found, then Little's law is applied.
\end{proof}

As shown in Fig.~\ref{fig:plot_ED_sim_vs_approx}, average delay of a sampling mix scales much better with $k$ (i.e., decays as $1/(k-1)$) compared to a batch mix (i.e., grows exponentially in $k$ beyond a value).
However, a sampling mix cannot provide a well-defined anonymity guarantee while a batch mix can (see Thm.~\ref{thm_batchmix_anonymity}).

\begin{theorem}
  All receivers of a target sender connected to a sampling mix can be identified with intersection attacks by an adversary that can infer the state of the mix.
\label{thm_samplingmix_statefulattack}
\end{theorem}
\begin{proof}
  Queues in the mix will empty out infinite number of times under stability. Suppose that the adversary can detect whenever the mix becomes empty.
  Firstly, assume $p \neq 0$. Given that a message from a target finds the mix empty, the arriving message will be forwarded with probability $p$ or no message will depart. If the message is immediately forwarded, a receiver of the target will revealed.
  Number of times repeating this attack required to identify a receiver is geometric with $p$, hence attack will be almost surely successful in finite time.

  Secondly, assume $p \neq 0$. Given that a message from a target finds more than one non-empty queue in the mix, the following departure may include messages going only to the receivers of the non-target senders. This reveals which receiver does not belong to the recipient set of the target. Eventually adversary will be left with the correct set of receivers.
\end{proof}

Sampling mix will empty out more frequently and cannot often build a state complex enough to hide the origin of the outgoing messages when $k$ is larger and/or arrival rate $\lambda$ is lower, hence intersection attacks with state knowledge will resolve faster.
Moreover, even simple intersection attacks that do not require state knowledge can deanonymize a target connected to a sampling mix if $p_a$ is not chosen carefully.

\begin{theorem}
  All receivers of a target sender connected to an $(n, k)$ sampling mix with $p_a \neq 1/n$ can be identified with intersection attacks that do not require state knowledge.
\label{thm_samplingmix_statelessattack}
\end{theorem}
\begin{proof}
  Once a message arrives to a queue in steady state, probability of a departure from any other queue is $q = p_o\rho = (1-p_a)/(n-1)$.
  Suppose $m=1$ and $p_a > q$ ($p_a < q$). Adversary can record the message departures per arrival from a target sender.
  By the law of large numbers, the greatest (smallest) number of departures will almost surely be observed on the correct receiver in the limit.
  Same idea applies when each sender communicates with multiple receivers.
  Finally, $p_a = q$ if and only if $p_a = 1/n$.

  In other words, in order to preserve anonymity, it is necessary to maximize  the uncertainty within the steady state probabilities of message departures from the queues.
  R\'enyi entropy is commonly used for measuring uncertainty and uniform distribution is known to maximize it, which is achieved by setting $p_a = 1/n$.
\end{proof}

\section{Conclusion}
\label{sec:sec_conclusion}
We proposed a queueing model for batch anonymity mixes and showed that batch a mix with a deterministic message dispatching policy ensures that no sender-receiver pair is exposed (referred to as anonymity in this paper) under intersection attacks. On the other hand, its incurred delay on message transfer grows quickly as the batch size gets close to the number of connected senders.
We introduced a sampling mix model that implements a randomized message dispatching policy. Sampling mix permits an exact delay analysis, which allowed us to show that randomization allows cutting the tail of delay immensely, however, at the cost of giving up on the anonymity guarantee implemented by its deterministic counterpart.
We hope to next use our proposed queueing model to understand the performance of mixes in terms of the time-to-deanonymize metric vs.\ the incurred message transfer delay.

\section*{Acknowledgments}
This research is based upon work supported by the National Science Foundation
under Grant No.~SaTC-1816404.

\bibliographystyle{IEEEtran}


\section{Introduction}

Since the pioneering work of Gage and Hamilton \cite{GageH} in 1986, the curve shortening problem has been studied widely in different settings. An important case is to investigate the problem with some boundary conditions.  Those kind of boundary problems are not only interesting from geometric viewpoint, but also important in the study of phase transition models. Also, they can be used to describe the interface propagation arising in some reaction-diffusion equations, such as the Allen-Cahn equation. Usually, there will be a contact angle condition
at the intersection points between the evolving curves (interfaces) and the domain boundaries (see, for example, \cite{ORS, RSK}). In particular, we focus on the following curvature flow
\begin{equation}
  V=H \qquad \text{on} \qquad\Gamma_t\subset \Omega, \label{mcf}
\end{equation}
in the band domain $\Omega:=\{(x,y)|-1\leq x\leq 1, y\in\R\}$ in $\R^2$, where, $\Gamma_t$ is a family of embedded curves
in $\Omega$ which contact the boundaries $\partial_\pm \Omega := \{\pm1\}\times\R $ of $\Omega$ with prescribed angles, $V$ and $H$ denote the normal velocity and the curvature of $\Gamma_t$, respectively. Indeed, Huisken \cite{Hui} considered this curvature flow in 1998. He discovered a distance comparison principle to show that the flow with fixed endpoints (Dirichlet boundary condition) exists globally and converges to a straight line linking two endpoints, and the flow having orthonormal boundary contact (homogeneous Neumann boundary condition) also exists globally but converges to a constant eventually.

In case, for each $t\geq 0$, $\Gamma_t$ is the graph of a function $y=u(x,t)$, we have
$$
  V=\frac{u_t}{\sqrt{1+u^2_x}},\qquad H=\frac{u_{xx}}{\left(1+u^2_x\right)^{3/2}},
$$
and the problem can be expressed as
\begin{equation}\label{A}
\left\{
 \begin{array}{ll}
\displaystyle u_t=\frac{u_{xx}}{1+u^2_x}, & -1<x<1,\ t>0,\\
u_x(\pm 1,t)=g_\pm ,& t>0,
 \end{array}
 \right.
\end{equation}
where $g_-, g_+$ denote the boundary contact conditions and will be made clear later.
 In 1993, Altschuler and Wu \cite{AW1} studied the inhomogeneous Neumann
problem, that is, the case where $g_+, -g_-$ are positive constants. They
proved that any global solution converges to a {\it grim reaper} (which is also called a
traveling wave by some authors).  In 1994, Altschular and Wu  \cite{AW2} extended their result
in \cite{AW1} to two dimension. In 2012, Cai and Lou \cite{CaiLou2} considered \eqref{A}
with $g_\pm $ being (almost) periodic functions of $u$, and proved that any solution
converges to a (almost) periodic traveling wave. Recently, Yuan and Lou \cite{YuanLou} considered the case
where $g_\pm = g_\pm (u)$ are asymptotic periodic functions as $u\to \pm \infty$.
They constructed some entire solutions connecting two periodic traveling waves.

Other works related to the mean curvature flow \eqref{mcf}, as well as
its anisotropic analogues, in domains in $\R^2$ with boundaries include \cite{LMN, MNL} for problems in band domain with undulating boundaries;
\cite{CGK, GGH, GH, Koh} for self-similar solutions in sectors; \cite{ChenGuo, GMSW}
for problems on the half plane, etc.

\medskip

In all of the above mentioned works, the boundary slopes are bounded, no matter they are linear or nonlinear.
In this paper, we consider a different case:
\begin{equation}\label{B}
\left\{
 \begin{array}{ll}
\displaystyle u_t=\frac{u_{xx}}{1+u^2_x}, &   -1<x<1,\ t>0,\\
u_x(\pm1,t)=\pm u(\pm1,t), &  t>0,\\
u(x,0)=u_0(x), & -\leq x\leq 1,
 \end{array}
 \right.
\end{equation}
which has {\it unbounded boundary slopes} when $u$ goes to infinity. This kind of boundary conditions are of
interest not only in geometry and partial differential equations but also in some applied fields.
In 2012, Chou and Wang \cite{CW} considered the equation in \eqref{B} with Robin boundary conditions:
$$
u_x(\pm 1,t) = \alpha_\pm u(\pm 1,t) +\beta_\pm,\quad t>0.
$$
They divided the parameters $\alpha_\pm$ and $\beta_\pm$ into several cases, and
studied the asymptotic behavior in each cases. For the cases $\min u\to \infty$
or $\max u\to -\infty$ (this happens, for example, when $\alpha_- <0 <\alpha_+$), however,
the authors did not obtain the convergence of the solution and left it as an open problem.
In this paper, we will show the convergence to a {\it grim reaper} for unbounded solutions to \eqref{B}.
Our approach (with slight modification) and our results remain valid for the problems with
other unbounded boundary slopes such as the Robin ones, $u_x (\pm 1,t) = \pm u^2 (\pm 1, t)$,
$u_x(\pm 1,t) = \pm e^{u(\pm 1,t)}$, etc.

The global well-posedness of the problem \eqref{B} is studied in a standard way.
For any time-global solution moving upward to infinity, $u_x$ is unbounded due to the boundary conditions.
This will be the main difficulty in studying the asymptotic behavior for such solutions. It
is natural to consider the convergence in $L^\infty_{loc} ((-1,1))$ topology in such cases.
This, however, also needs some {\it uniform} (in $t$) interior gradient
estimates.  The well known results in this field as in \cite{ES} are not applied here since
they depend on the boundedness of $u$. Instead, we will use the so-called zero number argument
(i.e., zero number diminishing properties, cf. \cite{Ang, Lou2}) for one dimensional parabolic equations to derive the uniform
bounds for the gradient of the solution in any interior domain (see details in sections 4 and 5).
Furthermore, as can be expected, the profile of the solution might converge to a traveling wave with
infinite slopes near the boundaries, which should be the {\it grim reaper} with span in $(-1,1)$, that is,
\begin{equation}\label{0-gr}
\varphi_0 (x) + \frac{\pi}{2} t \quad \mbox{ with } \quad \varphi_0 (x) := -\frac{2}{\pi} \ln \left[\cos \Big( \frac{\pi x}{2} \Big)\right]\ \ (x\in (-1,1)).
\end{equation}
Actually we will show the following result.

\begin{thm}\label{thm:main}
Assume $u_0(x)\in C^1([-1,1])$ with
\begin{equation}
  u_0(x)\geq 1\ (x\in [-1,1]),\qquad u'_0(\pm1)=\pm u_0(\pm 1) .\label{1.4}
\end{equation}
Then the problem \eqref{B} has a time-global solution $u(x,t)$. It moves upward to
infinity and
\begin{equation}
  u(x,t+s) - u(0,s) \to \varphi_0(x) + \frac{\pi}{2} t, \qquad \mbox{as}\ \  s\rightarrow \infty,\label{5}
\end{equation}
in the topology of $C^{2,1}_{loc}\left( (-1,1)\times\R\right).$
\end{thm}

This paper is arranged as follows. In section 2, as preliminaries, we present traveling waves ({\it grim reaper}s)
contacting the boundaries of $\Omega$ with various different constant slopes. In section 3, we give some rough a
priori estimates and then show the time-global existence for the solution of \eqref{B}.
In section 4, we consider symmetric solutions of \eqref{B}. First we present
precise estimates for $u_x$ by using the zero-number argument, and then show its convergence to the {\it grim reaper}.
Finally,  in section 5 we consider general solutions of \eqref{B} which are not necessary symmetric.
By some further uniform interior estimates we show its convergence to the {\it grim reaper}.

\section{Traveling Waves}

First we give a definition.
\begin{defn}
A function $\underline{u}(x,t)$ satisfying
$$
\left\{
 \begin{array}{l}
\displaystyle \underline{u}_t\leq\frac{\underline{u}_{xx}}{1+\underline{u}^2_{x}},\qquad -1<x<1,\ t>0,\\
\underline{u}_x(1,t)\leq \underline{u}(1,t),\quad\underline{u}_x(-1,t)\geq -\underline{u}(-1,t),\qquad t>0,
 \end{array}
 \right.
$$
is called a lower solution of (\ref{B}). A function $\overline{u}(x,t)$ satisfying the reversed inequalities is called an upper solution
of \eqref{B}.
\end{defn}

As preliminaries, we study the following problem
\begin{equation}\label{C}
\left\{
 \begin{array}{ll}
\displaystyle u_t=\frac{u_{xx}}{1+u^2_x}, & -1<x<1,\quad t\in\R,\\
u_x(\pm1,t)=\pm h, & t\in \R.
 \end{array}
 \right.
\end{equation}
For each $h>0$, a traveling wave of (\ref{C}) (also called a translating solution in \cite{AW1})
is a special solution of the form
$$
u(x,t)=\varphi(x;h)+c(h)t.
$$
Substituting this formula into (\ref{C}) we easily obtain
\begin{equation}
\varphi(x;h) :=-\frac{1}{c(h)}\ln \left[\cos\left(c(h)x\right)\right].
\qquad c(h) :=\arctan h,
\end{equation}
$\varphi$ is called a {\it grim reaper} in \cite{AW1}.  Note that
\begin{equation}
\varphi(\pm1;h)= \frac{\ln(1+h^2)}{2\arctan h},\qquad \varphi_x(\pm1;h)=\pm h.
\end{equation}
Hence, for any $M\in\R$, $\varphi(x;h)+c(h)t+M$ is a lower solution of (\ref{B}) when
$$
h\leq\varphi(1;h)+c(h)t+M.
$$
It is an upper solution if the reversed inequality holds.

Besides the traveling waves of \eqref{C}, we have another {\it grim reaper} $\varphi_0(x)+\frac{\pi}{2}t$
with $\varphi_0$ defined by \eqref{0-gr}. Note that the definition domain of $\varphi_0(x)$ is $(-1,1)$,
that is, this {\it grim reaper} lies completely in $\Omega$. In what follows we will use the above
{\it grim reaper}s to give a priori estimates for the solution of (\ref{B}).

\section{Global Well-posedness of (\ref{B})}

Assume $u(x,t)$ is a classical solution of (\ref{B}) in the time-interval $[0,T]$ for some $T>0$.
We first give its $L^{\infty}$ estimate.
\begin{lemma}\label{bound-est}
 There exist $C_1,C_2>0$ with $C_2$ depending on $T$ such that
 \begin{equation}\label{3.1}
 c(1) t -C_1 \leq u(x,t)\leq C_2(T),\qquad x\in[-1,1],\ t\in[0,T].
 \end{equation}
\end{lemma}

\begin{proof}
Assume, for some $M_0>0$,
\begin{equation}\label{3.2}
1\leq u_0(x)\leq M_0,\qquad x\in[-1,1].
\end{equation}
Then $$
\underline{u}(x,t):=\varphi(x;1)+c(1)t+1-\varphi(1;1)
$$
is a lower solution of (\ref{B}) and so, by the comparison principle we have
\begin{equation}\label{3.3}
\underline{u}(x,t)=\varphi(x;1)+c(1)t+1-\varphi(1;1)\leq u(x,t),\qquad x\in[-1,1],\  t>0,
\end{equation}
which leads to the first inequality in (\ref{3.1}).

Next we consider the upper bound. Since
$$
h-\varphi(\pm1;h)=h-\frac{\ln(1+h^2)}{2\arctan h}\rightarrow \infty\qquad \text{as} \qquad h\rightarrow \infty,
$$
there exists $h=h_T$ large such that
\begin{equation}\label{3.4}
h_T>\varphi(\pm1;h_T)+\frac{\pi}{2}T+M_0.
\end{equation}
Set
$$
\overline{u}(x,t):=\varphi(x;h_T)+c(h_T)t+M_0, \qquad x\in[-1,1],\ t>0.
$$
We verify that $\overline{u}$ is an upper solution of (\ref{B}) in the time interval $[0,T]$.
In fact, $\overline{u}$ satisfies the equation in (\ref{B}). Moreover, for $t\in[0,T]$,
$$
\overline{u}(1,t) <  \varphi(1;h_T)+\frac{\pi}{2}t + M_0<h_T = \overline{u}_x(1,t),
$$
by (\ref{3.4}), and
$$
\overline{u}(x,0)=\varphi(x;h_T)+M_0\geq M_0\geq u_0(x), \qquad x\in[-1,1].
$$
Hence, $\overline{u}$ is an upper solution of (\ref{B}). By comparison principle we have
$$
u(x,t)\leq \overline{u}(x,t)\leq \varphi(1;h_T)+M_0+\frac{\pi}{2}T,
\qquad x\in[-1,1],\ t\in[0,T].
$$
This proves the second inequality of (\ref{3.1}).
\end{proof}

\bigskip
Next, we give the gradient estimate.
\begin{lemma}
Let $u(x,t)$ be a solution of \eqref{B} in the time interval $[0,T]$. Then there exist $C_3(T)$ such that
$$
\left|u_x(x,t)\right|\leq C_3(T), \qquad x\in[-1,1],\ t\in[0,T].$$
\end{lemma}

\begin{proof}
From the above lemma, we see that
$$
\left|u_x(\pm1,t)\right|=\left|u(\pm1,t)\right| \leq C_1 + C_2(T) , \qquad t\in[0,T].
$$
Using the maximum principle for $u_x$ we see that
$$
\left|u_x(x,t)\right|\leq C_3(T):=\max\{\|u'_0\|_C, C_1 + C_2(T)\}, \qquad x\in[-1,1],\ t\in[0,T].
$$
This proves the lemma.
\end{proof}

With the above a priori estimates in hand, by using the standard parabolic theory
we obtain the time-global existence of the classical solution $u(x,t)$.
Its uniqueness is proved in the standard way by using the maximum principle.

\section{Symmetric Solutions}

In this section we consider symmetric solutions. More precisely, we consider the case where $u_0(x)\in C^1([-1,1])$
satisfies the following conditions:
\begin{equation}\label{ini-1}
  u_0(x)=u_0(-x),\quad 1\leq u_0(x)\leq M_0, \quad 0< u'_0(x)<\varphi'_0 (x) \mbox{\ for } x\in (0,1),
\quad u'_0 (1)= u_0(1),
\end{equation}
for some $M_0>0$.  In this case we easily have
$$
u(x,t)= u(-x,t),\quad u_x(x,t)>0 \quad \mbox{ for } x\in (0,1],\ t>0.
$$
To study the convergence of $u$ (actually, the convergence of $u(x,t+s)-u(0,s)$ as $s\to \infty$),
we need further estimate for $u_x$. We will do this by the so-called zero number argument.

\subsection{Finer lower gradient estimate}
In this part, we show that \lq\lq the gradient of $u$ is not too small" for $x\in (0,1)$.

Fix an $h_0>0$. For any $h\geq M_0+h_0$ we see that
$$
\underline{u} (x,t;h) := \varphi(x;h_0) +c(h_0) t +h
$$
is a lower solution of (\ref{B}). Denote the union of the graphes of $\underline{u}(x,t;h)$ as
$$
\mathcal{D}(t) := \left\{(x,y)\mid y\geq \underline{u} (x,t; M_0+h_0),\ x\in[-1,1] \right\}.
$$
Then $\mathcal{D}(t)$ is the upper half of $\Omega$ with bottom $\left\{\left(x, \underline{u} (x,t; M_0 +h_0)\right) \mid x\in[-1,1]\right\}$,
which moves upward with speed $c(h_0)$.

We now construct another lower solution below the real solution $u$ such that it moves faster than $\mathcal{D}(t)$, and so
pushes $u$ entering the domain $\mathcal{D}(t)$ for large $t$. Then, by considering the derivatives of
$u(x,t)$ and $\underline{u} (x,t;h)$ at their contact points, we can obtain the desired gradient estimate.
To construct another lower solution below $u$, we notice by (\ref{3.1}) that $u(x,t)$ moves up to infinity.
Hence for any given $h^0>h_0$, there exists $t^0$ large such that
$$
u(x,t^0)>\varphi(x; h^0)+h^0\geq h^0,\qquad x\in[-1,1].
$$
Thus,
$$
\underline{u}^* (x,t):=\varphi(x;h^0)+h^0+c(h^0)t
$$
is also a lower solution to (\ref{B}). By the comparison principle we have
\begin{equation}\label{4.1}
\underline{u}^* (x,t)\leq u(x,t+t^0), \qquad x\in[-1,1],\ t\geq 0.
\end{equation}

Since $c(h^0)>c(h_0)$ we see that for all large $t$ (to say, $t\geq T^0$), $\underline{u}^* (x,t)$ rushes
into the domain $\mathcal{D}(t)$, so dose $u(x,t+t^0)$ due to (\ref{4.1}).
In other words, when $t\geq t^0+T^0$, the graph of $u(x,t)$ is immersed in the solid $\mathcal{D}(t)$.
Therefore, for any $t_1 \geq t^0+T^0$ and any $x_1\in (0,1)$, there exists a unique $h\geq M_0+h_0$
such that $u(x_1,t_1) = \underline{u} (x_1,t_1;h)$. We now show that
\begin{equation}\label{underlineu-x<u-x}
0< \underline{u}_x (x_1,t_1;h) \leq u_x(x_1, t_1).
\end{equation}
In this sense, we say that \lq\lq the gradient of $u$ is not too small" for $x\in (0,1)$.

We first consider the case where $u(x,t_1) \leq \underline{u} (x,t_1;h)$ or
$u(x,t_1) \geq \underline{u} (x,t_1;h)$ for $x\in (0,1)$. In such cases, it is clear that
$\underline{u}_x (x_1,t_1;h) =u_x(x_1, t_1)$.
Next, we consider the case
\begin{equation}\label{max>0>min}
\max_{x\in [0,1]} [u(x,t_1)-\underline{u} (x,t_1;h)] >0 > \min_{x\in [0,1]} [u(x,t_1)-\underline{u} (x,t_1;h)].
\end{equation}
By continuity, these inequalities remain valid for any $t$ near $t_1$. Hence such $t$ form an interval.
Denote the largest interval containing $t_1$ by $(s_1, s_2)$, then \eqref{max>0>min} holds for any $t\in (s_1,s_2)$
and $s_1 >0$ since $\underline{u} (x,0;h) \geq h >M_0 \geq u(x,0)$. At $s_1$ we have the following claim:

\medskip
\noindent
{\bf Claim 1}. $\underline{u}(x,s_1;h)\geq u(x,s_1)$ in $[0,1]$, and $\underline{u}(1,s_1;h)= u(1,s_1)$.

Totally, there are four cases between $\underline{u}(x,s_1;h)$ and $u(x,s_1)$:
\begin{itemize}
\item[(a)] $\underline{u}(x,s_1;h) > u(x,s_1) $ in $[0,1]$;
\item[(b)] $\underline{u}(x,s_1;h)\leq u(x,s_1)$ in $[0,1]$;
\item[(c)] \eqref{max>0>min} holds at $s_1$;
\item[(d)] $\underline{u}(x,s_1;h) \geq u(x,s_1)$ in $ [0,1]$,
{\it equality} holds at some points and {\it strict inequality} holds at other points.
\end{itemize}
We now show that the cases (a), (b) and (c) are impossible, and Claim 1 holds in case (d).
In fact, in case (a), by continuity we have $\underline{u}(x,t;h)> u(x,t)$ for $t$ satisfying $s_1 <t\ll s_1 +1$.
In case (b), by comparison we have $\underline{u}(x,t;h) < u(x,t)$ for $x\in [0,1]$
and $t>s_1$. In case (c), there is a small time interval containing $s_1$ as its interior point
such that \eqref{max>0>min} holds in this interval.  Thus, all of these three cases contradict the
definition of $(s_1, s_2)$.  The only left possible case is (d). If $\underline{u}(1,s_1;h) > u(1,s_1)$,
then using the strong comparison principle in $[0,1]\times [s_1, t)$
for $s_1 < t \ll s_1 +1$ we conclude that $\underline{u}(x,t;h) > u(x,t)$ in this domain.
This again contradicts the definition of $(s_1, s_2)$. Therefore, $\underline{u}(1,s_1;h) = u(1,s_1)$,
and Claim 1 is proved.

\medskip
\noindent
{\bf Claim 2}. For any $t\in (s_1, s_2)$, $\underline{u}(\cdot,t;h)$ contacts $u(\cdot,t)$ at exactly one point
$x= X(t,h)$ in $(0,1)$ with $\underline{u}_x (X, t;h) <  u_x (X, t)$.

Since, by Claim 1,
\begin{equation}\label{eta-slope0}
\underline{u}_x (1,s_1;h) =h_0 < \underline{u} (1,s_1;h) = u(1,s_1) = u_x (1,s_1),
\end{equation}
we have $\underline{u} (x,s_1;h) > u(x,s_1)$ in $(x_* ,1)$ for some $x_*$
near $1$. Using Claim 1 and the comparison principle in the domain $\{(x,t)\mid 0<x<x_*,\ t>s_1 \mbox{ with } t-s_1 \ll 1\}$
we conclude that $\underline{u} (x,t;h) > u(x,t)$ in this domain.
Therefore, Claim 2 holds for $s_1 < t\ll s_1 +1$, and the unique contact point between
$\underline{u} (\cdot,t;h)$ and $u(\cdot,t)$ is near $1$.
If we can show that $\underline{u} (1,t;h) < u(1,t)$ for all $t\in (s_1, s_2)$,
then combining with the Neumann condition at $x=0$:
$$
\underline{u}_x (0,t;h) = u_x( 0,t) =0,
$$
we can derive Claim 2 by the zero number argument (\cite{Ang}).
Assume by contradiction that for some $s_3\in (s_1, s_2)$, $\underline{u} (1,t;h) < u(1,t)$ for $t\in (s_1, s_3)$,
but $\underline{u} (1,s_3;h) = u(1,s_3)$.  Then Claim 2 holds in $(s_1, s_3)$.  Due to
$\underline{u} (1,s_3;h) = u(1,s_3)$ and the analogue of \eqref{eta-slope0} at $s_3$,
there exists $x^* \in (0, 1)$ such that $\underline{u} (x,s_3;h) > u(x,s_3)$ for $x\in [x^*, 1)$.
By continuity, there exists a small $\epsilon>0$ such that $\underline{u} (x^*, t;h) > u(x^*,t)$
for $t\in  (s_3 -\epsilon, s_3 +\epsilon)$.  Then, in the time period $(s_3 -\epsilon, s_3)$, the unique contact point
$X(t,h)$ between $\underline{u}(\cdot,t;h)$ and $u(\cdot,t)$  lies in $(x^*, 1)$ since
$\underline{u} (x^*,t;h) > u(x^*,t)$ and $\underline{u} (1,t;h) < u(1,t)$ for $t\in (s_3 -\epsilon, s_3)$.
Therefore, $\underline{u} (x,t;h) > u(x,t)$ for $(x,t)\in [0,x^*] \times (s_3-\epsilon, s_3)$.
Using the strong maximum principle we see that $\underline{u} (x,s_3;h) > u(x,s_3)$ for $x\in [0,x^*]$.
Consequently, $\underline{u} (x,s_3;h) > u(x,s_3)$ in $[0,1)$.
This, however, contradicts the definition of $(s_1, s_2)$ and our assumption $s_3 \in (s_1, s_2)$.
Therefore, Claim 2 holds for all $t\in (s_1, s_2)$.

Using Claim 2 in particular at $t_1\in (s_1, s_2)$, we have $\underline{u}_x (X(t_1,h),t_1;h) < u_x(X(t_1,h),t_1)$, where
$X(t_1, h)$ is the unique contact point between $\underline{u} (\cdot,t_1;h)$ and $u(\cdot,t_1)$,
which is nothing but $x_1$. This proves \eqref{underlineu-x<u-x}. Since $t_1 \geq t^0+T^0$ and $x_1\in (0,1)$
are arbitrarily given, we actually obtain the following lemma.
\begin{lemma}\label{lem:small-gradient}
For any $t\geq t^0+T^0$ and any $x\in (0,1)$, there holds
$$
u_x( - x, t) \leq \varphi'(-x;h_0)<0,\qquad 0< \varphi'(x;h_0) \leq u_x(x, t).
$$
\end{lemma}

The inequalities in this lemma in fact can be improved to strict ones by using the strong maximum principle
and the zero number argument, which is omitted here, since the present version is enough for our approach below.

\subsection{Finer upper gradient estimate}
We now show that \lq\lq the gradient of $u$ is not too large" for $x\in (0,1)$.

By our assumption $u'_0(x)<\varphi'_0(x)$ for $x\in (0,1)$, we see that $u_0(x) \leq \varphi_0(x) + u_0(0)$ in $(-1,1)$, and so
\begin{equation}\label{4.6}
u(x,t)<\varphi_0(x) + u_0(0)  +\frac{\pi}{2}t,\qquad x\in(-1,1),\ t>0
\end{equation}
by the comparison principle. On the other hand, for any $r < u_0(0)$,
$\varphi_0(x) + r - u_0(x)$ has exactly two non-degenerate zeros.
Since $\zeta(x,t):=\varphi_0(x)+\frac{\pi}{2}t + r -u(x,t)$ satisfies a linear
parabolic equation whose coefficients are bounded in any compact interval of $(-1,1)\times (0,\infty)$,
and since
$$
\zeta(-1+0,t) = \zeta(1-0, t) =+ \infty,\quad t>0,
$$
we can use the zero number argument to conclude that, for any $t>0$,  either
\begin{itemize}
\item[(1)] $\zeta (\cdot, t)$ has two non-degenerate zeros $\pm\rho(t)$ with $\rho(t)\in(0,1)$; or
\item[(2)] $\zeta  (\cdot, t)$ has a unique degenerate zero $0$; or
\item[(3)] $\zeta  (\cdot, t)$ is positive, and has no zeros in $(-1,1)$.
\end{itemize}

Note that, for any $t>0$,  the graph of $u(\cdot,t)$ is immersed in
$$
\mathcal{E}(t) := \left\{\left. \left( x, \varphi_0(x) + r +\frac{\pi}{2}t\right) \right| x\in (-1,1),\ r\leq u_0(0)\right\}.
$$
Hence, for any $t_2>0$ and any $x_2 \in (0,1)$, there exists a unique $r = R(t_2,x_2)<u_0(0)$ such that
$\zeta (\cdot, t_2)$ with $r=R(t_2,x_2)$ has exactly two zeros $x= \pm x_2$. Consequently,
\begin{equation}\label{not too large}
\varphi'_0(-x_2) < u_x(- x_2,t_2) < 0,\qquad 0< u_x(x_2,t_2) < \varphi'_0(x_2).
\end{equation}
In this sense, we say that \lq\lq the gradient of $u$ is not too large" for $x\in (0,1)$.

\subsection{Convergence of the solution}

Assume $u(x,t)$ is a symmetric solution starting from an initial data satisfying \eqref{ini-1}.
Let $\{t_n\}$ be a time sequence with $t_n\rightarrow \infty\ (n\rightarrow\infty)$. Set
$$
u_n(x,t):=u(x,t+t_n)-u(0,t_n),\qquad x\in[-1,1],\ -t_n<t<\infty.
$$
For any given small $\varepsilon>0$, by the results in the previous two subsections we have
\begin{equation}\label{5.1}
\varphi' (x;h_0)<u_{nx}(x,t)<\varphi'_0(x),\qquad x\in (0,1-\varepsilon],\ n\gg 1.
\end{equation}
$h_0>0$ can be as large as possible, provided $n$ is sufficiently large and $t$ is bounded. Combining with Lemma \ref{bound-est} we have
$$
\left\|u_n(x,t) \right\|_{C^{1,0}\left([\varepsilon-1,1-\varepsilon]\times [-T,T]\right)}\leq C_1(\varepsilon,T),
$$
for any $T>0$. By the $L^p$ estimates, Sobolev embedding theorem and the Schauder estimate we have,
for any $\alpha\in(0,1)$,
$$
\left\|u_n(x,t) \right\|_{C^{2+\alpha,1+\frac{\alpha}{2}}\left([\varepsilon-1,1-\varepsilon]\times [-T,T]\right)}\leq C_2(\varepsilon,T).
$$
Therefore, for any $\beta\in (0,\alpha)$, there exist  a subsequence $\{u_{n_i}\}$ of $\{u_n\}$ and a function
$\mathcal{U}_{T,\varepsilon}\in C^{2+\alpha,1+\frac{\alpha}{2}} $ $\left([\varepsilon-1,1-\varepsilon]\times [-T,T]\right)$
such that
$$
\left\|u_{n_i}-\mathcal{U}_{T,\varepsilon}\right\|_{C^{2+\beta,1+\frac{\beta}{2}}
\left([\varepsilon-1,1-\varepsilon]\times [-T,T]\right)}\rightarrow 0\quad  (i\rightarrow \infty).
$$
Using the Cantor's diagonal argument, we see that there exist a function $\mathcal{U}\in C^{2+\alpha,1+\frac{\alpha}{2}}_{loc}
\left((-1,1)\times \R\right)$ and a subsequence of $\{u_n\}$ (denoted it again by $\{u_{n_i}\}$) such that
 $$
u_{n_i}\rightarrow \mathcal{U}\ \ (i\to \infty), \quad\text{in}\quad C^{2+\beta,1+\frac{\beta}{2}}_{loc}\left((-1,1)\times \R\right) \quad \text{topology}.
 $$
Moreover, $\mathcal{U}(x,t)$ is an entire solution of the equation in (\ref{B}) with $\mathcal{U}(0,0)=0$ and, by (\ref{5.1}),
$$
\varphi' (x;h_0)\leq \mathcal{U}_x(x,t)\leq \varphi'_0(x),\qquad  x\in [0,1),\ t\in \R.
$$
Since $h_0>0$ can be as large as possible and since
$$
\varphi' (x;h_0)\rightarrow \varphi'_0(x) \mbox{ for } x\in (-1,1),\qquad\quad \text{as}\quad h_0\rightarrow \infty,
$$
we conclude that
$$
\mathcal{U}_x(x,t) =  \varphi'_0(x),\qquad x\in (-1,1),\ t\in \R,
$$
and so
$$
\mathcal{U}(x,t)= \varphi_0(x)+C(t),\qquad x\in (-1,1),\ t\in \R.
$$
Finally, since $\mathcal{U}$ is an entire solution of the equation in (\ref{B}), we see that
$C'(t)=\frac{\pi}{2} $ and thus
$$
\mathcal{U}(x,t)=\varphi_0(x)+\frac{\pi}{2}t,\qquad x\in(-1,1),\ t\in \R.
$$

The above result implies that $\{u_{n_i}\}$ converges to the special {\it grim reaper} $\varphi_0(x)+\frac{\pi}{2}t$.
Since this {\it grim reaper} is unique and the time sequence $\{t_n\}$ is arbitrarily given, we actually obtain the following result.

\begin{thm}\label{thm:converge-symmetric}
Assume $u(x,t)$ is the time-global solution of \eqref{B} with initial data $u_0(x)$ satisfying \eqref{ini-1}. Then,
for any $\alpha\in (0,1)$,
\begin{equation}\label{5.2}
u(x,t+s)-u(0,t)\rightarrow \varphi_0(x)+\frac{\pi}{2}t, \qquad \text{as} \ s\rightarrow \infty,
\end{equation}
in the $C^{2+\alpha,1+\frac{\alpha}{2}}_{loc}\left((-1,1)\times \R \right)$ topology.
\end{thm}

\section{General Solutions}
The conclusion in the previous section holds only for symmetric solutions.
We consider general solutions in this section, that is, we assume $u_0$
satisfies \eqref{1.4} in this section.

\subsection{Interior estimates}
We choose a smooth, even function $\psi (x)$ with
$$
0< \psi'(x)<\varphi'_0(x),\quad 0<\psi(x) <u_0(x),\quad \psi''(x)>0 \quad \mbox{for } x\in (0,1],\quad \psi'(1)=\psi(1).
$$
Then $\psi$ satisfies \eqref{ini-1} except for $\psi(x)\geq 1$. One example is $\psi(x)= \delta \Big[ \frac{\sqrt{2}}{2}
- \frac{4}{\pi +4} \cos \frac{\pi x}{4} \Big]$ for small $\delta>0$. Denote the solution of (\ref{B})
with $u(x,0)=\psi(x)$ by $u(x,t;\psi)$. Then it is symmetric, $u_{xx}(x,t;\psi) >0$ due to $u_t(x,t;\psi)>0$,
and it satisfies all the conclusions in the previous section. Furthermore, it moves upward monotonically to infinity,
so we have
$$
u_0(x) < u(x,T;\psi),\qquad x\in[-1,1],
$$
for some positive $T$. Thus, by the comparison principle we have
\begin{equation}\label{4.9}
u(x,t;\psi) < u(x,t; u_0)< u(x,t+T;\psi) ,\qquad x\in [-1,1],\ t>0.
\end{equation}
This formula gives the $L^\infty$ estimate for $u(x,t; u_0)$.

In what follows, we want to present a uniform interior gradient estimate.
First we prove a lemma.

\begin{lemma}\label{lem4.3}
For any small $\varepsilon\in(0, \frac12)$ and any $t>0$, there hold
\begin{equation}\label{4.10}
\min_{1-2\varepsilon\leq x\leq 1-\varepsilon}\left|u_x(x,t;u_0)\right|<M_1 :=\varepsilon^{-1} \Big[\varphi_0(1-\varepsilon)+\frac{\pi}{2}T \Big],\qquad
\min_{\varepsilon-1 \leq x\leq 2\varepsilon-1}\left|u_x(x,t;u_0)\right|<M_1.
\end{equation}
\end{lemma}

\begin{proof}
We only prove the first inequality since the second one is proved similarly.
Assume by contradiction that, for some $t=t_0>0$,
$$
\left|u_x(x,t_0;u_0)\right| \geq M_1, \qquad x\in[1-2\varepsilon, 1-\varepsilon].
$$
Integrating it over $[1-2\varepsilon, 1-\varepsilon]$ we obtain
\begin{equation}\label{4.11}
u\left(1-\varepsilon, t_0; u_0\right)-u\left(1-2\varepsilon, t_0; u_0\right)
 \geq \varphi_0\left(1-\varepsilon\right)+\frac{\pi}{2}T.
\end{equation}

On the other hand, by (\ref{4.9}) we have
\begin{equation}\label{4.12}
 \begin{array}{lll}
 u \left(1-\varepsilon, t_0; u_0\right)-u\left(1-2\varepsilon, t_0; u_0\right)
 & < & u\left(1-\varepsilon, t_0+T; \psi\right)-u\left(1-2\varepsilon, t_0; \psi\right) \\
 & \leq & u\left(1-\varepsilon, t_0+T; \psi\right)-u\left(0, t_0; \psi\right)
 \end{array}
\end{equation}
since $u(\cdot, t; \psi)$ is convex and symmetric.
By \eqref{not too large}, $u(x,t_0;\psi)<\varphi_0(x)+u(0,t_0;\psi)$, and so by comparison we have
$$
u(x,t_0 +T;\psi)<\varphi_0(x)+\frac{\pi}{2}T+u(0,t_0;\psi).
$$
In particular, at $x=1-\varepsilon$ we have
$$
u\left(1-\varepsilon, t_0+T; \psi\right)-u\left(0, t_0; \psi\right)< \varphi_0(1-\varepsilon)+\frac{\pi}{2}T.
$$
This contradicts (\ref{4.11}) and (\ref{4.12}). This proves the lemma.
\end{proof}

Using this lemma we can prove the following interior gradient estimates.

\begin{lemma}\label{lem:interior-g-est}
For any small $\varepsilon>0$, there exists $T_\varepsilon >0$ such that
\begin{equation}\label{interior-g-est}
 |u_x (x,t; u_0)| \leq M_2, \quad   -1+2\varepsilon <x <1-2\varepsilon,\ t>T_\varepsilon,
 \end{equation}
where $M_2 := \max\{M_1, \|u_x(\cdot, T_\varepsilon; u_0)\|_{L^\infty}\}$ and $M_1$ is that in \eqref{4.10}.
\end{lemma}

\begin{proof}
Since $u(x,t;u_0) \to \infty$ as $t\to \infty$, there exists $ T' >0$ large such that $u(\pm 1, t; u_0)
> M_1$ for all $t>T' $. Denote $\zeta(x,t):= u_x(x,t; u_0)-M_1$, then $\zeta$ solves
$$
\zeta_t = \frac{\zeta_{xx}}{1+u^2_x} - \frac{2 u_x }{(1+u^2_x)} \zeta^2_x,\quad -1<x<1,\ t>0,
$$
and $\zeta (1,t)>0>\zeta(-1,t)$ for $t>T'$. Using the zero number properties (cf. \cite{Ang})
we conclude that, for some $T_\varepsilon  >T' $, the function $\zeta(\cdot, t)$ has only non-degenerate zeros for
$t\geq T_\varepsilon $. Denote the largest zero of $\zeta(\cdot,t)$ in $(-1,1)$ by $ \rho_+(t)$. Due to
the non-degeneracy of $\rho_+(t)$ we see that $x= \rho_+(t)$ is a continuous curve. Moreover,
\eqref{4.10} in the previous lemma indicates that $\rho_+(t) > 1- 2\varepsilon$.
In a similar way we can find another continuous curve $x=\rho_-(t)$ for $t>T_\varepsilon $ ($T_\varepsilon $ can be chosen larger if necessary),
with $\rho_-(t) \in (-1, -1 + 2\varepsilon)$, such that $u_x(\rho_-(t),t) = - M_1$ for $t >T_\varepsilon$.
Then, using the maximum principle for $u_x$ in the domain $D(T_\varepsilon ) := \{(x,t) \mid \rho_-(t)<x<\rho_+(t),\ t>T_\varepsilon \}$
we conclude that $|u_x (x,t;u_0)|\leq M_2$ in $D(T_\varepsilon)$. The estimate \eqref{interior-g-est} then follows
from the fact $\rho_-(t)<-1+2\varepsilon < 1-2\varepsilon < \rho_+(t)$ for $t>T_\varepsilon$.
\end{proof}

\subsection{Convergence of general solutions}

Let $\{t_n\}$ be any time sequence with $t_n\to \infty$, we consider the solution sequence $\{ u(x,t+t_n; u_0)-u(0,t_n;\psi)\}$.

For any given small $\varepsilon>0$ and any $\tau >0$, let $T_\varepsilon$ be as that in Lemma \ref{lem:interior-g-est},
then \eqref{4.9} and \eqref{interior-g-est} imply that, for all large $n$, $u(x,t+t_n; u_0)-u(0,t_n;\psi)$ is bounded
in $C^{1,0}([2\varepsilon-1, 1-2\varepsilon]\times [-\tau, \tau])$ norm, and the bounds are independent of $n$.
Using the standard parabolic theory we can even have the $C^{2+\alpha,1+\frac{\alpha}{2}}
([2\varepsilon-1, 1-2\varepsilon] \times [-\tau, \tau])$ (for any $\alpha\in (0,1)$) bounds for the solution sequence, also uniform in $n$, and
so we can find a convergent subsequence.
Taking $\varepsilon \to 0$ and $\tau \to \infty$, and using the Cantor's diagonal argument we conclude that,
there is a subsequence of $\{t_n\}$, denoted is again by $\{t_n\}$, such that
\begin{equation}\label{u-n-to-W}
u(x,t+t_n; u_0)-u(0,t_n;\psi)\rightarrow \mathcal{W}(x,t) \qquad \text{in}\quad
C^{2+\alpha,1+\frac{\alpha}{2}}_{loc}\left((-1,1)\times \R \right)\quad \text{topology},
\end{equation}
for some entire solution $\mathcal{W}$ of the equation in (\ref{B}).
On the other hand, as a consequence of Theorem \ref{thm:converge-symmetric} we have, as $n\to \infty$,
$$
u(x,t+t_n;\psi)-u(0, t_n;\psi)\rightarrow \varphi_0(x)+\frac{\pi}{2}t, \qquad u(x,t+T+t_n;\psi)
-u(0,t_n;\psi)\rightarrow \varphi_0(x)+\frac{\pi}{2}(t+T).
$$
Hence, we conclude from (\ref{4.9}) that
\begin{equation}\label{5.6}
\varphi_0(x)+\frac{\pi}{2}t\leq\mathcal{W}(x,t)\leq \varphi_0(x)+\frac{\pi}{2}t+\frac{\pi}{2}T,\qquad x\in(-1,1),\ t\in \R.
\end{equation}

Denote
$$
\theta(x,t):= \arctan u_x(x,t),\quad x\in [-1,1],\ t>0.
$$
Then $\theta$ satisfies
$$
\left\{
 \begin{array}{ll}
 \theta_t = \cos^2 \theta \cdot \theta_{xx}, & -1<x<1,\ t>0,\\
 \theta(\pm 1,t)=  \pm \arctan u(\pm 1,t), & t>0.
 \end{array}
 \right.
$$
The global existence of $u$ implies that $\theta$ exists for all $t>0$. For any $T>0$, by Lemma \ref{bound-est} we have
\begin{equation}\label{bound-theta}
|\theta(x,t)|\leq \arctan [C(T)],\quad x\in [-1,1],\ t\in [0,T],
\end{equation}
and $\theta(\pm 1,t) = \pm \arctan u(\pm 1,t) \to \pm \frac{\pi}{2}$ as $t\to \infty$.
Moreover, by \eqref{u-n-to-W} we have
$$
\theta(x,t_n +t) \to \Theta(x,t) := \arctan [\mathcal{W}_x (x,t)] \mbox{ in } C^{2+\alpha, 1+\alpha/2}_{loc} ((-1,1)\times \R)
\mbox{ topology}.
$$

1. We claim that $\Theta(x,t)$ is a stationary solution of $\theta_t = \cos^2 \theta \cdot \theta_{xx}$. Without loss of generality
we only need to prove that
$\Theta(x,0)$ is a stationary one. Denote $a:= \Theta_x (0,0),\ b:= \Theta(0,0)$ and $v(x):= ax+b$,
then $v(x)$ is a stationary solution of $\theta_t = \cos^2 \theta \cdot \theta_{xx}$, and the function
$\xi (x,t) := \theta (x,t) -v(x)$ satisfies one of the following boundary conditions:
$$
\begin{array}{ll}
 \mbox{in case }  v(-1) >-\frac{\pi}{2},\ v(1) \geq \frac{\pi}{2}: & \xi(-1,t) < 0,\ \xi (1,t)<0 \mbox{ for all large } t;\\
  \mbox{in case }  v(-1) >-\frac{\pi}{2},\ v(1) < \frac{\pi}{2}:  & \xi(-1,t) < 0,\ \xi(1,t)>0 \mbox{ for all large } t;\\
  \mbox{in case }   v(-1) \leq -\frac{\pi}{2},\ v(1) \geq \frac{\pi}{2}: & \xi(-1,t) > 0,\ \xi(1,t)<0 \mbox{ for all large } t;\\
  \mbox{in case }   v(-1) \leq -\frac{\pi}{2},\ v(1) < \frac{\pi}{2}: & \xi(-1,t) > 0,\ \xi(1,t)>0 \mbox{ for all large } t.
\end{array}
$$
Moreover, $\xi$ satisfies the linear equation $\xi_t = \cos^2 \theta\cdot \xi_{xx}$.
Hence the zero number argument is applied to conclude that $\xi(x,t)$ has no degenerate zeros for all large $t$.
Then using a similar argument as in the proof of \cite[Lemma 2.6]{DM} we conclude that the $\omega$-limit $\Theta(x,t)-v(x)$
of $\theta (x,t)-v(x)$ either satisfies (1) $\Theta(x,t)\equiv v(x)$, or satisfies (2) $\Theta(x,t)\not \equiv v(x)$ and
$\Theta(x,t)-v(x) $ has no degenerate zeros for each $t\in \R$. Case (2), however,
contradicts the definition of $v(x)$. Therefore, only case (1) holds. This prove our claim.

2. Next we show that $\Theta (x,t) \equiv \frac{\pi}{2}x$. In the previous step, we have shown that
$\Theta(x,t)\equiv v(x)\equiv ax +b$ for some  $a,b\in \R$.
If $v(\pm 1)= \pm \frac{\pi}{2}$, then the conclusion is proved.
On the contrary, we assume, without loss of generality that, for some $x_0\in (0,1)$ and some small $\delta >0$, there holds
$\Theta(x,t)\equiv v(x) <\frac{\pi}{2}x- \pi \delta$ in $x\in [x_0,1)$.
For any given small $\epsilon>0$, since
$|\theta(x,t_n) - \Theta(x,0)|\leq \frac12 \pi \delta$ in $[-1+\epsilon, 1-\epsilon]$
when $n$ is sufficiently large, we have
$$
 \arctan u_x (x,t_n) = \theta (x,t_n) \leq \Theta(x,0) + \frac12 \pi \delta = v(x) + \frac12 \pi \delta
 < \frac{\pi}{2} (x - \delta),\quad x\in [-1+\epsilon, 1-\epsilon].
$$
Integrating the above inequalities over $[x_0,1-\epsilon]$ we have
$$
u(1-\epsilon,t_n) \leq u(x_0,t_n) + \int_{x_0}^{1-\epsilon} \tan \Big[\frac{\pi}{2} (x -\delta) \Big] dx= u(x_0, t_n) +
\varphi_0 (1-\epsilon-\delta)-\varphi_0(x_0-\delta).
$$
Taking limit as $n\to \infty$ we have
$$
\mathcal{W}(1-\epsilon,0)\leq \mathcal{W}(x_0,0) + \varphi_0 (1-\epsilon-\delta)-\varphi_0(x_0-\delta).
$$
Combining with \eqref{5.6} we have
$$
\varphi_0(1-\epsilon) \leq \mathcal{W}(x_0,0) + \varphi_0 (1-\epsilon-\delta)-\varphi_0(x_0-\delta).
$$
This is a contradiction as $\epsilon\to 0$. This proves our claim, and so $\Theta(x,t)\equiv \frac{\pi x}{2}$.

3. As an $\omega$-limit of $\theta$, the function $\Theta(x,t)\equiv \frac{\pi}{2}x$ is unique. Hence, $\theta(x,t)\to \frac{\pi}{2}x$ as
$t\to \infty$. Equivalently, $\mathcal{W}_x (x,t)\equiv \tan [\frac{\pi}{2}x]$, and so
$$
\mathcal{W}(x,t) = \varphi_0(x) + K(t).
$$
Since $\mathcal{W}(x,t)$ is an entire solution, by its equation we have
$\mathcal{W}_t(x,t)\equiv \frac{\pi}{2}$, and so $K(t) = \frac{\pi}{2}t +K_0$
for some $K_0 \in \R$. By \eqref{5.6} we actually have $ K_0 \in [0, \frac{\pi}{2}T]$.
In summary,
$$
\mathcal{W}(x,t) = \varphi_0(x)+\frac{\pi}{2}t+ K_0, \qquad x\in(-1,1),\ t\in \R.
$$

\medskip
\noindent
{\it Proof of Theorem} \ref{thm:main}. From above discussion we see that, as $n\to \infty$,
$$
u(x,t+t_n; u_0) - u(0,t_n;\psi) \to \varphi_0(x) +\frac{\pi}{2} t + K_0,
$$
and so
\begin{eqnarray*}
u(x,t+t_n;u_0) - u(0,t_n;u_0) & = & [u(x,t+t_n;u_0) - u(0,t_n;\psi)] - [u(0,t_n;u_0) - u(0,t_n;\psi)]\\
 & \to & \varphi_0(x) + \frac{\pi}{2} t .
\end{eqnarray*}
Since $\{t_n\}$ is an arbitrarily chosen time sequence tending to infinity, we conclude that
$$
u(x,t+s;u_0) - u(0,s;u_0) \to \varphi_0 (x) +\frac{\pi}{2} t,\quad \mbox{ as } s\to \infty.
$$
This proves Theorem \ref{thm:main}. \hfill $\Box$

\begin{remark}\rm
We see from Sections 4 and 5 that the zero number argument plays a key role in the {\it uniform}
interior gradient estimates, which does not rely on the boundedness of $u$ and can be
extended to other problems with unbounded solutions (but only for one dimensional problems).
\end{remark}