Convergence in Hartigan-Wong k-means method and other algorithms
H-W's algorithm, from the 1979 paper, takes as input initial clusters. However, the authors suggest a method for obtaining them in their last section. They write that it is guaranteed that no cluster will be empty after the initial assignment in the subroutine. It goes as follows:


Compute the overall mean $\bar{x}$.
Order the observations with respect to their distance to $\bar{x}$, that is $||x_i - \bar{x}||_2$ (in ascending order I guess?).
Take the points in position $\{ 1 + (L-1) [M/K] \}$, where $L=1, \dots, K$, as initial centroids. ($[\ \cdot\ ]$ most probably refer to the floor function, hence the $1$ at the beginning.)


As for the main algorithm, it is described in a paper called Hartigan's K-Means Versus Lloyd's K-Means-Is It Time for a Change? by N Slonim, E Aharoni, K Crammer, published in 2013 by AJCAI. Note that this version simply uses a random initial partition. It goes as follows.

For the vectors $x \in \mathcal{X}$ and a target number of clusters $K$,


Set $\mathcal{C}$ to be a random partition of $\mathcal{X}$ into $K$ clusters and compute the centroid vector associated to each $C \in \mathcal{C}$, denote them $v_C$.
Scan $\mathcal{X}$ in a random order, and for all $x \in \mathcal{X}$

2.1 Set the stopping indicator to $s = 1$

2.1. Tentatively remove $x$ from its cluster $C$, letting $C^{-} = C \setminus \{ x \}$.

2.2 Find
\begin{align*}
C^+ = \Big\{ \mathrm{argmin}_{C^* \in (\mathcal{C} \setminus C) \cup C^{-}}\ \frac{1}{n} d(x,v_C^*) + \frac{1}{n} \sum_{y \in C^*} [d(y,v_{C^* \cup x}) - d(y,v_{C^*})] \Big\} \cup \{ x \} 
\end{align*}

2.3 If $C^{+} \neq C$, then reset $C \leftarrow C^{-}$ and $C^* \leftarrow C^{+}$ and update their respective vectors $v_{C}$ and $v_{C^*}$. Also reset $s \leftarrow 0$.
If $s=0$, go back to 2.


Here I make a slight abuse of notation letting $C^*$ be the solution of the $\mathrm{argmin}$ in 2.2. This step, 2.2, simply computes the change in the loss obtained if $x$ is added to $C^*$, as opposed to it staying alone in its own cluster (once removed). $d$ stands for distance. Finally, the particular ways to compute the updated vectors $v_C$, $v_{C^* \cup \{x \}}$ and so on is provided in the paper cited above. They can be done in linear time.

I think the answers to all your questions are implicit in the above algorithm... However, I still have to make sure this implementation of the algorithm is standard. In particular if it is the one implemented in R. Any comments/edit are welcomed.