\section{Introduction}

	Sequence segmentation, as an important task in natural language understanding (NLU), partitions a sentence into multiple segments. The first two rows of Table \ref{tab:Example Demo} show a case from a syntactic chunking dataset. The input sentence is a sequence of tokens and the output segments are multiple labeled phrases. These segments are nonoverlapping and fully cover the input sentence.

	In previous works, there are two dominant approaches to sequence segmentation. The most common is to regard it as a sequence labeling problem with resorting to IOB tagging scheme~\citep{huang2015bidirectional,akbik2018contextual,liu-etal-2019-gcdt}. This method is simple yet very effective, providing tons of state-of-the-art performances. For example, \newcite{huang2015bidirectional} present Bidirectional LSTM-CRF for named entity recognition (NER) and POS tagging, which adopts BiLSTM~\citep{hochreiter1997long} to read the input sentence and CRF~\citep{lafferty2001conditional} to decode the label sequence. An alternative method employs a transition-based system to incrementally segment and label an input token sequence~\citep{zhang2016transition,zhang2018simple}. For instance, \newcite{zhang2016transition} present a transition-based model for Chinese word segmentation that exploits not only character embedding but also token embedding. This type of method enjoys a number of attractive properties, including theoretically lower time complexity and capturing non-local features.

	\begin{table}[t]
		\centering

		\setlength{\tabcolsep}{1.1mm}{}
		\begin{tabular}{c|c}

			\hline
			\multirow{2}{*}{Sentence} & Tangible capital will be  about \\
			& \$ 115 million . \\

			\hline
			\multirow{2}{*}{Segments} & (Tangible capital, NP), (will be, VP), \\
			& (about \$ 115 million, NP), (., O)  \\
			\hline

			\multirow{2}{*}{IOB Tags} & B-NP I-NP B-VP I-VP \\

			& B-NP I-NP I-NP I-NP O  \\

			\hline
			& SHIFT SHIFT REDUCE-NP \\

			Transition &SHIFT SHIFT REDUCE-VP \\

			Actions	& SHIFT SHIFT SHIFT SHIFT \\

			& REDUCE-NP OUT \\

			\hline
		\end{tabular}

		\caption{The first two rows show an example extracted from CoNLL-2000 dataset~\citep{sang2000introduction}. The last two rows are two types of token-level labels commonly used to represent the segments.}
		\label{tab:Example Demo}
	\end{table}

	\begin{figure*}
		\centering
		\includegraphics[width=0.85\textwidth]{images/1.pdf}

		\caption{This case illustrates how our method segments a natural language sentence.}
		\label{fig:Approach}
	\end{figure*}

	The above two approaches are essentially at token level, where a single segment is represented by multiple token-level labels (e.g., transition actions). In spite of the adequacy, the labels used to model the relation among output segments are far more than the segments themselves. As demonstrated in Figure \ref{tab:Example Demo}, modeling the transition between the two segments, ``will be" and ``about \$ 115 million", consumes 6 IOB tags or 8 transition actions. This ill-posed design certainly limits the full potential of segmentation models to capture the long-term dependencies among segments.

	Previously, \citet{kong2015segmental} attempted to develop segment-level models, which define a joint probability distribution over the partition of an input sentence and the labeling of the segments. Such an approach circumvents using token-level labels. However, we find that, in experiments, it underperforms current token-level models. Moreover, its use of dynamic programming (DP) incurs quadratic running time, which is too slow for both training and inference (see Section \ref{sec:Running Time Analysis}).

	In this paper, we introduce a novel framework that incrementally segments a sentence at segment level. The segmentation process is iterative and incremental. At each iteration, the proposed framework determines the leftmost segment of the remaining unprocessed sequence. Under this scheme, we don't resort to token-level labels and enjoy linear time complexity. The implementation contains two stages. Firstly, we utilize LSTM-minus~\citep{wang2016graph,cross-huang-2016-span} technique to construct the representations for all the phrases. Secondly, we adopt LSTM to model the iterative segmentation process, which captures the strong correlation among segments. At every step, the input consists of the previous segment and the remaining unprocessed sequence, and the output is the leftmost segment.

	Figure \ref{fig:Approach} depicts how our framework segments the sentence in Table \ref{tab:Example Demo}. The output segments are obtained in an iterative and incremental manner. At each iteration, the leftmost segment of the remaining sequence is extracted and labeled. Compared with token-level models, we take much fewer steps to complete the segmentation process.

	Extensive experiments have been conducted on syntactic chunking and Chinese part-of-speech (POS) tagging across 3 datasets. The proposed framework has obtained new state-of-the-art performances on all of them. Besides, qualitative study and the results on segmenting long-length sentences confirm its effectiveness in capturing long-term dependencies.

	Our contributions are as follows:
	\begin{itemize}

		\item we present a novel framework that incrementally segments a natural language sentence at segment level. In comparison, previous approaches are mostly at token-level;

		\item we have notably outperformed previous baselines and established new state-of-the-art results on the 3 datasets of syntactic chunking and Chinese POS tagging. Experiments also show that our model is competitive with strong NER baselines;

		\item compared with prior methods, our model well captures the long-term dependencies among segments. This is strongly verified by qualitative study and the experiment on segmenting long-length sentences.

	\end{itemize}

	The source code of this work is available at https://github.com/LeePleased/LeftmostSeg.

\section{Architecture}

	We denote a $n$-length input sentence as $\mathbf{x} = [x_1, x_2, \cdots, x_n]$, where $x_i$ is a token (such as a word or a character). The output segments are represented as $\mathbf{y} = [y_1, y_2, \cdots, y_m]$, where $m$ is the amount of segments. Every segment $y_k$ is denoted as a triple $(i_k, j_k, l_k)$. $(i_k, j_k)$ is the span of the segment, corresponding to the phrase $\mathbf{x}_{i_k, j_k} = [\mathbf{x}_{i_k}, \mathbf{x}_{i_k + 1}, \cdots, \mathbf{x}_{j_k}]$. $l_k$ is from a predefined label space $\mathcal{L}$ and specifies the label of the segment. These segments are non-overlapping and fully cover the input sentence.

	\begin{figure*}
		\centering
		\includegraphics[width=0.8\textwidth]{images/2.pdf}

		\caption{An example to show the incremental process to recognize leftmost segments.}
		\label{fig:Transduction}
	\end{figure*}

	The example in Table \ref{tab:Example Demo} is represented as
	\begin{equation}\nonumber
	\begin{aligned}
	\mathbf{x} =  [& \mathrm{Tangible}, \mathrm{capital}, \mathrm{will}, \mathrm{be}, \mathrm{about}, \mathrm{\$}, \mathrm{115}, \\
	& \mathrm{million}, \mathrm{.}] \\
	\mathbf{y} = [& (1, 2, \mathrm{NP}), (3, 4, \mathrm{VP}), (5, 8, \mathrm{NP}), (9, 9, \mathrm{O})]
	\end{aligned}.
	\end{equation}

\subsection{Phrase Representation Construction}

	Our goal here is to construct the representation for every phrase $\mathbf{x}_{i,j}$. Later, we will use the phrase representation to embed the segment $y_{k-1}$ and the unprocessed sequence $\mathbf{x}_{i_k, n}$.

	Firstly, each token $x_k$ is represented as
	\begin{equation}
	\mathbf{e}_k = \mathbf{E}^{\text{t}}(x_k) \oplus \mathrm{CharCNN}(x_k),
	\end{equation}
	where $\mathbf{E}^\text{t}$ is a token embedding matrix and $\oplus$ is the column-wise vector concatenation. Following previous works~\citep{ma-hovy-2016-end,liu-etal-2019-gcdt}, we use $\mathrm{CharCNN}$ to extract the character-level representations.

	Secondly, we utilize bidirectional LSTMs $\overrightarrow{f}^{\text{e}}$ and $\overleftarrow{f}^{\text{e}}$ to compute the context-sensitive representation for each token $x_k$:
	\begin{equation}
	\left\{\begin{aligned}
	\overrightarrow{\mathbf{h}}^{\text{c}}_k & = \overrightarrow{f}^{\text{e}}(\overrightarrow{\mathbf{h}}^{\text{c}}_{k-1}, \mathbf{e}_k) \\ \overleftarrow{\mathbf{h}}^{\text{c}}_k & = \overleftarrow{f}^{\text{e}}(\overleftarrow{\mathbf{h}}^{\text{c}}_{k+1}, \mathbf{e}_k)  \\
	\mathbf{h}^{\text{c}}_k & = \overrightarrow{\mathbf{h}}^{\text{c}}_k \oplus \overleftarrow{\mathbf{h}}^{\text{c}}_k
	\end{aligned}\right..
	\end{equation}

	Inspired by previous works~\citep{wang2016graph,gaddy-etal-2018-whats} in syntactic analysis, we integrate LSTM-minus features into phrase representations. Specifically, the representation $\mathbf{h}^\text{p}_{i, j}$ for a phrase $\mathbf{x}_{i,j}$ is computed as the concatenation of the difference of LSTM hidden states:
	\begin{equation}
	\mathbf{h}^\text{p}_{i, j} = \mathbf{h}^{\text{c}}_j \oplus (\mathbf{h}^{\text{c}}_j - \mathbf{h}^{\text{c}}_i) \oplus \mathbf{h}^{\text{c}}_i.
	\end{equation}

	Inside algorithm~\citep{lari1990estimation} is another method to extracting phrase representation. Its advantage is to incorporate the potential hierarchical structure of natural language without using treebank annotation. For example, \citet{drozdov-etal-2019-unsupervised-latent} utilize inside algorithm to recursively compute the content representations. Despite the attractive property, its time complexity $\mathcal{O}(n^3)$  is too inefficient to practice use.

\subsection{Leftmost Segment Determination}

	Figure \ref{fig:Transduction} demonstrates how our method iteratively and incrementally segments the sentence in Table \ref{tab:Example Demo}. LSTM is used as the backbone to model the iterative process. At every step, the input consists of the prior segment and the unprocessed sequence, and the output is the predicted segment.

	Firstly, we embed the previous segment $y_{k-1}$ and the unprocessed sequence $\mathbf{x}_{i_{k},n}$. The previous segment is represented as
	\begin{equation}
		\mathbf{h}^\text{s}_{k-1} = \left\{\begin{array}{cc}
		\mathbf{h}^\text{p}_{i_{k - 1}, j_{k-1}} \oplus \mathbf{E}^\text{l}(l_{k-1}) & k > 1 \\
		\mathbf{v} & k = 1
		\end{array}\right.,
		\label{equ:Equation 7}
	\end{equation}
	where $\mathbf{E}^\text{l}$ is a label embedding matrix and $\mathbf{v}$ is a trainable vector. The unprocessed sequence is embedded as $\mathbf{h}^\text{p}_{i_k,n}$.

	At each iteration $k$, we use another LSTM $f^\text{d}$ to model the dependency among segments:
	\begin{equation}
	\mathbf{h}^\text{d}_k = f^\text{d}(\mathbf{h}^\text{d}_{k-1}, \mathbf{h}^\text{s}_{k-1} \oplus \mathbf{h}^\text{p}_{i_k,n}).
	\label{equ:Equation 8}
	\end{equation}
	During training, $i_k$ is known since the ground truth segments $\mathbf{y}$ is reachable. At evaluation time, we set $i_k = j_{k-1} + 1$.

	\begin{algorithm*}
		\caption{Inference Procedure}
		\label{algo:inference}

		\KwIn{The representations for all the phrases, $\mathbf{h}^\text{p}_{i,j}, 1 \le i \le j \le n$.}
		\KwOut{The sequence of predicted segments, $\hat{\mathbf{y}} = [\hat{y}_1, \hat{y}_2, \cdots, \hat{y}_{\hat{m}}]$.}

		Set a list $\mathbf{y}$ as $[]$ and set a counter as $k = 1$. \\
		Denote the remaining input token list as $\mathbf{x}$. \\
		Initialize the representation for previous segment as $\mathbf{v}$. \\
		Initialize the representation for unprocessed sequence as $\mathbf{h}^\text{p}_{1,n}$. \\

		\While{$\mathbf{x}$ is not empty}{
			Get LSTM hidden state $\mathbf{h}^\text{d}_k$ by using Equation \ref{equ:Equation 8}. \\
			Predict a segment $\hat{y}_k$ by using the Equations from \ref{equ:Equation 10} to \ref{equ:Equation 12}. \\
			Append the new segment $\hat{y}_k$ into list $\hat{\mathbf{y}}$. \\
			Reset the representation for previous segment as $\mathbf{h}^\text{p}_{\hat{i}_{k}, \hat{j}_{k}} \oplus \mathbf{E}^\text{l}(\hat{l}_{k})$. \\
			Reset the representation for unprocessed sequence as $\mathbf{h}^\text{p}_{\hat{j}_{k} + 1,n}$. \\
			Pop the tokens $[x_{\hat{i}_k}, x_{\hat{i}_k + 1}, \cdots, x_{\hat{j}_k}]$ from the remaining tokens $\mathbf{x}$. \\
			Increase the counter by 1: $k = k + 1$.
		}

		The amount of predicted segments: $\hat{m} = k$.

	\end{algorithm*}

	Then, we separately predict the span and the label of a segment. We define a set $\mathcal{S}_k$ containing all valid span candidates for prediction:
	\begin{equation}
		\mathcal{S}_k = \{(i_k, i_k), (i_k, i_k + 1), \cdots, (i_k, n)\}.
	\end{equation}
	The probability of a span $(i,j) \in \mathcal{S}_k$ is
	\begin{equation}
	Q^\text{s}_{k,i,j} \propto \exp \big((\mathbf{h}^\text{d}_k)^{\top}\mathbf{W}^\text{s}\mathbf{h}^\text{p}_{i, j} \big).
	\label{equ:Equation 10}
	\end{equation}
	The probability of a label $l \in \mathcal{L}$ for a span $(i, j)$ is
	\begin{equation}
		Q^\text{l}_{k,i,j,l} \propto \exp\big(\mathbf{E}^\text{l}(l)^{\top} \mathbf{W}^\text{l} (\mathbf{h}^\text{p}_{i,j} \oplus \mathbf{h}^\text{d}_k ) \big).
		\label{equ:Equation 11}
	\end{equation}
	The matrices $\mathbf{W}^\text{s}$ and $\mathbf{W}^\text{l}$ are learnable.

	Finally, the leftmost segment is obtained as
	\begin{equation}
	\left\{\begin{aligned}
	(\hat{i}_k, \hat{j}_k) & = \mathop{\arg\max}_{(i,j) \in \mathcal{S}_k} Q^\text{s}_{k,i,j} \\
	\hat{l}_k & = \mathop{\arg\max}_{l \in \mathcal{L}} Q^\text{l}_{k, i, j, l} \\
	\hat{y}_k & = (\hat{i}_k, \hat{j}_k, \hat{l}_k)
	\end{aligned}\right..
	\label{equ:Equation 12}
	\end{equation}

	The iterative process ends when the remaining sequence $\mathbf{x}_{j_k + 1, n}$ is empty (i.e., $j_k = n$).

\subsection{Training and Inference}

	During training, we use teacher forcing where every segment $y_k$ is predicted using its previous ground-truth segments $[y_1, y_2, \cdots, y_{k-1}]$. A hybrid loss is induced as
	\begin{equation}
	\mathcal{J} = - \sum_{y_k \in \mathbf{y}} (\log Q^\text{s}_{k,i_k,j_k} + \log Q^\text{l}_{k,i_k,j_k,l_k} ).
	\label{equ:Equation 13}
	\end{equation}

	At test time, every segment $\hat{y}_k$ is inferred in terms of the previous predicted segments $[\hat{y}_{1}, \hat{y}_2, \cdots, \hat{y}_{k-1}]$. Algorithm \ref{algo:inference} demonstrates how our proposed framework makes inference. Note that Algorithm \ref{algo:inference} uses greedy search to get $\hat{y}$ because it is both fast in speed and effective in accuracy in our experiments, although beam search may be better in accuracy.

\section{Experiments}

	\begin{table*}
		\centering

		\setlength{\tabcolsep}{7.1mm}{}
		\begin{tabular}{c|c|c}

				\hline
				\multicolumn{2}{c|}{Approach} & CoNLL-2000 \\

				\hline
				\multicolumn{2}{c|}{Segmental RNN~\citep{kong2015segmental}} & $95.08$ \\
				\multicolumn{2}{c|}{Bi-LSTM + CRF~\citep{huang2015bidirectional}} & $94.46$ \\
				\multicolumn{2}{c|}{Char-IntNet-5~\citep{xin-etal-2018-learning}} & $95.29$ \\
				\multicolumn{2}{c|}{GCDT~\citep{liu-etal-2019-gcdt}} & $95.17$ \\

				\cdashline{1-3}
				\multicolumn{2}{c|}{Flair Embedding~\citep{akbik2018contextual}} &$96.72$  \\
				\multicolumn{2}{c|}{Cross-view Training~\citep{clark-etal-2018-semi}}  & $97.00$ \\
				\multicolumn{2}{c|}{GCDT w/ BERT~\citep{liu-etal-2019-gcdt}} & $96.81$ \\

				\hline
				\multirow{2}{*}{This Work} & Our Model & $\mathbf{96.13}$ \\
				& Our Model w/ BERT & $\mathbf{97.05}$ \\

				\hline

		\end{tabular}

		\caption{The performances of the baselines and our models on CoNLL-2000 dataset.}
		\label{tab:Experiment on Syntactic Chunking}
	\end{table*}

	\begin{table*}
		\centering

		\setlength{\tabcolsep}{7.1mm}{}
		\begin{tabular}{c|c|cc}

			\hline
			\multicolumn{2}{c|}{Approach} & PTB9 & UD1 \\

			\hline
			\multicolumn{2}{c|}{Segmental RNN~\citep{kong2015segmental}} & $92.16$ & $90.01$ \\
			\multicolumn{2}{c|}{Bi-RNN + CRF (single)~\citep{shao-etal-2017-character}} & $91.89$ & $89.41$ \\
			\multicolumn{2}{c|}{Bi-RNN + CRF (ensemble)~\citep{shao-etal-2017-character}} & $92.34$ & $89.75$ \\
			\multicolumn{2}{c|}{Lattice LSTM~\citep{zhang-yang-2018-chinese}} & $92.13$ & $90.09$ \\
			\multicolumn{2}{c|}{Glyce + Lattice LSTM~\citep{meng2019glyce}} & $92.38$ & $90.87$ \\

			\cdashline{1-4}
			\multicolumn{2}{c|}{BERT~\citep{devlin-etal-2019-bert}} & $92.29$ & $94.79$  \\
			\multicolumn{2}{c|}{Glyce + BERT~\citep{meng2019glyce}} & $93.15$ & $96.14$ \\

			\hline
			\multirow{2}{*}{This Work} & Our Model & $\mathbf{92.56}$ & $\mathbf{91.65}$ \\
			& Our Model w/ BERT & $\mathbf{93.38}$ & $\mathbf{96.43}$ \\

			\hline

		\end{tabular}
		\caption{The results on the two datasets of Chinese POS tagging.}
		\label{tab:Experiment on Chinese POS tagging}
	\end{table*}

	Extensive experiments have been conducted on syntactic chunking and Chinese POS tagging across 3 datasets. Firstly, our models have obtained new state-of-the-art performances on all the datasets. Then, we have investigated ablation studies to understand the importance of each component. Lastly, case study and the results on segmenting long-length sentences confirm the effectiveness of the proposed framework in capturing the long-term dependencies among segments.

\subsection{Settings}

	Syntactic chunking segments a word sequence into multiple labeled groups of words. We use CoNLL-2000 dataset~\citep{sang2000introduction}, which defines 11 syntactic chunk types (NP, VP, PP, etc.). Standard data includes a training set and a test set. Following~\citet{xin-etal-2018-learning}, we randomly sample 1000 sentences from the training set as the development set. Chinese POS tagging converts a Chinese character sequence into a token sequence and associates every word with a POS tag. We use Penn Chinese Treebank 9.0 (CTB9)~\citep{xue2005penn} and Universal Dependencies 1.4 (UD1)~\citep{nivre2016universal}. CTB9 contains the source text in various genres, covering its previous versions (e.g., CTB6). We use the Chinese section of UD1. We follow the same format and partition of the two datasets as~\citet{shao-etal-2017-character}.

	We use the same neural network configuration for all 3 datasets. The dimensions of token embedding and label embedding are respectively set as $300$ and $50$. The hidden unit sizes for the encoder and decoder are $256$ and $512$, respectively. The layers of two LSTMs are both $2$. L2 regularization is set as $1 \times 10^{-6}$ and dropout ratio is set as $0.4$ for reducing overfit. The above setting is obtained by grid search. We adopt Adam~\citep{kingma2014adam} as the optimization algorithm and adopt the suggested hyper-parameters. For CoNLL-2000 dataset, the cased, 300d Glove~\citep{pennington2014glove} is used to initialize token embedding. $\mathrm{CharCNN}$ is not used in Chinese tasks. The batch size is set as 16. All our models in experiments are running on NVIDIA Tesla P100.

	At test time, following previous literature, we convert the prediction of our model into IOB format and use the standard conlleval script\footnote{https://www.clips.uantwerpen.be/conll2000/chunking/co\\nlleval.txt.} to get the F1 score.  We select the model that works the best on development set, and then evaluate it on test set. In all the experiments, the improvements of our models over the baselines are statistically significant with $p < 0.05$ under t-test.

\subsection{Results on Syntactic Chunking}
\label{sec:Results on Syntactic Chunking}

	\begin{table*}
		\centering

		\begin{tabular}{c|ccc}

			\hline
			Approach & CoNLL-2000 & PTB9 & UD1 \\

			\hline
			Our Model & $\mathbf{96.13}$ & $\mathbf{92.56}$ & $\mathbf{91.65}$ \\

			\hline
			w/o $\mathrm{CharCNN}$ & $95.81$ & - & - \\
			w/o LSTM-minus, w/ Inside Algorithm & $\mathbf{96.35}$ & $\mathbf{92.71}$ & $\mathbf{91.87}$ \\

			\cdashline{1-4}
			w/o LSTM Decoder $f^d$, w/ MLP & $94.91$ & $91.28$ & $90.05$ \\

			w/o Phrase Representation $\mathbf{h}^p_{i_{k - 1}, j_{k-1}}$ in Equation \ref{equ:Equation 7} & $95.78$ & $92.09$ & $91.27$ \\

			w/o Label Representation $\mathbf{E}^l(l_{k-1})$ in Equation \ref{equ:Equation 7} & $95.82$ & $91.92$ & $91.04$ \\

			\cdashline{1-4}
			w/o Greedy Search, w/ Beam Search & $\mathbf{96.22}$ & $\mathbf{92.77}$ & $\mathbf{91.72}$  \\

			\hline
		\end{tabular}

		\caption{The results of ablation experiments on all three datasets.}
		\label{tab:Abalation Studies}
	\end{table*}

	Our models are compared with two groups of baselines. One of them is trained without any external resources besides the training data:
	\begin{description}

		\item[Segmental RNN] It's a segment-level model that defines a joint probability distribution over the partition of an input sequence and the labeling of the segments;

		\item[Bi-LSTM + CRF] It utilizes bidirectional LSTM to read the input sentence and CRF to decode the label sequence;

		\item[Char-IntNet-5] It is a funnel-shaped CNN model with no down-sampling which learns a better internal structure for tokens;

		\item[GCDT] It deepens the state transition path at each position in a sentence and assigns each token with a global representation learned from the entire sentence.

	\end{description}

	The other uses extra unlabeled corpora or fine-tunes on a pre-trained language model:
	\begin{description}

		\item[Flair Embedding] It firstly pre-trains a character-level language model on a large corpus, and then uses a sequence labeling model (e.g., Bi-LSTM + CRF) to fine-tune on it;

		\item[Cross-view Training] It designs a LSTM based sentence encoder to facilitate semi-supervised learning. The model can benefit from massive unlabeled corpora;

		\item[GCDT w/ BERT] It adopts BERT as additional token embeddings to improve GCDT.

	\end{description}

	We adopt most of the results of baselines as reported in~\citet{huang2015bidirectional,akbik2018contextual,xin-etal-2018-learning}. Since the evaluation method of GCDT is not standard (see the experiment setup in \citet{luo2020hierarchical}), we correct its source code\footnote{https://github.com/Adaxry/GCDT.} to retest the performance. The result for Segmental RNN is from our re-implementation.

	Table \ref{tab:Experiment on Syntactic Chunking} demonstrates that we have notably outperformed previous methods and achieved new state-of-the-art results on CoNLL-2000 dataset. When not using external resource, we obtain the F1 score of $96.13$, which outperforms Segmental RNN by $1.10\%$, Bi-LSTM + CRF by $1.77\%$, Char-IntNet-5 by $0.88\%$, and GCDT by $1.01\%$. Note that Segmental RNN, a segment-level model, underperforms Char-IntNet-5, a token-level model, by $0.22\%$. To make a fair comparison with the baselines using additional unlabeled corpora, we also use BERT~\citep{devlin-etal-2019-bert}, a powerful pretrained language model, to replace our token embedding. In this way, we achieve the F1 score of $97.05$, which outnumbers Flair Embedding by $0.34\%$, Cross-view Training by $0.05\%$, and GCDT w/ BERT by $0.25\%$. All these results verify the effectiveness of our framework.

\subsection{Results on Chinese POS Tagging}
\label{sec:Results on Chinese POS Tagging}

	We categorize the baselines into two types. The methods without using external resources:
	\begin{description}

		\item[Segmental RNN] Also described in Section \ref{sec:Results on Syntactic Chunking};

		\item[Bi-RNN + CRF] It utilizes bidirectional RNN and CRF to model joint word segmentation and POS tagging. Ensemble learning is also used to improve the results;

		\item[Lattice LSTM] It's a lattice-structured LSTM that encodes Chinese characters as well as all potential tokens that match a lexicon;

		\item[Glyce + Lattice LSTM] It incorporates Chinese glyph information into Lattice LSTM.

	\end{description}

	Others using a pre-trained language model:
	\begin{description}

		\item[BERT] It is a language model pre-training on a large corpus. It uses the representations from the last layer to predict IOB tags;

		\item[Glyce + BERT] It integrates Chinese glyph information into BERT Tagging model.

	\end{description}

	We take most of the performances of baselines from \citet{meng2019glyce}. The results for Segmental RNN are from our re-implementation.

	Table \ref{tab:Experiment on Chinese POS tagging} shows that our models have achieved state-of-the-art results on the two datasets, PTB9 and UD1. When BERT is not used, we have obtained the F1 scores of $92.56$ and $91.65$, which outperform Glyce + Lattice LSTM by $0.19\%$ and $0.86\%$ and Bi-RNN + CRF (ensemble) by $0.24\%$ and $2.12\%$. Note that Segmental RNN, a segment-level model, underperforms Lattice LSTM, a token-level model, by $0.09\%$ on UD1. When using BERT, even without incorporating Chinese glyph information, we still obtain the F1 scores of $93.38$ and $96.43$, which outperform Glyce + BERT by $0.25\%$ and $0.30\%$. These results further confirm the effectiveness of our proposed framework.

\subsection{Ablation Studies}

	\begin{table*}
		\centering

		\begin{tabular}{c|cccc|c}

			\hline
			Approach &  1-22 (711) & 23-44 (1030) & 45-66 (248) & 67-88 (23) & Overall \\

			\hline
			Bi-LSTM + CRF & $94.01$ & $94.76$ & $92.98$ & $87.09$ & 94.23 \\

			GCDT & $94.95$ & $95.52$ & $93.77$ & $87.14$ & 95.17 \\

			\hline
			Our Model & $\mathbf{96.01}$ & $\mathbf{96.55}$ & $\mathbf{94.82}$ & $\mathbf{91.07}$ & $\mathbf{96.13}$ \\

			\hline
		\end{tabular}

		\caption{The F1 scores for the sentences of different length ranges.}
		\label{tab:Segmenting Long-length Sentences}
	\end{table*}

	\begin{table*}
		\centering

		\setlength{\tabcolsep}{1.1mm}{}
		\begin{tabular}{c|c|c}

			\hline
			\multicolumn{2}{c|}{\multirow{2}{*}{Input Sentence}} & Other antibodies sparked by the preparation are of a sort \\
			\multicolumn{2}{c|}{} & rarely present in large quantities in infected or ill individuals \\

			\hline
			\multirow{6}{*}{Output Segments} & \multirow{3}{*}{GCDT} & (Other antibodies, NP) (sparked, VP) (by, PP) (the preparation, NP) \\

			& & (are, VP) (of, PP) (a sort, NP) \textbf{(rarely present, VP)} (in, PP)\\

			& & (large quantities, NP) (in, PP) (infected or ill individuals, NP) \\

			\cline{2-3}
			& \multirow{3}{*}{Our Model} & (Other antibodies, NP) (sparked, VP) (by, PP) (the preparation, NP) \\

			& & (are, VP) (of, PP) (a sort, NP) (rarely present, ADJP) (in, PP)\\

			& & (large quantities, NP) (in, PP) (infected or ill individuals, NP) \\

			\hline

		\end{tabular}
		\caption{The case is from CoNLL-2000 dataset. The predicted segments of our model is correct.}
		\label{tab:Case Study}
	\end{table*}

	As shown in Table \ref{tab:Abalation Studies}, we conduct ablation studies to explore the impact of every component.

	\paragraph{Effect of Representation Learning.} Following prior works~\citep{ma-hovy-2016-end,liu-etal-2019-gcdt}, we employ $\mathrm{CharCNN}$ to incorporate character information into word representations. By removing it, the F1 score on CoNLL-2000 decreases by $0.33\%$. Inside algorithm is another technique to construct phrase representations. After using it to replace LSTM-minus, the results on the three datasets are slightly improved by $0.23\%$, $0.16\%$, and $0.24\%$. Our implementation of inside algorithm is the same as described in \citet{drozdov-etal-2019-unsupervised-latent}. Despite the slight improvements, its time complexity $\mathcal{O}(n^3)$ is too slow for both training and inference. Empirically, we find that the running time of inside algorithm is about $7$ times slower than that of LSTM-minus.

	\paragraph{Effect  of Modeling the Dependencies Among Segments.} LSTM decoder $f^d$ models the long-term dependencies among segments. The prediction of every segment $y_k$ is conditional on prior segments $y_{k'}, 1 \le k' < k$. By replacing it with multilayer perceptron (MLP), the performances fall by $1.29\%$, $1.40\%$, and $1.78\%$ on the three datasets. We use both phrase representation $\mathbf{h}^p_{i_{k - 1}, j_{k-1}}$ and label representation $\mathbf{E}^l(l_{k-1})$ to embed the previous segment $y_{k-1}$. After removing the phrase representations, the F1 scores decrease by $0.37\%$, $0.51\%$, and $0.42\%$ on the three datasets. By removing the label representations, the results also drop by $0.32\%$, $0.70\%$, and $0.67\%$.

	\paragraph{Effect of Inference Algorithm} Beam search is a widely used technique in language generation tasks, like machine translation~\citep{bahdanau2014neural,NIPS2017Vaswani,li2020rewriter} and data-to-text generation~\citep{shen-etal-2020-neural,li-etal-2020-slot,li2020interpretable}. We have attempted to use beam search (with the beam size being $5$), instead of greedy search, for inference. By doing so, the F1 scores of our models increased by $0.09\%$ on CoNLL-2000, $0.23\%$ on PTB9, and $0.08\%$ on UD1. While obtaining slightly better performances, the time costs for inference become intolerably high (increase about $10$ times). Therefore, we adopt greedy search as the default inference algorithm.

\subsection{Segmenting Long-length Sentences}
\label{sec:Segmenting Long-length Sentences}

	\begin{table*}
		\centering

		\begin{tabular}{c|c|cc}

			\hline
			Approach &  Time Complexity & Training Time & Evaluation Time \\

			\hline
			Segmental RNN & $\mathcal{O}(n^2|\mathcal{L}|^2)$ &16m17s & 3m01s \\

			Bi-LSTM + CRF & $\mathcal{O}(n|\mathcal{L}|^2)$ & 3m28s & 0m26s \\

			\hline
			Our Model & $\mathcal{O}(n|\mathcal{L}|)$ & 2m36s & 0m20s \\

			\hline
		\end{tabular}

		\caption{Running time comparisons on CoNLL-2000 dataset.}
		\label{tab:Running Time Comparisons}
	\end{table*}

	\begin{table*}
		\centering

		\begin{tabular}{c|c|c}

			\hline
			Approach & CoNLL-2003 & OntoNotes 5.0 \\

			\hline
			BiLSTM-CNN-CRF~\citep{chen2019grn} & $91.21$ & $87.05$ \\

			GRN~\citep{chen2019grn} & $91.44$ & $87.67$ \\

			HCR~\citep{luo2020hierarchical}  & $\mathbf{91.96}$ & $\mathbf{87.98}$ \\

			\hline
			Our Model & $91.42$ & $87.74$ \\

			\hline
		\end{tabular}
		\caption{The results on two NER datasets.}
		\label{tab:Experiment on NER}
	\end{table*}

	Compared with token-level methods, our framework better captures the long-term dependencies among segments. Therefore, our model should be more accurate in segmenting long-length sentences. To verify this, we test the baselines and our model on the sentences of different lengths.

	The study is conducted on CoNLL-2000 dataset. Bi-LSTM + CRF and GCDT are very strong baselines and have open source implementations. We use toolkit NCRFPP\footnote{https://github.com/jiesutd/NCRFpp.} to reproduce the performances of Bi-LSTM + CRF. We use the reproduction in Section \ref{sec:Results on Syntactic Chunking} as the results of GCDT. Table \ref{tab:Segmenting Long-length Sentences} shows the experiment results. Each column name denotes the sentence length range and the case number. Our model notably outperforms prior methods in terms of long sentence length ranges. For length range 45-66, our model obtains the F1 score of $94.82$, which outperforms Bi-LSTM + CRF by $1.98\%$ and GCDT by $1.12\%$. For length range 67-88, our model outperforms Bi-LSTM + CRF by $4.57\%$ and GCDT by $4.51\%$.

\subsection{Case Study}

	In Table \ref{tab:Case Study}, we present an example extracted from the test set of CoNLL-2000. Given an input sentence, we show the prediction from a strong baseline, GCDT, and our model. The output segments from our model are consistent with the ground truth segments. The segment in bold is the one incorrectly produced by GCDT.

	Understanding the structure of this sentence is very hard because the main constituents ``Other antibodies", ``are", and ``rarely present" locate far apart. By removing the two middle phrases ``sparked by the preparation" and ``of a sort", the original sentence can be simplified to ``Other antibodies are rarely present in large quantities in infected or ill individuals", which is very clear. Therefore, correctly predicting the segment, (rarely present, ADJP), implies that our model well captures the long-term dependencies among the segments. For GCDT, a token-level segmentation model, it mistakes ``rarely present" for the verb phrase of ``a sort". This may result from the following two causes:
	\begin{itemize}

		\item The adjacent segments labeled with (NP, VP) frequently appear in training data;

		\item the token ``present" acting as a verb is much more common than as a adjective.

	\end{itemize}
	Both potentials indicate that token-level models can't capture the long-term dependencies well.

\subsection{Running Time Analysis}
\label{sec:Running Time Analysis}

	Table \ref{tab:Running Time Comparisons} demonstrates the running time comparison among different methods. The last two columns are respectively the running times for training (one epoch) and evaluation. We set the batch size as 16 and run all the models on 1 GPU. From the table, we can draw the following two conclusions. Firstly, Segmental RNN, a segment-level model, is very slow for both training and inference due to the high time complexity. For instance, its training and testing are respectively $6.26$ and $9.05$ times slower than ours. Secondly, our framework is very efficient. For example, training our model for one epoch is $1.33$ times faster than training Bi-LSTM + CRF, a token-level model.

\subsection{Results on NER}

	While our model is tailored for sequence segmentation tasks, we have also tested its performances on two widely used NER datasets, CoNLL-2003~\citep{sang2003introduction} and OntoNotes 5.0~\citep{pradhan2013towards}. Note that, in NER, the label correlations among adjacent segments are very weak. This seems bad for our model.

	Table \ref{tab:Experiment on NER} diagrams the comparison of our model and strong NER baselines. The results of GRN and HCR are copied from \citet{chen2019grn,luo2020hierarchical}. For BiLSTM-CNN-CRF, its scores on CoNLL-2003 and OntoNotes 5.0 are respectively from \citet{chen2019grn} and our re-implementation. From the table, we can see that our model is competitive with prior methods. For example, our model underperforms HCR by only $0.27\%$ on OntoNotes 5.0. In particular, our F1 scores are notably higher than those of BiLSTM-CNN-CRF by $0.23\%$ and $0.79\%$ on the two datasets. This experiment shows that our model is applicable to general sequence labeling tasks.

\section{Related Work}

	There are two mainstream approaches to sequence segmentation. One of them treat sequence segmentation as a sequence labeling problem by using IOB tagging scheme~\citep{huang2015bidirectional,xin-etal-2018-learning,clark-etal-2018-semi,akbik2018contextual,liu-etal-2019-gcdt,li-etal-2020-handling}. Each token in a sentence is labeled as B-tag if it’s the beginning of a segment, I-tag if it is inside but not the first one within the segment, or O otherwise. This method is extensively studied by prior works and provides tons of state-of-the-art results.  \newcite{xin-etal-2018-learning} introduce a funnel-shaped convolutional architecture that learns a better internal structure for the tokens. \newcite{akbik2018contextual} propose an efficient character-level framework that uses pretrained character embedding. \newcite{clark-etal-2018-semi} adopt semi-supervised learning to train LSTM encoder using both labeled and unlabeled corpora. Despite the effectiveness, these models are at token level, relying on multiple token-level labels to represent a single segment. This limits their full potential to capture the long-term dependencies among segments.

	The other uses transition-based systems as the backbone to incrementally segment and label an input sequence~\citep{qian2015transition,zhang2016transition,zhang2018simple}. For example, \newcite{qian2015transition} design special transition actions to jointly segment, tag, and normalize a sentence. These models have many advantages, such as theoretically lower time complexity and capturing non-local features. However, they are still token-level models, which predict transition actions to shift a token from the buffer to the stack or assign a label to a span.

	Recently, there is a surge of interest in developing span-based models, such as Segmental RNN~\citep{sarawagi2004semi,kong2015segmental} and LUA~\citep{li2020segmenting}. These methods circumvent using token-level labels and directly label the phrases in a sentence. Span-based models also enjoy great popularity in language modeling~\citep{li2020span}, NER~\citep{yu-etal-2020-named,li2021empirical}, and constituent parsing~\citep{cross-huang-2016-span,stern-etal-2017-minimal}. However, its underperforms current token-level models (see Section \ref{sec:Results on Syntactic Chunking} and Section \ref{sec:Results on Chinese POS Tagging}) and is very slow in terms of running time (see Section \ref{sec:Running Time Analysis}).

\section{Conclusion}

	In this work, we present a novel framework to sequence segmentation that segments a sentence at segment level. The segmentation process is iterative and incremental. For every step, it determines the leftmost segment of the remaining sequence. Implementations involve LSTM-minus to extract phrase representations and RNN to model the iterations of leftmost segment determination. Extensive experiments have been conducted on syntactic chunking and Chinese POS tagging across 3 datasets. We have achieved new state-of-the-art performances on all of them. Case study and the results on segmenting long-length sentences both verify the effectiveness of our framework in modeling long-term dependencies.

\section*{Acknowledgments}

	 This work was done when the first author did internship at Ant Group. We thank anonymous reviewers for their kind and constructive suggestions.


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

Magnetic fields dominate the dynamics in the Sun's upper atmosphere,
the solar corona. On the solar surface, i.e., the photosphere,
magnetic fields are seen to change continuously; magnetic flux
emergence brings new flux from the solar interior into the atmosphere,
and meanwhile the flux is advected and dispersed by surface motions
such as granulation, differential rotation and meridional
circulation. Consequently, the coronal field evolves in response to
(or driven by) the changing of the photospheric field, and thus
complex dynamics occur ubiquitously in the corona, including the
interaction of newly emerging field with the pre-existing one,
twisting and shearing of the magnetic arcade fields, magnetic
reconnection, and magnetic explosions which are manifested as flares
and coronal mass ejections (CMEs).

As we are still not able to measure directly the three-dimensional
(3D) coronal magnetic fields, numerical modelling has long been
employed to reconstruct or simulate the coronal magnetic fields based
on different assumption of the magnetohydrodynamic (MHD) equations,
such as the most-frequently used force-free field
model~\citep{Wiegelmann2012solar}, which has been developed for
over four decades. However, the force-free field
assumption is only valid for the equilibrium of the corona, and it
cannot be used to follow a continual and dynamic evolution of the
magnetic fields. In recent years, with vector magnetogram data in the photosphere
measured routinely with high resolution and cadence, data-driven
modelling is becoming a viable tool to study coronal magnetic field
evolution, which can self-consistently describe both the quasi-static
and the dynamic evolution phases~\citep[e.g.,][]{Wu2006, FengX2012,
  Cheung2012, Jiang2016NC, Leake2017, Inoue2018, Hayashi2019,
  GuoY2019, Pomoell2019}. Still, due to the limited constraint from
observation, data-driven models are developed with very different
settings from each other~\citep{Toriumi2020}. For simplicity, some used the
magneto-frictional model~\citep{Cheung2012, Pomoell2019}, in which the
Lorentz force is balanced by a fictional plasma friction force. As
such, the magnetic field evolves mainly in a quasi-static way~\citep{Yang1986},
although in some case an eruption can be reproduced but it evolves in
a much slower rate than the realistic one~\citep{Pomoell2019}, because the dynamic is
strongly reduced by the frictional force. Some used the so-called
zero-$\beta$ model~\citep{Inoue2016Review, Inoue2018, GuoY2019}, in which the gas
pressure and gravity are neglected. The zero-$\beta$ model might fail
when there is fast reconnection in the field, in which the thermal
pressure could play an important role in the dynamics in weak-field region of magnetic
field dissipation. Therefore, it is more realistic and {accurate} to
solve the full MHD equations to deal with the nonlinear interaction of magnetic fields with plasma.

To drive the full MHD model, one needs to specify all the eight
variables (namely plasma density, temperature, and three components of
velocity and magnetic field, respectively) in a self-consistent way at
the lower boundary. Previously, with only the magnetic field obtained
from observations, we employed the projected-characteristic
method~\citep{Nakagawa1980, Nakagawa1987, Wu1987, Hayashi2005, Wu2006}
to specify the other variables according to the information of
characteristics based on the wave-decomposition principle of the full
MHD system~\citep[e.g.,][]{Jiang2016ApJ,Jiang2016NC}. Specifically,
the full MHD equations are a hyperbolic system that can be
eigen-decomposed into a set of characteristic wave equations (i.e.,
compatibility equations), which are independent with each other; on
the boundary, these waves may propagate inward or outward of the
computational domain because the wave speeds (the eigenvalues) of both
signs generally exist at the boundary. If the wave goes out of the
computational domain through the boundary, it carries information
from the inner grid to the boundary, and thus the corresponding compatibility
equation should be used to constrain the variables on the boundary
surface. Thus if there are five waves going out, and with three
components of magnetic field specified by observed data, the eight
variables are fully determined. Ideally and in principle, the
projected-characteristic method is the most self-consistent one with
inputting data that is partially available (for example, in our case,
only the data of magnetic field).  However, such conditions, i.e.,
exactly five waves going out, are often not satisfied in the whole
boundary and other assumptions are still necessary. Even worse, in
areas near the magnetic polarity inversion line (PIL) where the normal magnetic field component is small,
the Alfv{\'e}n wave information goes mainly in the transverse
direction rather than the normal one, and the projected-characteristic
approach may fail. Another shortcoming of the method lies in its
difficulty in code implementation, which needs to perform
eigen-decomposition and solve a linear system of the compatibility
equations on every grid point on the boundary to recover the primitive
variables from the characteristic ones.

In this paper we test another way of specifying the bottom boundary
conditions, which uses the surface velocity derived by the DAVE4VM
method~\citep{Schuck2008} to drive the MHD model. The differential
affine velocity estimator (DAVE) was first developed for estimating
velocities from line-of-sight magnetograms, and was then modified to
directly incorporate horizontal magnetic fields to produce a
differential affine velocity estimator for vector magnetograms
(thus called DAVE4VM). It is generally accepted that the coronal magnetic fields are
line tied in the dense photosphere and are advected passively by the surface
motions of the photosphere, such as the shearing, rotational, and converging flows~\citep{Priest2002}.
Thus, with the surface velocity in hand, we can solve the
magnetic induction equation on the bottom boundary to update
self-consistently the magnetic field to implement such a line-tying boundary condition,
which is a much simpler way than
using the projected-characteristic method. To illustrate the approach,
we take the solar AR 12158 as an example to simulate its two-day
evolution from 2014 September 8 to 10 during its passage of the solar
disk. This AR is of interest since it produced an X1.6 eruptive flare
accompanied with a fast CME with speed of $\sim 1300$~km~s$^{-1}$,
which is well documented in the literature~\citep{Vemareddy2016,
  ZhaoJ2016, ZhouG2016}. The AR is well isolated from neighboring ARs,
thus is suited for our modelling focused on a single AR. Within a few days prior to this
major eruption, the AR developed from a weakly sheared magnetic arcade
into a distinct sigmoidal configuration, indicating a
continual injection of non-potential magnetic energy into the corona through the
photospheric surface motion. Indeed, prior to the eruption, the major
sunspot of the AR showed a significant rotation of over $200^{\circ}$
in five days~\citep{Vemareddy2016}. Based on the sigmoidal hot coronal loop
seen immediately before the flare, many authors have interpreted it as a
pre-existing magnetic flux rope which erupted and resulted in the flare
and CME~\citep[e.g.,][]{Vemareddy2016, ZhaoJ2016, ZhouG2016}, while
some nonlinear force-free field extrapolations appear not to support
this~\citep{DuanA2017, Lee2018}. Thus a fully data-driven MHD
simulation can provide valuable insight in addressing this issue by
following the dynamic evolution of the coronal magnetic field,
although the main objective of this paper is to describe the methods.

In the following we first describe our model equation in
Section~\ref{sec:model}. The data-driven simulation consists of three
steps, and first we constructed an MHD equilibrium based on a single
vector magnetogram observed for the start time of our simulation,
which is described in Section~\ref{initial}. Then in
Section~\ref{flow} we calculated the surface flow field using the
DAVE4VM code with the time series of vector magnetograms. Finally, we
input the flow field in the model to drive the evolution of the MHD
system, as described in Section~\ref{evol}. Summary and discussion are
given in Section~\ref{concl}.

\begin{table}[htbp]
  \centering
  \begin{tabular}{lll}
    \hline \hline
    Variable &  Expression & Value \\
    \hline
    Density        & $\rho_{s}$ =  $n m $ & $2.29 \times 10^{-15}$~g~cm$^{-3}$ \\
    Temperature    & $T_{s}$       & $1\times 10^{6}$~K \\
    Length         & $L_{s}$  = $16$~arcsec  &  $11.52$~Mm\\
    Pressure       & $p_{s}=2nk_{B}T_{s}$   &  $2.76\times 10^{-2}$~Pa \\
    Magnetic field & $B_{s}=\sqrt{\mu_{0}p_{s}}$ & $1.86$~G\\
    Velocity       & $v_{s}=\sqrt{p_s / \rho_{s} }$ & $110$~km~s$^{-1}$ \\
    Time           & $t_{s}=L_{s}/v_{s}$ & $105$~s \\
    Gravity        & $g_{s}=v_{s}/t_{s}$ & $1.05$~km~s$^{-2}$\\
    \hline
  \end{tabular}
  \caption{Parameters used for non-dimensionalization. $n$ is a
    typical value of electron number density in the corona given by
    $n=1\times 10^{9}$~cm$^{-3}$ and $m$ is the mean atomic mass.}
  \label{tab:table1}
\end{table}

\section{MHD Equations}
\label{sec:model}

We numerically solve the full MHD equations in a 3D Cartesian geometry
by an advanced conservation element and solution element (CESE)
method~\citep{Jiang2010}.  Before describing the model equations in
the code, it is necessary to specify the quantities used for
non-dimensionalization. Here we use typical values at the base of the
corona for non-dimensionalization as given in
Table~\ref{tab:table1}. In the rest of the paper all the variables and
quantities are written in non-dimensionalized form if not mentioned
specially.

In non-dimensionalized form, the full set of MHD equations are given as
\begin{eqnarray}
  \label{eq:MHD}
  \frac{\partial \rho}{\partial t}+\nabla\cdot (\rho\vec v) =
  -\nu_{\rho}(\rho-\rho_0),\nonumber \\
  \rho\frac{D\mathbf{v}}{D t} = - {\bf \nabla }  p+\vec J\times \vec B+\rho\vec
  g + \nabla\cdot(\nu\rho\nabla\mathbf{v}),\nonumber\\
  \frac{\partial \vec B}{\partial t} =
   {\bf \nabla }  \times (\vec v \times \vec B), \nonumber\\
  \frac{\partial T}{\partial t}+\nabla\cdot (T\vec v) =
  (2-\gamma)T\nabla\cdot\vec v.
\end{eqnarray}
where $\vec J = \nabla \times \vec B$, $\nu$ is the kinetic viscosity,
and $\gamma$ is the adiabatic index.

Note that we artificially add a source term $-\nu_{\rho}(\rho-\rho_0)$
to the continuity equation, where $\rho_0$ is the density at the
initial time $t=0$ (or some prescribed form), and $\nu_{\rho}$ is a
prescribed coefficient. This term is used to avoid a ever-decreasing
of the density in the strong magnetic field region, which we often encounter
in the very low-$\beta$ simulation. It can maintain
the maximum Alfv{\'e}n speed in a reasonable level, which may
otherwise increase and make the iteration time step smaller and
smaller and the long-term simulation unmanageable. {Specifically, this source
term is a Newton relaxation of the density to its initial value by a
time scale of
\begin{equation}
 \tau_{\rho} = \frac{1}{\nu_{\rho}} = 20 \tau_{\rm A},
\end{equation}
where $\tau_{\rm A} = 1/v_{\rm A}$ is the Alfv{\'e}n time with length of 1 (the length unit) and the Alfv{\'e}n speed $v_{\rm A} = B/\sqrt{\rho}$. Thus it is
sufficiently large to avoid its influence on the fast dynamics of
Alfv{\'e}nic time scales. As a result, we used $ \nu_{\rho} = 0.05 v_{\rm A} $ in all the
simulation in this paper.}

No explicit resistivity is used in the magnetic induction equation,
but magnetic reconnection is still allowed through numerical diffusion
when a current layer is sufficiently narrow such that its width is
close to the grid resolution. In the energy (or temperature) equation,
we set $\gamma=1$ for simplicity, such that the energy equation
describes an isothermal process. \footnote{{Although in this case we can simply discard the energy equation by setting the temperature as a constant, we still keep the full set of equations in our code which can thus describe either the isothermal or adiabatic process by choose different value of $\gamma$.}} The kinetic viscosity $\nu$ will be given
with different values when needed, which is described in the
following sections.

\begin{figure*}[htbp]
  \centering
 \includegraphics[width=0.8\textwidth]{fig1}
 \caption{Comparison of the original HMI vector magnetogram
   for time of 00:00~UT on 2014 September 8 (top
   panels), the preprocessed one (middle) and the final smoothed
   one (bottom). From left to right are shown for magnetic field
   components $B_x$, $B_y$, and $B_z$, respectively.}
  \label{1}
\end{figure*}

\begin{figure}[htbp]
  \centering
 \includegraphics[width=0.45\textwidth]{fig2}
 \caption{Evolutions of the total magnetic energy, kinetic energy, and
   residual of magnetic field in the process of the constructing the
   initial MHD equilibrium. The black curves show the results for the
   first evolution phase, and the blue ones show that for the second,
   `deeper' evolution phase.}
  \label{2}
\end{figure}

\begin{figure*}[htbp]
  \centering
 \includegraphics[width=0.7\textwidth]{fig3}
 \caption{3D magnetic field of the MHD equilibrium based on the
   magnetogram taken for 00:00~UT on 2014 September 8. The thick lines
   are magnetic field lines and are pseudo-colored by value of the
   force-free parameter $\alpha$. The background is shown with the
   magnetic flux distribution (i.e., $B_z$) on the bottom surface of
   the model.}
  \label{3}
\end{figure*}

\section{Construction of an initial MHD equilibrium}
\label{initial}
We first constructed an initial MHD equilibrium based on a single
vector magnetogram taken for time of 00:00~UT on 2014 September 8 by
SDO/HMI. Such an equilibrium is assumed to exist when the corona is
not in the eruptive stage, and is crucial for starting our subsequent
data-driven evolution. The vector magnetogram is preprocessed to
{reduce} the Lorentz force and is further smoothed to filter out the
small-scale structures that are not sufficiently resolved in our
simulation. Reduction of the Lorentz force is helpful for reaching a
more force-free equilibrium state~\citep[e.g.,][]{Wiegelmann2004}, and
smoothing is also necessary to mimic the magnetic field at the coronal
base rather directly the photosphere, because the lower boundary of
our model is placed at the base of the corona~\citep{JiangC2020}. The
preprocessing is done using a code developed by~\citet{Jiang2014Prep},
which is originally designed for NLFFF extrapolation using the HMI
data~\citep{Jiang2013NLFFF}. {Specifically, the total Lorentz force and torque are quantified by two surface integrals associated with three components of the photospheric magnetic field, and an optimization method is employed to minimize these two functions by modifying the magnetic field components within margins of measurement errors.}
Then Gaussian smoothing with FWHM of
8~arcsec is applied to all the three components of the magnetic
field. {Figure}~\ref{1} compares the original HMI vector magnetogram with
the preprocessed one as well as the final smoothed version. Using the
vertical component $B_z$ of this preprocessed and smoothed
magnetogram, a potential field is extrapolated and used as the initial
condition of the magnetic field in the MHD model.

In addition to the magnetic field, an initial plasma as the
background atmosphere is also needed to start the MHD simulation. We
used an isothermal plasma in hydrostatic equilibrium. It is stratified
by solar gravity with a density $\rho=1$ at the bottom and a uniform
temperature of $T=1$. Using this typical coronal plasma, we did not
directly input the magnetic field into the model but multiplies it
with a factor of $0.05$ such that the maximum of magnetic field
strength in normalized value is approximately $50\sim 100$ in the
model. If using the original values of magnetic field, its strength
(and the characteristic Alfv{\'e}n speed) near the lower surface is
too larger, which will be a too heavy burden on computation as the
time step of our simulation is limited by the CFL condition. With the
reduced magnetic field, we further modified the value of solar gravity
to avoid an unrealistic large plasma $\beta$ (and small Alfv{\'e}n speed) in the corona. This is
because if using the real number of the solar gravity
($g_{\odot} = 274$~m~s$^{-2}/g_{s} = 0.26$), it results in a pressure
scale height of $H_{p} = 3.8$, by which the plasma pressure and
density decay with height much slower than the magnetic field. With
the weak magnetic field strength we used, the plasma $\beta$ will
increase with height very fast to above $1$, which is not realistic in
the corona. To make the pressure (and density) decrease faster in the
lower corona, we modified the gravity as
\begin{equation}
  \label{eq:gravity}
  g = \frac{1.5H_{p}g_{\odot}}{(1 + 0.15z)^{2}}.
\end{equation}
By this, we got a plasma with $\beta<1$ mainly within $z<10$ and the
smallest value is $5\times 10^{-3}$.

Then we input the transverse field of the smoothed magnetogram to the
model. This is done by modifying the transverse field on the bottom
boundary incrementally using linear extrapolation from the potential
field to the vector magetogram in time with a duration of 1 $t_s$
until it matches the vector magnetogram. This will drive the coronal
magnetic field to evolve away from the initial potential state, since
the change of the transverse field will inject electric currents and
thus Lorentz forces, which drive motions in the computational
volume. Note that in this phase all other variables on the bottom
boundary are simply fixed, thus the velocity remaining zero. This is
somewhat un-physical since the Lorentz force will introduce nonzero
flows on the bottom boundary, but it provides a simple and `safe' way
(avoiding numerical instability) to bring the transverse magnetic
field into the model. Once the magnetic field on the bottom surface is
identical to that of the vector magnetogram, the system is then left
to relax to equilibrium with all the variables (including the magnetic
field) on the bottom boundary fixed. To avoid a too large velocity in
this phase such that the system can relax faster, we set the kinetic
viscosity coefficient as $\nu = 0.5\Delta x^{2}/\Delta t$, where
$\Delta x$ is the local grid spacing and $\Delta t$ the local time
step, determined by the CFL condition with the fastest magnetosonic
speed. Actually this is the largest viscosity one can use with given
grid size $\Delta x$ and time step $\Delta t$, because the CFL
condition for a purely diffusive equation with diffusion coefficient
$\nu$ requires $\Delta t \le 0.5 \Delta x^{2}/\nu$.

For the purpose of minimizing the numerical boundary influences
introduced by the side and top boundaries of the computational volume,
we used a sufficiently large box of
$(-32, -32, 0) < (x, y, z) < (32, 32, 64)$ embedding the field of
view of the magnetogram of $(-8.75, -8.25) < (x, y) < (8.75, 8.25)$.
The full computational volume is resolve by a non-uniform
block-structured grid with adaptive mesh refinement (AMR), in which
the highest and lowest resolution are
$\Delta x = \Delta y = \Delta z = 1/16$ (corresponding to $1$~arcsec
or 720~km, matching the resolution of the vector magnetogram) and
$1/2$, respectively. The AMR is controlled to resolve with the
smallest grids the regions of strong magnetic gradients and current
density. {The magnetic field outside of the area of the magnetograms on the lower boundary is given as zero for the vertical component and simply fixed as the potential field for the transverse components.} On the side and top boundaries, we fixed the plasma density,
temperature and velocity. The horizontal components of magnetic field
are linearly extrapolated from the inner points, while the normal
component is modified according to the divergence-free condition to
avoid any numerical magnetic divergence accumulated on the boundaries.

In {Figure}~\ref{2}, the curves colored in black show the evolution of the
magnetic and kinetic energies integrated for the computational volume,
and the residual of the magnetic field of two consecutive time
steps which is defined as
\begin{equation}
  {\rm Res}_B = \sqrt{ \frac{1}{3} \sum_{\delta=x,y,z} \frac{\sum\nolimits_{i}
  \left( B_{i\delta}^{k}-B_{i\delta}^{k-1}\right) ^{2} }{\sum\nolimits_{i}\left( B_{i\delta}^{k} \right) ^{2} } },
\end{equation}
where the indices $k$ and $k-1$ refer to the two consecutive time
steps and $i$ goes through all the mesh points. It can be seen that
the magnetic energy increases sharply in a few time
units\footnote{Note that at the very beginning the magnetic energy
  actually decreases shortly for about $0.5~t_s$, which is unphysical
  as the potential field energy is in principle the lowest energy with
  a given magnetic flux distribution on the bottom boundary. Such an
  unphysical evolution is a result of the fact that we modified
  directly the transverse magnetic field on the bottom boundary in an
  unphsyical way.}, reaching $\sim 1.4$ of the potential field energy
$E_0$ (here $E_{0} = 1.2 \times 10^{33}$~erg when scaled to the
realistic value in the corona), and then keeps almost constant during
the relaxing phase with bottom boundary fixed. Very similar, the
kinetic energy first increases and later keeps on the level of
$3\times 10^{-3}$ of $E_0$. The residual of the magnetic field
increases in the first $t_s$ as we continually modified the transverse
field at the bottom boundary which drives quickly the evolution of the
field in the corona. Then it decreases to below $10^{-5}$ with a time
duration of $20~t_s$, which indicates that the magnetic field reaches
a quasi-equilibrium state.

To make the field even closer to equilibrium, we carried out a
`deeper' relaxation by running the model again, which is started with
the relaxed magnetic field obtained at $t=20~t_s$ and the initially
hydrostatic plasma. Now we reduce the kinetic viscosity to
$\nu = 0.05\Delta x^{2}/\Delta t$, i.e., an order of magnitude smaller
than the previously used one, which will let the magnetic field relax
further. Furthermore, the magnetic field at the boundary boundary is
allowed to evolve in a self-consistent way by assuming the bottom
boundary as a perfectly line-tying and fixed (i.e., $\vec v = 0$)
surface of magnetic field lines. However, such a line-tying condition
does not indicate that all magnetic field components on the
boundary are fixed, because even though the velocity $\vec v$ is given
as zero on the bottom boundary, it is not necessarily zero in the
neighboring inner points. To self-consistently update the magnetic
field, we solve the magnetic induction equation on the bottom
boundary. Slightly different from the one in the main
equations~(\ref{eq:MHD}), the induction equation at the bottom surface
is given as
\begin{equation}
\label{photo_B_equ}
  \frac{\partial \vec B}{\partial t} =   {\bf \nabla }  \times (\vec v \times \vec B) +
  \eta_{\rm stable} \nabla_{\perp}^{2} \vec B,
\end{equation}
where we added a surface diffusion term defined by using a surface
Laplace operator as
$\nabla_{\perp}^{2} =\frac{\partial^{2}}{\partial
  x^{2}}+\frac{\partial^{2}}{\partial y^{2}}$
with a small resistivity for numerical stability near the PIL
$\eta_{\rm stable} = 1\times 10^{-3} e^{-B_z^2}$, since the magnetic
field often has the strongest gradient on the photosphere {around} the
main PIL. The surface induction {Equation}~(\ref{photo_B_equ}) in the code is
realized by second-order difference in space and first-order forward
difference in time. Specifically, on the bottom boundary (we do not
use ghost cell), we first compute $\vec v\times \vec B$, and then use
central difference in horizontal direction and one-sided difference
(also 2nd order) in the vertical direction to compute the convection
term $ {\bf \nabla }  \times (\vec v \times \vec B)$. The surface Laplace
operator is also realized by central difference.

The curves colored in blue in {Figure}~\ref{2} show the evolution of
different parameters during this relaxation phase~\footnote{{Note that the initial values of the blue curves (the deeper relaxation phase) do not equal to the values of the black curves at $t = 20 t_s$ although the deeper relaxation phase starts from that time point. This is because in the deeper relaxation process, there are three ways different from the initial relaxation one: (1) the velocity is reset to zero and the plasma is reset to hydrostatic state, thus the kinetic energy is reset to zero; (2) the viscosity is abruptly reduced by an order of magnitude; (3) the boundary condition is changed. Furthermore, not data of every time step in the run is recorded, thus small difference is shown in the magnetic energy of the initial value of the blue line (which is not exactly for the time of $t=20 t_s$) from the $t = 20 t_s$ of the black line.}}. Initially one can
see a fast decrease of the magnetic energy because the magnetic field
becomes more relaxed. As the viscosity is reduced significantly, the
kinetic energy first increases, as driven by the residual Lorentz
force of the magnetic field, to almost $5\times 10^{-3}~E_0$. Then as
the magnetic field relaxed, the kinetic energy decreases fast to
eventually less than $10^{-3}~E_0$, which is a very low level. The
residual of magnetic field of two consecutive time steps also
decreases to $10^{-5}$. These values show that the magnetic field
reaches an excellent equilibrium.

{Figure}~\ref{3} shows the 3D magnetic field lines of the final relaxed
MHD equilibrium. Note that the field lines are false-colored by the
values of the force-free factor defined as
$\alpha = \vec J\cdot \vec B/B^2$, which indicates how much the field
lines are non-potential. For a force-free field, this parameter is
constant along any given field line. As can be seen, the magnetic
field is close to a perfect force-free one since the color is nearly
the same along any single field line. In the core of the
configuration, the field lines are sheared significantly along the
PIL, and thus have large values of $\alpha$ and current density. On
the other hand, the overlying field is almost current-free or
quasi-potential field with $\alpha \sim 0$, which plays the role of
strapping field that confines the inner sheared core. Such a
configuration is typical for eruption-productive
ARs~\citep[e.g.,][]{Schrijver2008a, Sun2012, Jiang2016ApJ,
  Toriumi2019, DuanA2019}.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.48\textwidth]{fig4}
  \caption{Evolution of Poynting flux (black line) and the accumulated
    energy (red line) which is time integration of the Poynting
    flux. The gray bars denote the data gaps of the vector
    magnetograms for which a simple linear interpolation in time is
    used to fill.}
  \label{4}
\end{figure}

\begin{figure*}[htbp]
  \centering
  \includegraphics[width=0.8\textwidth]{fig5}
  \caption{Four snapshots of the surface flows (the final smoothed
    version) as derived using the DAVE4VM code.}
  \label{5}
\end{figure*}

\section{Derive the surface flow field}
\label{flow}

Based on the time sequence of vector magnetograms, it is
straightforward to derive the surface velocity by employing the
DARE4VM code developed by~\citet{Schuck2008}. The differential affine
velocity estimator (DAVE) was first developed for estimating
velocities from line-of-sight magnetograms, and was then modified to
directly incorporate horizontal magnetic fields to produce a
differential affine velocity estimator for vector magnetograms
(DAVE4VM). We use the SHARP data with cadence of 12 minutes and pixel
size of 1~arcsec (by rebinning the original data with pixel size of
0.5~arcsec), and first of all, we fill the data gap using a linear
interpolation in the time domain to generate a complete time series of
three days from 00:00~UT on 2014 September 8 to 24:00~UT on September
10 with cadence of 12 minutes. Then we input the time series of vector
magnetogram (rebinned as 1~arcsec per pixel) in the DAVE4VM code. A
key parameter needed by the DAVE4VM code is the window size, and here
we use 19~pixels, following~\citet{Liu2012} and
\citet{LiuC2014}. After obtaining the surface velocity, we first reset
those in the weak-field region (with total magnetic field strength
below 100~G) as zero, because there are large errors and unresolved
scales in these regions. {Figure}~\ref{4} shows evolution of the Poynting
flux ${\rm d}E/{\rm d}t$, which is defined by
\begin{equation}\label{poy}
  \frac{{\rm d} E}{{\rm d} t} = \frac{1}{\mu_0}\int_S \left[(B_{x}^{2}+ B_{y}^{2})v_{z} - ( B_{x}v_{x}+B_{y}v_{y} )B_{z}\right] {\rm d}x{\rm d}y,
\end{equation}
where $S$ is the photospheric surface, and its time accumulation
$E_{\rm acc}$, as computed by using the surface flow $(v_x, v_y, v_z)$
and the magnetic field. It can be seen that, except the data gap
intervals\footnote{During the data gaps the Poynting flux is found to
  become negative adruptly, which might be a result of our simple
  linear interpolation in filling the gaps of the magnetograms. More
  optimized method for filling the data gaps will be considered to
  recover a more consistent evolution of Poynting flux.}, the magnetic
energy is continually injected in the corona through the photosphere,
and in the three days, it gains $\sim 5\times 10^{32}$~erg.

Before being input into the model, the flow data are also needed to be
smoothed. We smoothed the time series of flow maps in both the time
and space domains, with a Gaussian FWHM of 6 hours (i.e., 30 time
snapshots) and 8~arcsec, respectively, which is finally input to the
data-driven model. {Figure}~\ref{5} shows 4 snapshots of the surface
velocity after smoothing. The speed of the flow is generally a few
hundreds of meters per second and the main feature is a clear and
persistent rotation of the main sunspot. Note that during the three
days the basic configuration of the photospheric magnetic flux
distribution is rather similar with only somewhat dispersion. So the
magnetic energy injection should come mainly from the transverse
rotational flows. In addition to the rotational flow, we can see very
evident diverging flow existing persistently near the boundary of the
sunspot.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.45\textwidth]{fig6}
  \caption{Magnetic energy (top) and kinetic energy (bottom) evolution
    of the MHD system as driven by the photospheric flows. The black lines represent the values of the volume-integrated energies. The blue
    line in the top panel shows the time integration of the total
    Poynting flux on the bottom surface from the simulation. The red
    line shows the accumulated magnetic energy from directly the
    observation data, i.e., the same one shown in {Figure}~\ref{4}.}
  \label{6}
\end{figure}

\begin{figure*}[htbp]
  \centering
  \includegraphics[width=0.9\textwidth]{fig7}
  \caption{Evolution of the magnetic field lines (shown in different
    view angles in the left and middle columns) and their comparison
    with observed coronal loops observed in EUV wavelength of 94~{\AA}
    by SDO/AIA (right column). The field lines are shown in the same
    format in {Figure}~\ref{3}. An animation of the magnetic field
    evolution is attached for this figure.}
  \label{7}
\end{figure*}

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.49\textwidth]{fig8}
  \caption{Evolution of magnetic flux distribution on the bottom
    surface in the simulation, shown for 4 snapshots {of the same times shown
    in {Figure}~\ref{5}}.}
  \label{8}
\end{figure}

\section{Data-driven simulation}
\label{evol}

We input the surface velocity at the bottom boundary to drive the
evolution of the model, by starting the simulation from the solution
obtained from the time point of $t=58$ in the relaxation phase (see {Figure}~\ref{2}) as
described in Section~\ref{initial}. Here to save computing time, the
cadence of the input flow maps, which is originally 12 min, was
increased by $68.6$ times when inputting into the MHD model. As a
result, an unit of time in the simulation, $t_s$, corresponds to
actually $t_s \times 68.6=7200$~s, i.e., 2 hours, in the HMI
data. Compressing of the time in HMI data is justified by the fact
that the speed of photospheric flows is often a few
$100$~m~s$^{-1}$. So in our model settings, the evolution speed of the
boundary field, even enhanced by a factor of $68.6$, is still smaller
by two orders of magnitude than the coronal Alfv{\'e}n speed (on the
order of $10^3$~km~s$^{-1}$), and the quick reaction of the coronal
field to the slow bottom changes should not be affected. The
implementation of the bottom boundary conditions is the same as that
for the deeper relaxation phase described in
Section~\ref{initial}. That is, on the bottom surface, we solved the
equation~(\ref{photo_B_equ}) to update all the three components of
magnetic field with the flow field prescribed by those derived in
Section~\ref{flow}, while the plasma density and temperature are
simply fixed. In the driven-evolution phase, the kinetic viscosity is
also used as the smaller one $\nu = 0.05\Delta x^{2}/\Delta t$, which
corresponds to a Reynolds number of $10$ for the length of a grid cell
$\Delta x$.

We show an approximately two-day (a duration of $26~t_s$ or $52$ hours
in reality) evolution of the MHD system as an example, while the
further evolution associated with an eruption and the physical
mechanisms will be left for future study. {Figure}~\ref{6} presents the
global energy evolutions. It can be seen the magnetic energy (the
black line in the top panel) increases monotonously as driven by the
surface flows. By the end of our simulation, it reaches approximately
$1.8$ times of the initial potential energy, and thus the total
magnetic energy obtained in this driving phase is
$(1.8-1.4) \times E_0 \approx 4.8 \times 10^{32}$~erg. On the other
hand, the kinetic energy (bottom panel of {Figure}~\ref{6}) keeps below
the level of $1\times 10^{-3} E_0$ with a mild increase, which
indicates that the system remains a stable, quasi-static evolution. In
the top panel of {Figure}~\ref{6} we also plots the time integration of
total Poynting flux (the blue line), using the magnetic field and the
flow field on the bottom boundary of the simulation, which is the
energy injected into the volume from the bottom boundary through the
surface flow. If our boundary condition is accurately implemented, the
energy injected from the bottom surface should match the magnetic
energy obtained in the computational volume, since other energies are
negligible if compared with the magnetic energy. As can be seen, the
trend of magnetic energy evolution (the black line) matches rather
well that of the energy input by the surface flow (the blue line),
{albeit a small numerical error that increases slowly the total
magnetic energy, which is also seen in the relaxation phase as shown
in {Figure}~\ref{2}}. It is worthy noting that the magnetic energy
evolution is also in good agreement with the accumulated magnetic
energy as derived from directly the observation data (the red line, {which is the same value shown in {Figure}~\ref{4} but multiplied by a factor of $0.05^2$ because the field input to the model is multiplied by $0.05$}) as
calculated in Section~\ref{flow}, which suggests that our data-driven
model can reliably reproduce the magnetic energy injection into the
corona through the photospheric motions. {The small difference between the observation-derived and simulated energies might be due to the smoothing of the velocity since it filters out the small-scale flows that also contribute to the total Poynting flux.}

{Figure}~\ref{7} (and its attached animation) shows the 3D magnetic field
lines and their evolution in comparison with SDO/AIA image of coronal
loops in the 94~{\AA} wavelength which highlights the hot, core
coronal loops in the AR. Overall, we can see a slow stressing of the
field lines mainly through the sunspot rotational flow. It renders the
core field lines to form a more and more coherent S shape, which
resembles the observed sigmoid structure in the core of the AR. The
increase of non-potentiality can also be seen from the increase of the
force-free factor $\alpha$ in the core region. {Figure}~\ref{8} shows
evolution of the surface magnetic flux distribution (see also the
animation attached to {Figure}~\ref{7}). {The four snapshots are taken
for the same times given in {Figure}~\ref{5}, and thus comparing {Figure}~\ref{8} with {Figure}~\ref{5}
shows the difference between the simulated magnetograms and the observed ones.}
In addition to the rotation of
the main sunspot, an evident feature is the enhancement of the field
strength along the PIL. This is owing to the divergence flow from the
edge of the sunspot, which continually convect the magnetic flux to
the PIL. Such pileup of magnetic flux near the PIL, however, is not
seen in the observed magnetograms ({Figure}~\ref{5}), which rather show
field decaying. Such decaying of magnetic flux is likely due to the
global turbulent diffusion of photospheric magnetic field by granular
and supergranular convection~\citep{Wang1989} and other small-scale
turbulence and flux cancellations in the photosphere, which is not
being recovered by the DAVE4VM code and thus not reproduced in our
simulation.

\section{Conclusions}
\label{concl}

This paper is devoted to the description of a new approach of
data-driven modelling of solar AR magnetic field evolution,
{in which, we have for the first time utilized directly the photospheric velocity field from
 DAVE4VM to drive a full-MHD model.}
To setup
the initial conditions, we used a special MHD relaxation approach with
sufficiently small kinetic viscosity to construct a true MHD
equilibrium based on a single vector magnetogram. Then we derived the
photospheric surface flows from a time series of observed magentograms
based on the DAVE4VM method. The surface flows were finally inputted,
again in time sequence, to the bottom boundary of the MHD model to
self-consistently update the magnetic field at every time step, which
is implemented by solving directly the magnetic induction equation at the
bottom boundary  using finite difference method.

We applied this data-driven model to study the magnetic field
evolution of AR~12158 with SDO/HMI vector magnetograms. The initial
MHD equilibrium is calculated using magnetogram observed for 00:00~UT
on 2014 September 8, and a two-day duration of the AR evolution is
then simulated using the data-driven MHD model. Overall, the evolution
is characterized by a slow stress of the field lines mainly through
the rotational flow of the AR's leading sunspot, which makes the core
field lines to forms an coherent S shape consistent with the sigmoid
structure as seen in the AIA images. Such evolution proceeds in a
quasi-static way since the kinetic energy of the system remains less
than its magnetic energy by three orders of magnitude, while the
magnetic energy increases monotonously as driven by the surface flow,
and reaches approximately double of the initial potential energy by
the end of the simulation. The magnetic energy obtained in the
simulation during the surface driving phase matches closely the
accumulated magnetic energy as calculated directly from the original
vector magnetogram with the DAVE4VM derived flow field.

With the surface flow specified at the bottom boundary, the magnetic
field can be updated self-consistently by solving the induction
equation at the surface boundary. However, our simulation shows that
discrepancy arises between the simulated magnetic field at the bottom
surface with the original magnetograms. That is, in the simulation,
magnetic flux are significantly enhanced along the PIL, owing to the
divergence flow from the edge of the sunspot, whereas in the observed
magnetograms such pileup of magnetic flux near the PIL is not seen,
which rather shows field decaying. {In the next step, we will
  try to use some ad-hoc flux cancellation near the PIL, for which a
  straightforward way is to increase the value of $\eta_{\rm stable}$
  in equation~(\ref{photo_B_equ}), such that the flux diffusion speed
  is comparable to that of the flux pileup.}  In the future, we will
consider to include the global diffusion of photospheric magnetic
field by granular and supergranular convection to simulate more
realistically the magnetic field evolution. {On the other hand,
  the discrepancy might result from errors in the vector magnetograms,
  since these errors can introduce spurious flows with the DAVE4VM
  code and thus influences our simulation. To elucidate this, we will
  test our model with error-free magnetograms from recent convective
  flux-emergence
  simulations~\citep[e.g.,][]{ChenF2017,Cheung2019,Toriumi2019S} as
  the ground-truth data.}

{Our ultimate purpose is to use the model} to
study how the coronal field, as driven by the slow photospheric motions,
reaches a unstable state and runs into eruption. From this, we will be able to
investigate in details the topology of the evolving magnetic field leading to the eruption,
to see whether it forms a magnetic flux rope and becomes ideally unstable~\citep[e.g.,][]{Kliem2006, Aulanier2010}, or, a simply sheared
 arcade with a interally-formed current sheet to trigger flare reconnection~\citep{Moore2001}, which
would be helpful for resolving the long-standing debates on the triggering mechanism of
solar eruptions~\citep{ChenP2011}.

\acknowledgments This work is jointly supported by National Natural
Science Foundation of China (NSFC 41822404, 41731067, 41574170,
41531073), the Fundamental Research Funds for the Central Universities
(Grant No.HIT.BRETIV.201901) and Shenzhen Technology Project
JCYJ20190806142609035. Data from observations are courtesy of NASA
SDO. The computational work was carried out on TianHe-1(A), National
Supercomputer Center in Tianjin, China. {We are grateful to the reviewers for
helping to improve the paper.}