\section*{Abstract}
We address the relationship between membrane microheterogeneity and anomalous subdiffusion in cell membranes by carrying out Monte Carlo simulations of two-component lipid membranes. We find that near-critical fluctuations in the membrane lead to transient subdiffusion, while membrane--cytoskeleton interaction strongly affects phase separation, enhances subdiffusion, and eventually leads to hop diffusion of lipids. Thus, we present a minimum realistic model for membrane rafts showing the features of both microscopic phase separation and subdiffusion.
\newline\newline

\noindent\textbf{PACS:} 87.14.Cc, 87.15.ak, 87.16.dj, 87.16.dt.
\newline

\noindent\textbf{Keywords:} Lipid membranes, DMPC, DSPC, phase separation, membrane rafts, subdiffusion.

\newpage

\section*{Introduction}

Anomalous subdiffusion in cell membranes is an intriguing phenomenon whose molecular origins are still a subject of debate \cite{Saxton, Destainville, Saxton_response}. From the general viewpoint, the phenomenon is related to (dynamic) microheterogeneities in the properties of the cell membrane creating a rugged energy landscape for the diffusing protein or lipid molecule. The current understanding of membrane microdomains embraces the concepts of lipid rafts \cite{Simons} and cytoskeleton-based picket fence \cite{Fujiwara}. The general concept of membrane rafts \cite{Pike} implies a dynamic interplay between the membrane local composition and phase, and local membrane-protein interactions, which can result in local lipid demixing and phase separation. This view is supported by the recent experimental evidence demonstrating that local lipid demixing and phase separation can be induced by crosslinking \cite{Hammond}, change in the local membrane curvature \cite{Roux}, and interaction with a cytoskeleton \cite{Liu}. Additionally, the presence of the membrane-associated actin network leads to a shift in the membrane phase transition temperature and, potentially, to broadening of the phase transition \cite{Liu}.

The efficiency of external perturbations in changing the local properties of the membrane should be strongly enhanced in the vicinity of the membrane critical point. This is indeed supported by recent experimental observations \cite{Sorre}. Moreover, it was suggested that dynamic microdomains in the cell membrane are nothing but near-critical fluctuations in the local composition and phase of the membrane  \cite{Veatch}. All this is in agreement with the observation that the composition of the cell membrane is adjusted to keep it above the phase transition temperature, and is regulated in response to environmental changes to maintain this condition \cite{Mouritsen_chapter}. In cases where this mechanism fails, cold shock damage takes place \cite{Drobnis}.

In this paper, we address the relationship between membrane microheterogeneity and anomalous subdiffusion in cell membranes by carrying out Monte Carlo simulations of two-component lipid membranes on experimentally relevant spatial scales ($\sim$1 $\mu$m) and time intervals ($\sim$1 s), with a special emphasis on interactions with the model membrane skeleton.

We demonstrate that (i) near-critical fluctuations in a free lipid membrane can lead to transient anomalous subdiffusion, (ii) phase separation in two-component (and, hence, multicomponent) lipid membranes can be strongly affected by interaction with the membrane skeleton, which, depending on the temperature and membrane composition, can either lead to precipitation of highly dynamic membrane domains (rafts), or prevent large-scale phase separation, and (iii) interaction with the membrane skeleton enhances anomalous subdiffusion and eventually leads to hop-diffusion of lipids. Thus, we construct a minimum realistic model for membrane rafts showing the features of both microscopic phase separation and anomalous subdiffusion.

The binary lipid system consisting of 1,2-di\-myr\-istoyl-\textit{sn}-glycero-3-phos\-pho\-cho\-line (DMPC) and 1,2-distearoyl-\textit{sn}-glycero-3-phosphocholine (DSPC) is chosen as a generic model of a
two-component membrane. The particular choice of this system is motivated by the fact that it is well studied both in vitro \cite{Shimshick, Mabrey, Lentz, Dijck, Wilkinson, Nibu, Sugar, Hac} and in silico (see \cite{Sugar, Hac, Heimburg_book} and refs. therein). Even though the exact temperature range of phase separation in this membrane is too high for most living organisms, the main qualitative conclusions of the study related to behavior of cell membranes are not affected by the particular thermodynamic parameters of the lipid mixture.

From the physical perspective, it is assumed that living organisms try to maintain the composition of cell membranes to keep them above the phase transition temperature \cite{Mouritsen_chapter}. Therefore, to be biologically relevant, in the present paper we mostly focus on the upper part of the phase diagram.

\section*{Materials and Methods}

\subsubsection*{Experimental}

The saturated phospholipids DMPC and DSPC with the melting temperatures 297 and 328 K were purchased from Avanti Polar Lipids (Alabaster, AL). Multilamellar vesicle (MLV) suspensions were obtained by rehydration of a dry lipid film. Excess heat capacity curves of DMPC/DSPC MLV suspensions in a 10 mM Hepes buffer, pH 7.4, were obtained using a VP-DSC calorimeter (MicroCal, Northampton, MA) at a scan rate of 2--3 K/h.

The empirical phase diagram based on the temperatures of the onset and completion of the phase transition determined from the scanning calorimetry data as described \cite{Sugar, Hac} are in agreement with previously reported results (see Supporting Material).

\subsubsection*{Monte Carlo simulations}

Our approach to lattice-based Monte Carlo (MC) simulations of a two-component membrane is generally similar to the one described previously \cite{Sugar, Hac}. To facilitate efficient simulations on experimentally relevant spatial scales ($\sim$1 $\mu$m) and time intervals ($\sim$1 s), where a particular type of lipid packing and fine molecular details should be of little importance, we further simplified the model and represented the membrane as a square lattice, each node of which represents a molecule of one of the two lipid types, which can be in either gel or fluid conformational state. An elementary MC step consists of two substeps: (i) an attempt to change the state of a randomly chosen lipid and (ii) an attempt of position exchange of a randomly chosen next-neighbor pair of lipids. As in \cite{Hac}, a rate function was introduced for the next-neighbor exchange step to ensure $\sim$40 times slower lipid diffusion in the gel phase. For an $L \times L$ lattice, an MC cycle consists of a chain of $L^2$ elementary MC steps. For every lipid composition and temperature studied, the membrane was first equilibrated, if required, and simulations were run for $(6-20) \times 10^{6}$ MC cycles to collect the necessary data. More details on the simulation procedure are given in the Supporting Material.

Simulations were carried out on an $L \times L = 600 \times 600$ ($400 \times 400$ when modeling the effects of the membrane skeleton) square lattice with periodic boundary conditions. By assuming the average lipid headgroup size of 0.8 nm, this corresponds to a membrane with an experimentally relevant size of $0.48 \times 0.48$ $\mu$m$^2$ ($0.32 \times 0.32$ $\mu$m$^2$). By comparing the DMPC diffusion coefficient obtained in our simulations of pure DMPC in the fluid phase at 304 K with its experimental values at the same temperature ($3 - 6$ $\mu$m$^2$/s \cite{Kuo, Vaz, Dolainsky}), we found that one MC cycle corresponds to $\sim$50 ns, and thus, our simulations cover processes on timescales up to $\sim$0.3--1 s.

\subsubsection*{Model for the membrane skeleton}

The membrane skeleton was modeled by a random Voronoi tessellation satisfying the periodic boundary conditions. For simulations with $L = 400$, random tessellations with $N = 36$ compartments were used, which gives the average linear size of the compartment $\ell = L N^{-1/2} = 66.6$ lattice units $\approx 53$ nm. The generated filament meshwork was projected onto the square lattice thereby creating a set of pixels representing the locations of the filaments. Each of these locations can be assigned to be a cytoskeleton pinning site. To mimic the membrane--cortical skeleton interaction, a simple rule inspired by experimental data on lipid interactions with transmembrane proteins \cite{Kinnunen} was followed: a lipid located at a filament pinning site is forced to assume the gel conformation with no explicit restrictions on its mobility. Thus, the pinning site does not present an obstacle for lipid diffusion: its effect on diffusion is indirect and takes place solely due to a lower lipid mobility in the gel-state local environment. The effect of the varying strength of the membrane-skeleton interaction was modeled by randomly assigning a fraction of filament position pixels to be filament pinning sites. Simulations were carried out with the filament pinning density set to 25\%, 50\%, and 100\%; in these cases the total number of pinning sites amounted to approximately 1\%,  2\%, and 4\% of the total membrane area.

Our choice of immobile pinning sites is justified by experimental observations demonstrating that band 3 and ankyrin strongly bound to the membrane skeleton show a very low diffusion coefficient of $\sim 10^{-4} - 10^{-3}$ $\mu$m$^2$/s over time intervals up to tens of seconds \cite{Tomishige, Cairo}.

It should be pointed out that the approach used in the present work to account for interactions of lipid molecules with membrane proteins is conceptually similar, though not identical, to the one used in several previous MC simulation-based studies \cite{Sperotto, Heimburg1996, Sabra, Gil}. In these works, it was assumed that proteins, preferentially wetted by one of the membrane phases, could freely diffuse in the membrane, which results in their accumulation in this phase and, in case of two-component membranes, in fact enhances large-scale phase separation. (The situation is more complicated in case of active proteins -- for details, see \cite{Gil} and refs. therein.) It was found out that, under certain conditions, interaction of a single-component membrane with small transmembrane proteins or peptides can even lead to emergence of a closed loop of gel--fluid coexistence with a lower critical point \cite{Zhang}. It should be pointed out, however, that these works focused only on structural properties of protein-loaded membranes, and did not address the issues of diffusion in the membrane.

What sets the present work apart from the above-mentioned studies, is that here we consider interaction of the membrane with the \textit{immobile} cytoskeleton and study its effects on phase separation and diffusion of lipids in the membrane.

\subsubsection*{Analysis of lipid diffusion data}

Positions of a small fraction of lipid molecules (150 for simulations with $L = 600$ and 50 for simulations with  $L = 400$, which amounted to 0.04\% and 0.03\%, respectively) were recorded, and the time- and ensemble-averaged mean-square displacement (MSD) was determined as follows:
\begin{equation}
\textit{MSD}(\tau) = \frac{1}{t_\text{max}-\tau} \sum_{t=1}^{t_\text{max}-\tau} \left\{ \left[x(t) - x(t+\tau) \right]^2
													                	              + \left[y(t) - y(t+\tau) \right]^2 \right\},
\label {Eq_MSD}
\end {equation}
where $\tau$, $t$, and $t_\text{max}$ are times measured in units of MC cycles;
here, $\tau$ denotes the time lag, and $t_{\text{max}}$ is the total length of the lipid molecule trajectory. No difference between time- and ensemble-averaged MSD (Eq. \ref{Eq_MSD}) and time-only- and ensemble-only-averaged MSDs beyond the statistical error level was observed.

In case of normal diffusion, the MSD grows in a linear fashion with time: $\textit{MSD}(\tau) = 4 D \tau$ in 2D, where $D$ is the translational diffusion coefficient. In case of anomalous subdiffusion, the MSD shows a slower sublinear power-law growth $\textit{MSD}(\tau) \sim \tau^{\beta}$ with $0<\beta<1$ (see, e.g., \cite{Metzler}). An alternative description of diffusion showing deviations from the normal behavior (also in cases where it cannot be described in terms of subdiffusion) can be provided using an effective time-dependent diffusion coefficient $D(\tau)$.

Therefore, to characterize the behavior of the MSD curves, the local exponent of the mean-square displacement
\begin{equation}
\beta_{\text{MSD}}(\tau) = d\log {\textit{MSD}(\tau)} /d \log \tau,
\label {Eq_beta_MSD}
\end {equation}
and the effective time-dependent diffusion coefficient
\begin{equation}
D(\tau) = \frac{1}{4} d \textit{MSD}(\tau)/d \tau
\label {Eq_D}
\end {equation}
were calculated.

\subsubsection*{Simulation of FCS experiments}

To simulate fluorescence correlation spectroscopy (FCS) \cite{FCS_review} measurements,
the tracked particles (as above, in the amount of 150 for simulations with $L = 600$ and 50 for simulations with  $L = 400$) were assumed to be fluorescent. Fluorescence intensity fluctuations $\delta F(t) = F(t) - \langle F \rangle$ about the mean intensity $\langle F \rangle$ in a 2D Gaussian detection spot $\exp(-2r^2/r_{0}^{2})$ were recorded, and their autocorrelation function $G(\tau) = \langle \delta F(t) \delta F( t+\tau ) \rangle / \langle F \rangle ^{2}$ was calculated. The detection spot size was set to $r_0 = 31$ lattice units $\approx 25$ nm, the size experimentally achievable using the stimulated emission depletion (STED) FCS technique \cite{Kastrup}. We additionally note that this detection spot is much smaller than the lattice size ($r_0/L \approx 0.05$ for $L = 600$ and $r_0/L \approx 0.08$ for $L = 400$), which allowed us to avoid artifacts in the fluorescence autocorrelation function and additionally ensured that our simulations are experimentally relevant.

FCS curves were averaged over nine different positions on the lattice; when studying the effects of the membrane skeleton, the results were additionally averaged over five random realizations of the filament network.

Simulated FCS curves were analyzed using the model
\begin{equation}
G(\tau) = G(0)/\left[ 1+( \tau /\tau _\text{D}) ^{\beta _{\text{FCS}}}\right].
\label {Eq_FCS}
\end {equation}

For $\beta_{\text{FCS}}=1$ this expression corresponds to normal diffusion, while for $0 < \beta_{\text{FCS}} < 1$ it provides a simple way to describe anomalous subdiffusion in FCS \cite{FCS_review}. In case of normal diffusion, $\tau_\text{D}$ is related to the diffusion coefficient of fluorescent particles $D$ and detection spot size $r_0$: $\tau_\text{D} = r_0^2/(4D)$.

Notice that since FCS is sensitive not only to the $\textit{MSD}(\tau)$, but also to the higher moments of the distribution of displacements of a diffusing particle, the connection between $G( \tau )$ and $\textit{MSD}(\tau)$ is straightforward only in case of a Gaussian distribution of displacements \cite{FCS_review}. As a result, in case of anomalous diffusion of particles, generally $\beta _{\text{FCS}} \neq  \beta_{\text{MSD}}$.

\section*{Results and Discussion}

\subsubsection*{Phase and component separation in the membrane}

With appropriate tuning of lipid interaction parameters, the empirical heat capacity-based phase diagram obtained from our MC simulation is in agreement with our experimental data (Fig. \ref{Fig1}), as well as with previously published experimental and simulation data on the same lipid system (see Supporting Material). At the same time, the empirical phase diagram constructed on the basis of heat capacity data may differ from the real phase diagram of the system and thus not provide an insight into the microscopic structure and the dynamics of the membrane.

\begin{figure}[!t]
\begin{center}
\includegraphics{Fig1.pdf}
\caption{Left-hand panel: Component and lipid state separation phase diagram of DMPC/DSPC lipid mixtures. Phase transition temperatures as determined from experimental ($\triangle$) and simulated ($\circ$) excess heat capacity curves. Lipid state spinodal (- - -), lipid state binodal (-- -- --), and lipid demixing curves (---). Right-hand panel: representative snapshots of equilibrium membrane configurations at temperatures and membrane compositions corresponding to filled squares in the left-hand panel. Lattice size: $600 \times 600$; scale bar: 200 lattice units $\approx 160$ nm.}
\label{Fig1}
\end{center}
\end{figure}

More details on the phase diagram can be obtained based on binodal and spinodal curves of the system. The binodal curves, also known as coexistence curves, are the boundaries of the region in the phase diagram in which the equilibrated system shows a complete separation of the two phases. The spinodal encloses the region where the mixture is unstable with respect to small local fluctuations of the composition and always lies inside the area enclosed by the binodal, with the exception of a critical point, where the binodal and spinodal touch.

To reconstruct the binodals and spinodals of the DMPC/DSPC lipid mixture, we analyzed the static structure factors \cite{Fisher} for lipids $S_\text{L}(k)$ and lipid states $S_\text{S}(k)$ (see Supporting Material), reflecting the character of spatial fluctuations of the membrane composition and state.

We found that outside the phase coexistence region, the Ornstein--Zernike (OZ) approximation \cite{Fisher} $S^{\text{OZ}}(k) = S(0)/[ 1 + ( k \xi )^{2}] $, where $k$ is a wavenumber, and $\xi$ is a correlation length, provided an excellent description of the lipid-state structure factors $S_{\text{S}}(k)$. The temperature dependences of $S_{\text{S}}(0)$ and $\xi_{\text{S}}$ were used to estimate spinodal and binodal curves. In particular, the spinodals were determined by extrapolating $1/S_{\text{S}}(0)$ dependence to zero crossing \cite{Fisher, Strobl}, and binodals were determined from the condition $d(1/\xi_{\text{S}})/dT = 0$ \cite{Strobl}. Interestingly, in some regions of the phase diagram, the binodal deviates quite strongly from the phase coexistence boundary estimated from the excess heat capacity data (Fig. \ref{Fig1}). This discrepancy indeed shows that the analysis of heat capacity data does not necessarily recover the real phase diagram of the system. Here, this behavior reflects the continuous character of the phase transition in the membrane, whose proper description requires a combined component and state phase diagram \cite{Michonova-Alexova}.

What process is reflected by the excess heat capacity curves where they fail to describe the fluid--gel phase separation? An analysis of the structure factors $S_{\text{L}}(k)$ for the lipid species helps to answer this question. It appears that, generally, two OZ components are required to describe these data: $S_{\text{L}}(k) = S_{\text{L}1}^{\text{OZ}}(k) + S_{\text{L}2}^{\text{OZ}}(k)$. Parameters of component 1 only weakly depend on the temperature and reflect demixing of lipids in the same state. Parameters of the second component $S_{\text{L}2}(0)$ and $\xi_{\text{L}2}$ show strong temperature dependences similar to those of $S_{\text{S}}(0)$ and $\xi_{\text{S}}$, and thus describe appearance of dynamic microscopic domains, which cannot be treated as distinct thermodynamically stable phases \cite{Heimburg_book}. We therefore define the lipid demixing curves as a set of points on the phase diagram satisfying the condition $S_{\text{L}1}(0) = S_{\text{L}2}(0)$. The fact that the lipid demixing curves are in excellent agreement with the calorimetry-based empirical phase diagrams, supports the above reasoning, as do the snapshots presented in Fig. \ref{Fig1}. Notice that in this region the membrane undergoes near-critical fluctuations. On the other hand, in the coexistence region, equilibrated membranes show complete phase separation and form large-scale lipid domains (Fig. \ref{Fig1}).

The above analysis of the structure factors $S_{\text{L}}(k)$ and $S_{\text{S}}(k)$ and equilibrated membrane snapshots shows that in the part of the phase diagram where the lipid demixing curve closely approximates the binodal, the phase transition from the fluid state to the fluid--gel coexistence has a quasi-abrupt character. On the other hand, in the region where the lipid demixing curve strongly deviates from the binodal, the membrane shows characteristic near-critical fluctuations (Fig. \ref{Fig1}), and the transition becomes continuous as the system approaches the critical point (DMPC/DSPC $\approx$ 20:80, $T_{\text{c}} = 320.5 \pm 0.2$ K -- see \cite{NJP} for details). Remarkably, this behavior is qualitatively similar to recent experimental observations on three-component lipid membranes \cite{Veatch}: depending on the membrane composition, the transition to the two-phase coexistence takes place either in an abrupt manner -- when the membrane does not pass through a critical point, or via critical fluctuations -- when the membrane does pass through a critical point. A more detailed consideration of the phase diagram of the DMPC/DSPC system based on MC simulations will be published elsewhere \cite{NJP}.

In the fluid--gel phase coexistence region, large-scale phase separation takes place. Since no explicit or implicit penalties are imposed \cite{Baumgart2003, Semrau, Ursell} on the domain size in our model, and domain growth is driven by minimization of the line tension energy, phase separation results in formation of a single circular-shaped domain of the minority phase (Fig. \ref{Fig1}).

\subsubsection*{Effect of the cytoskeleton on phase separation in the membrane}

In the region of near-critical fluctuations the membrane is expected to be very sensitive to external perturbations, including the interaction with the membrane skeleton. It appears that in the vicinity of the critical point the interaction with the membrane skeleton leads to immediate condensation of the gel phase on the skeleton filaments, and formation of membrane domains, in a good agreement with results of experiments on lipid bilayers with a reconstituted actin skeleton \cite{Liu}. These domains dynamically change their shape, but nevertheless stay pinned to the filaments (Fig. \ref{Fig2}). Notice that the effect of the filaments is remarkably robust with respect to the filament pinning density. We strongly believe that these domains represent the minimal model of membrane rafts \cite{Pike}. The minimum character of this model stems from the fact that in the present scenario domain formation is thermodynamically driven, and does not require any active processes like chemical cross-linking of membrane components or lipid recycling.

\begin{figure}[!t]
\begin{center}
\includegraphics{Fig2.pdf}
\caption{Effect of the membrane skeleton on the phase separation in a DMPC/DSPC membrane. Representative snapshots of membrane configurations are shown for the free membrane (\textit{first row}) and membrane interacting with a network of filaments at 100$\%$ (\textit{second row}), 50$\%$ (\textit{third row}) and 25$\%$ pinning density (\textit{fourth row}). Snapshots for the free membrane, as well as for the membrane interacting with filaments at $T = 321$ and 322 K, represent fully equilibrated configurations; snapshots for the membrane interacting with filaments at $T = 310$ and 317.5 K correspond to equilibration time of $6 \times 10^6$ MC cycles (see text for discussion). For presentation purposes, the filaments are drawn thicker than they are in reality. Lattice size: $400 \times 400$; scale bar: 125 lattice units $\approx 100$ nm.}
\label{Fig2}
\end{center}
\end{figure}

If the membrane is abruptly cooled down from the all-fluid state to a temperature in the fluid--gel coexistence region of the phase diagram, phase separation takes place, and domains of the fluid and gel phase start to nucleate and coarsen with time. In a free membrane, domains grow according to the power law $R(t) \sim t^n$ with the growth exponent $n$  depending on the particular domain growth mechanism \cite{Furukawa}. In our simulations for a free membrane, depending on the lipid composition and temperature, $n$ takes values from $1/4$ to $1/3$ (Fig. \ref{Fig3}), consistent with general expectations for the domain growth in 2D \cite{Furukawa} and in agreement with experimental results (see, e.g., \cite{Jensen}) (for a more detailed discussion, see \cite{NJP}).

In a finite-size system with a linear dimension $L$, domain coarsening eventually stops, and a single circular-shaped domain is produced (Fig. \ref{Fig1}) with the radius $R_\infty \simeq (X/\pi)^{1/2}L$, where $0 < X < 1/2$ is the fraction of the minority phase.

Remarkably, we found that the presence of the membrane skeleton strongly inhibits or even eventually prevents large-scale phase separation in the phase coexistence region (Fig. \ref{Fig2}). As a result, the radius of the membrane domains in this case is largely determined by the characteristic compartment radius $R_{\text{comp}} \simeq \ell/2$ of the filament network. The way we account for interaction of lipid molecules with filaments at their pinning sites is similar to a theoretical model (the random-field Ising model \cite{RFIM}) implying the presence of static random position-dependent perturbations in a 2D system of spins. This model predicts that the initial power-law domain growth $R(t) \sim t^n$ is strongly slowed down at intermediate stages, and crosses over to an extremely slow logarithmic growth $R(t) \sim \log t$; eventually the domains are expected to reach the perturbation strength-dependent maximum size $R' < R_\infty$ in an exponential time to create an equilibrium disordered state \cite{RFIM_theory, RFIM_simulations}. In our system, the perturbation strength is determined by the filament pinning density and the average compartment radius $R_{\text{comp}}$. Therefore, if the interaction of the membrane with the cytoskeleton is strong enough, and ${R_\infty}^2 \gg R_{\text{comp}}^2$, one can expect that the domain growth stops when domains reach the characteristic size $R_{\text{comp}} \lesssim R' < R_\infty$.

\begin{figure}[!t]
\begin{center}
\includegraphics{Fig3.pdf}
\caption{Effect of the membrane skeleton on the the domain growth in DMPC/DSPC 50:50 membrane abruptly cooled from the all-fluid state down to $T = 310$ K in the fluid--gel phase coexistence region. Kinetics of domain growth when the membrane is free (\textit{upper curve}) and in the presence of membrane skeleton with the filaments pinning density of 50$\%$ (\textit{lower curve}). Black solid line shows the power law dependence $R(t) \sim t^n$ with $n = 0.32$; dashed line shows the stage of the slow logarithmic growth $R(t) \sim \ln t$ in the presence of membrane skeleton. Representative membrane configurations obtained in our MC simulations at time instants $10^4$, $10^5$, and $10^6$ MC cycles are shown at the corresponding curves. Lattice size: $400 \times 400$.}
\label{Fig3}
\end{center}
\end{figure}

We indeed observed an extreme slowing down of the domain growth in the presence of the membrane skeleton and cross-over to the slow logarithmic growth (Fig. \ref{Fig3}). As is evident from the Figure, at the end of a simulation run of $6 \times 10^6$ MC cycles, the domain sizes are indeed $R(t) \sim R_{\text{comp}}$. Although some growth is still observed at this stage, it is so slow, that at present it is unclear, whether the system will evolve toward a disordered equilibrium state featuring a number of small domains, as suggested in \cite{RFIM_theory, RFIM_simulations}, or the slow logarithmic growth will continue further to produce eventually (after an extremely long time) a single domain with the size $R_\infty$. What is clear, though, that the time required for complete equilibration is so long, that from the practical viewpoint one can state that indeed the presence of the membrane skeleton prevents large-scale phase separation in the membrane. This becomes especially clear if one takes into account that a membrane in a live cell is not in the equilibrium state, and a number of other processes affecting the membrane state and composition, e.g., lipid recycling, take place in parallel on different spatial and time scales.

Thus, the cytoskeleton-induced inhibition of large-scale phase separation can serve as one of the possible explanations why micrometer-scale membrane domains are observed in giant plasma membrane vesicles \cite{Plasma_membrane_vesicles}, but not in cell membranes.

In addition, our observations suggest that the established cryoprotective role played by the cytoskeleton \cite{Clark} may consist in delaying phase separation-induced cold shock damage of living cells.

\subsubsection*{Diffusion of lipids in the free membrane}

At higher temperatures, away from the phase transition and the near-critical fluctuations region, the membrane is in the homogeneous all-fluid state, and, not unexpectedly, lipid diffusion is normal (Fig. \ref{Fig4}).

This picture changes upon approaching the critical point in the region of near-critical fluctuations where  interpenetrating fluctuating domains form in the membrane. The analysis of the immediate environment of the DMPC and DSPC lipids shows that, while the DSPC lipid does not have a pronounced preference for the state of its local environment and therefore is evenly distributed between fluid and gel domains, the DMPC lipid shows a strong preference for the fluid local environment and is thus predominantly partitioned into fluid domains. Under these conditions, the DSPC lipid shows no significant deviation from the normal diffusion behavior (data not shown). By contrast, quite unexpectedly, the DMPC lipid was found to demonstrate very pronounced anomalous subdiffusion (Fig. \ref{Fig4}).

\begin{figure}[!t]
\begin{center}
\includegraphics{Fig4.pdf}
\caption{Effect of the proximity to the phase transition on diffusion of DMPC lipids in a DMPC/DSPC 20:80 membrane. Mean-square displacement $\textit{MSD}(\tau)$ (a), local exponent $\beta_{\text{MSD}}(\tau)$ (b), time-dependent diffusion coefficient $D(\tau)$ (c), and FCS autocorrelation $G(\tau)/G(0)$ (d). Panels (a) -- (c) top to bottom: $T = 328$, 322, and 321 K. Panel (d): $T = 328$ (\textit{left}) and 321 K (\textit{right}). For clarity, data for 322 K are omitted in panel (d). Panel (d) additionally shows fits to the FCS diffusion model Eq. (\ref{Eq_FCS}) giving $\beta_{\text{FCS}} = 1.01$ at 328 K and $\beta_{\text{FCS}} = 0.86$ at 321 K (-- -- --). For comparison, a fit of 321 K data with fixed $\beta_{\text{FCS}} = 1.0$ is shown (- - -).}
\label{Fig4}
\end{center}
\end{figure}

At a first glance, this subdiffusive behavior is rather surprising, especially in the light of the seminal work by Kawasaki \cite{Kawasaki} demonstrating that the diffusion coefficient of a two-component mixture does not vanish upon approaching the critical point, which implies normal diffusion, at least in the long-time asymptotic regime. A closer look at the data shows, however, that the subdiffusive behavior is in fact transient, though it covers several orders of magnitude in time, and at longer times the cross-over to the normal diffusion takes place. One can also notice that at very short time lags, diffusion is normal as well. This behavior is also evident from the plots of the local exponent $\beta_{\text{MSD}}(\tau)$ and effective time-dependent diffusion coefficient $D(\tau)$ extracted from the $\textit{MSD}(\tau)$ dependences (Fig. \ref{Fig4}b, c).

How this cross-over from normal diffusion to a subdiffusive behavior and back to normal diffusion at long times can be explained? By recalling that the DMPC lipid is preferentially partitioned into fluid domains and having a look at the corresponding membrane snapshots (Fig. \ref{Fig1}), we realize that DMPC molecules diffuse on a dynamically rearranging fractal-like fluid phase pattern. Therefore, at early times, when the DMPC molecule explores its immediate environment, normal diffusion should take place. Later, though at times shorter than a characteristic time of rearrangement of this dynamic fractal structure, the motion of DMPC lipid molecules is subdiffusive, similar to what is expected in case of diffusion on fractal-like percolation clusters \cite{ben-Avraham_Havlin}. At much longer times, evolution of these clusters leads to dynamic percolation behavior (see, e.g., \cite{dynamic_percolation}), which results in a crossover from the subdiffusive motion to normal diffusion. In contrast to diffusion on static percolation clusters, where cross-over to normal diffusion occurs only above the percolation transition \cite{ben-Avraham_Havlin}, in the case of lipid diffusion in a near-critical membrane, cross-over to normal diffusion will always take place irrespectively of whether the fluid phase is above or below the percolation threshold.

This subdiffusive behavior is also clearly observed in our simulations of FCS experiments (Fig. \ref{Fig4}d). Upon a closer approach to the critical point of the system, the autocorrelation functions of fluorescence intensity fluctuations $G(\tau)$ progressively stronger deviate from the dependence expected for normal diffusion (i.e., Eq. (\ref{Eq_FCS}) with $\beta_{\text{FCS}} = 1$). A good description of $G(\tau)$ data in this case can only be obtained with $\beta_{\text{FCS}} < 1$. For example, we find that for DMPC/DSPC 20:80 mixture at $T = 321.0$ K the best fit is obtained with $\beta_{\text{FCS}} = 0.86$ (Fig. \ref{Fig4}d), and a smaller $\beta_{\text{FCS}} = 0.79$ is required to fit $G(\tau)$ at $T = 320.7$ K, which is closer to the critical point (data not shown).

A new important result of the present work is that transient anomalous subdiffusion spanning several orders of magnitude in time can be observed in a system close to its critical point. We emphasize that the appearance of the transient subdiffusive behavior in the region of near-critical fluctuations is not related to any specific properties of the system and should be observed close to criticality in various systems independent of their origin, including, of course, multicomponent lipid membranes.

It is well known that phase separation in lipid bilayers produces domains with typical sizes ranging from a few to several tens of micrometers (provided that there are no restrictions on the domain growth). In agreement with that, we observe in our simulations that within the coexistence region an equilibrated membrane always shows complete phase separation resulting in a single circular-shaped domain of the minority phase. Therefore, in this case, the only reasonable way to carry out FCS measurements requires parking the detection spot into the bulk of the majority phase away from the inter-phase boundary, or into the center of the minority phase domain. Exactly this approach is used in experimental FCS studies on membranes showing large-scale phase separation (see, e.g., \cite{Bacia, Kahya, Chiantia}). Not unexpectedly, this approach results in normal diffusion in both phases with phase-dependent diffusion coefficients. This behavior is also observed in our simulations (data not shown).

In contrast to that and quite surprisingly, an MC simulation study of Hac et al. \cite{Hac} reported FCS curves strongly deviating from the normal diffusion model exactly in the phase coexistence region, i.e., where large-scale phase separation takes place. These results were attributed by the authors of \cite{Hac} as resulting from subdiffusive motion of lipids, which is in contradiction with the macroscopic phase separation in the system. This seeming contradiction is resolved upon careful examination of the approach to simulations in \cite{Hac}. There, the diameter of the FCS detection spot was approximately equal to the simulation box size. In this case, the FCS detection spot covers effectively the whole simulated membrane, and, even in the presence of large-scale phase separation, motion of lipids in both fluid and gel phases contributes to the FCS results, which explains deviations of FCS curves from the normal diffusion model observed in \cite{Hac} in the fluid--gel phase coexistence region and rule out the interpretation of these results in terms of subdiffusion. Moreover, the situation simulated in \cite{Hac} is very unlikely for experiments on single lipid bilayers within the phase coexistence region, where membrane domains typically have radii of a few to several tens of micrometers, whereas the typical FCS detection spot size is $\sim 200$ nm for the standard FCS \cite{FCS_review}, and down to 25 nm for STED FCS \cite{Kastrup}.

Here we should also mention a previous lattice-based simulation study \cite{Sugar2005main} addressing the character of lateral diffusion of lipid molecules in DMPC/DSPC 50:50 membranes. There, in agreement with our results, diffusion is normal at $T > 320$ K, i.e. in the all-fluid membrane state (cf. Fig. \ref{Fig1}). Deviations from the normal diffusion for the DMPC/DSPC 50:50 mixture found in \cite{Sugar2005main} within the temperature range $T = 300 - 320$ K (i.e. exactly where the large-scale phase separation for this lipid composition takes place -- see Fig. \ref{Fig1}) were interpreted in terms of subdiffusive motion of lipids. As in \cite{Hac}, the spatial scales on which deviations from normal diffusion were observed in \cite{Sugar2005main} approach the system size. Hopping of lipid molecules between the macroscopic gel and fluid domains definitely leads to deviations from normal diffusion at these large observation scales, although these deviations can hardly be interpreted in terms of anomalous subdiffusion.

We emphasize that our results presented in this paper unambiguously show that local component and phase fluctuations in a near-critical membrane lead to the transient anomalous subdiffusion spanning several orders of magnitude in time. Under suitable experimental conditions, we believe, this subdiffusive behavior can be observed experimentally using, e.g., the (STED) FCS technique.

\subsubsection*{Effects of membrane--cytoskeleton interaction on lipid diffusion}

Interaction with the membrane skeleton strongly affects the character of diffusion in the region of near-critical fluctuations (Fig. \ref{Fig5}). With increasing the filament pinning density, the interaction of the membrane with filaments becomes more pronounced: it first enhances the anomalous diffusion and eventually leads to hop-diffusion of lipid molecules. In the presence of filaments, the faster diffusion process clearly corresponds to Brownian motion within compartments defined by the membrane skeleton, whereas the slower diffusion is due to hopping between the compartments, and therefore it strongly depends on the filament pinning density. Notice that gel phase condensation at the membrane skeleton substantially increases the effective thickness of filaments and thus enhances their influence on lipid diffusion.

\begin{figure}[!t]
\begin{center}
\includegraphics{Fig5.pdf}
\caption{Effect of the membrane skeleton on diffusion of DMPC lipids in a DMPC/DSPC 20:80 membrane at 322 K. Mean-square displacement $\textit{MSD}(\tau)$ (a), local exponent $\beta_{\text{MSD}}(\tau)$ (b), time-dependent diffusion coefficient $D(\tau)$ (c), and FCS autocorrelation $G(\tau)/G(0)$ (d). Panels (a) -- (c) top to bottom, and panel (d) left to right: filament pinning density equals 0$\%$ (free membrane), 25$\%$, 50$\%$, and 100$\%$.}
\label{Fig5}
\end{center}
\end{figure}

Remarkably, in this case our results for the mean-square displacement $\textit{MSD}(\tau)$ (Fig. \ref{Fig5}a) and effective time-dependent diffusion coefficient $D(\tau)$ (Fig. \ref{Fig5}c) are in good qualitative agreement with results of single-molecule tracking experiments in cell membranes \cite{Murase}, where results were interpreted in terms of a picket-fence model. Our $\textit{MSD}(\tau)$ data also qualitatively agree with a theoretical model for diffusion on an infinite periodic square meshwork of penetrable barriers \cite{Wieser}, except for in our case the short-time motion of a molecule within a compartment is subdiffusive due to near-critical fluctuations. We again emphasize that even a relatively low pinning density of the cytoskeleton filaments leads to appearance of rather strong barriers for diffusing lipids due to filament-induced local condensation of the gel phase.

Our simulations of STED FCS experiments with the experimentally achievable detection spot size \cite{Kastrup, Eggeling} show (Fig. \ref{Fig5}d) that, while at early times all FCS curves exhibit behavior characteristic of subdiffusion, interaction with the membrane skeleton leads to appearance of a second slow diffusion component. The contribution of this slow component increases with the filament pinning density. Clearly, this component is related to the large-scale diffusion that involves crossing the filament-induced barriers on the membrane. This two-component character of $G(\tau)$ is in qualitative agreement with the recent STED FCS measurements of lipid diffusion in cell membranes \cite{Eggeling}. It would be very interesting to follow the dependence of the FCS autocorrelation functions and the respective apparent diffusion coefficients on the spatial scale using the variable detection spot FCS technique \cite{Masuda, Wawreziniek}, as has been done in \cite{Eggeling}. This work is currently in progress.

Here, we point out once again that in our model the pinning sites do not represent obstacles for diffusing lipids, and affect lipid diffusion only indirectly, via the lower lipid mobility in the local gel-state environment (see Materials and Methods). We anticipate that in case where pinning sites would be presented by immobile particles preferentially wetted by gel-phase lipids, the effects of the membrane--filament interaction on diffusion of lipids should become even more pronounced. This issue will be addressed in our future work.

\section*{Conclusion}
In conclusion, in this paper, we carried out Monte Carlo simulations of model two-component lipid membranes on experimentally relevant spatial scales and time intervals. This allowed us to better understand the details of the dynamics of phase and component separation, and by this means to address the relationship between membrane microheterogeneity and anomalous subdiffusion in cell membranes.

We observed that interaction of the near-critical membrane with the cytoskeleton strongly affects phase separation: for the membrane in the state of near-critical fluctuations it leads to local phase separation and formation of dynamic domains with a size of a few tens of nanometers, which conform to the current understanding of membrane rafts \cite{Pike}.

We find that in a membrane showing near-critical fluctuations lipids show subdiffusive behavior covering several orders of magnitude in time. The interaction of the membrane with the cytoskeleton enhances subdiffusion and eventually leads to hop-diffusion of lipids.

In the fluid--gel phase coexistence region of the phase diagram, interaction of the membrane with the cytoskeleton is found to tremendously slow down large scale phase separation in the membrane, which may explain the established connection between the cytoskeleton and cold shock resistance of organisms \cite{Clark}.

The concepts of lipid rafts \cite{Simons} and cytoskeleton-related picket fence \cite{Fujiwara} are frequently discussed as the alternative viewpoints on the origin of anomalous diffusion in cell membranes. In our paper, we bring these two concepts together to show that not only they do not contradict one another, but rather work in synergy, resulting in formation of cytoskeleton-induced dynamic lipid domains in the near-critical membrane. By this means we construct what, we believe, is a minimum raft model, since the domain formation is driven solely by thermodynamics and does not require either chemical cross-linking of membrane components, or lipid recycling.

\section*{Acknowledgements}
We acknowledge inspiring discussions with H. Rigneault and C. Favard at the early stage of the project.

The work was supported by the Deutsche Forschungsgemeinschaft via Research Group FOR 877.

\newpage


\section{Introduction} \label{sec_intro}
\setcounter{equation}{0}

In this paper we study numerical approximations of the electric field integral equation (EFIE)
on a surface $\Gamma = \partial\Omega$, where $\Omega \subset {\field{R}}^3$ is
a Lipschitz polyhedron. The EFIE is a boundary integral formulation of a
boundary value problem for the time-harmonic Maxwell equations
in the domain exterior to $\Omega$.
It models the scattering of electromagnetic waves at a perfectly conducting body.

For the numerical solution of the EFIE we use
the Galerkin boundary element method (BEM).
It employs ${\mathbf H}({\mathrm div}_{\G})$-conforming families of surface elements,
namely Raviart-Thomas (RT) and Brezzi-Douglas-Marini (BDM) elements, to discretise
the variational formulation of the EFIE (called Rumsey's principle).
This approach is referred to as the natural BEM for the EFIE.
As in finite element methods, the convergence may be achieved either
by keeping polynomial degrees fixed and refining the mesh ($h$-version),
or by fixing the mesh and increasing polynomial degrees ($p$-version),
or by simultaneous $h$-refinement and $p$-enrichment ($hp$-version).

The Galerkin BEM is widely used in the engineering practice for the
simulation of electromagnetic scattering. Moreover, it had been used
long before a rigorous theoretical analysis of the method became available.
Error analyses of different boundary element schemes for the eddy current problem on smooth obstacles are given by MacCamy and Stephan in \cite{McCamyS_83_BET,McCamyS_84_SPT}.
Despite these early results the BEM-analysis for non-smooth obstacles started much later.
In fact, the convergence and a priori error analysis of the $h$-BEM for the EFIE
on piecewise smooth (open or closed) surfaces
has been developed within the last decade
(see \cite{BuffaCS_02_BEM,HiptmairS_02_NBE,BuffaHvPS_03_BEM,BuffaC_03_EFI,BuffaH_03_GBE}),
and the corresponding results for high-order methods ($p$- and $hp$-BEM)
are very recent (see \cite{BespalovH_NpB,BespalovH_10_hpA,BespalovHH_Chp}).
In particular, \cite{BespalovH_NpB} and \cite{BespalovHH_Chp} establish quasi-optimal convergence
of the natural $hp$-BEM on meshes of shape-regular elements.
An essential ingredient for the proofs in these papers are the projection-based
interpolation operators developed by Demkowicz and co-authors
\cite{DemkowiczB_03_pIE,DemkowiczB_05_HHH,Demkowicz_08_PES}.
These operators also facilitate a priori error analysis, where one
needs $hp$-approximation theory in specific trace spaces.
In \cite{BespalovH_10_hpA} we prove an a priori error estimate
of the $hp$-BEM on quasi-uniform meshes of affine elements under the assumption
that the regularity of the solution to the EFIE is given in Sobolev
spaces on $\Gamma$. The latter result states that the
method converges with the same rate $r + \frac 12$ in both $h$ and $p^{-1}$ (here,
$r$ denotes the Sobolev regularity order, and $p$ is assumed to be large enough).
However, it is evident from the numerical results reported in \cite{Leydecker_Thesis} that
the $p$-BEM converges faster than the $h$-BEM for the EFIE.
This is similar to what was observed and proved for the $p$-BEM and the $h$-BEM
applied to elliptic problems in three-dimensions (cf.
\cite{SchwabS_96_OpA,HeuerMS_99_ECB,BespalovH_08_hpB,BespalovH_10_hpW}).

With this paper we fill a gap in the theory of the BEM for the EFIE on polyhedral
surfaces by proving a precise error estimate for the $hp$-BEM on quasi-uniform meshes.
Similarly to the elliptic case we make use of explicit expressions for singularities in
electromagnetics fields (these expressions are available from
\cite{CostabelD_00_SEF,BespalovH_NpB}). The established error
estimate shows that the convergence
rates in $h$ and $p$ depend on the strongest singularity exponent, and that the $p$-BEM
converges twice as fast as the $h$-BEM on quasi-uniform meshes.
This extends the results of \cite{BespalovH_NpB}, where the
analysis was restricted to the $p$-BEM on a plane open surface.

It is now well known that appropriate decompositions of vector fields
(on both the continuous and discrete level) are
critical for the analysis of electromagnetic problems and their approximations.
In the case of the EFIE the main idea is to isolate the kernel of the ${\mathrm div}_{\G}$-operator
such that the complementary field possesses an enhanced smoothness
(cf. \cite{BuffaC_03_EFI,Buffa_05_RDS,BespalovHH_Chp}).
The corresponding decompositions of singular vector fields greatly facilitate
our analysis in this paper as well. In particular, the complementary vector fields are singular
vector functions belonging to ${\mathbf H}^{1/2}(\Gamma)$. Then, in the case of the BDM-based BEM,
these vector functions can be approximated component-wise by continuous piecewise polynomials,
and the desired $hp$-error estimates are derived by using the corresponding
results in \cite{BespalovH_08_hpB} for scalar singularities (belonging to $H^{1/2}(\Gamma)$).
However, this simple approach does not work for the RT-based BEM because
the dimension of the underlying RT-space on the reference element is smaller
than the dimension of the BDM-counterpart. That is why, the results
of \cite{BespalovH_08_hpB} cannot be applied directly, and additional technical
arguments are needed (see the proof of Lemma~\ref{lm_approx_v}
in Section~\ref{sec_proof_v}).

The rest of the paper is organised as follows.
In the next section we formulate the EFIE in its variational form
and recall the typical structure of the solution to this model problem.
In Section~\ref{sec_BEM} we introduce the $hp$-version of the BEM
for the EFIE and formulate the main result of the paper (Theorem~\ref{thm_main_hp})
stating convergence rates of the method. This result follows from
the general approximation theorem (Theorem~\ref{thm_gen_approx}) established
in Section~\ref{sec_gen_approx}. The proof of Theorem~\ref{thm_gen_approx}
relies, in particular, on two technical lemmas which are proved in
Section~\ref{sec_proofs}.

Throughout the paper, $C$ denotes a generic positive constant
independent of $h$ and $p$.

\section{Formulation of the problem} \label{sec_problem}
\setcounter{equation}{0}

Let $\Gamma$ be a Lipschitz polyhedral surface in ${\field{R}}^3$.
Throughout the paper we will use exactly the same
notation as in \cite{BespalovH_10_hpA} for all involved differential and
boundary integral operators as well as for Sobolev spaces of scalar functions
and tangential vector fields on $\Gamma$ (all essential definitions are given
in \cite[Section~3.1]{BespalovH_10_hpA}). In particular, we use boldface symbols
for vector fields, and the spaces (or sets) of vector fields are denoted in
boldface as well.

For a given wave number $\kappa > 0$, we denote by $\Psi_\kappa$ (resp., $\hbox{\mathversion{bold}$\Psi$}_\kappa$)
the scalar (resp., the vectorial) single layer boundary integral operator
on $\Gamma$ for the Helmholtz operator $-\,\Delta - \kappa^2$
(see \cite[Section~5]{BuffaH_03_GBE}).
The variational formulation for the EFIE will be posed in the following
Hilbert space of tangential vector fields on $\Gamma$:
\[
  {\mathbf X} = {\mathbf H}^{-1/2}({\mathrm div}_{\G},\Gamma) :=
        \{{\mathbf u} \in {\mathbf H}^{-1/2}_{\|}(\Gamma);\; {\mathrm div}_{\G}\,{\mathbf u} \in H^{-1/2}(\Gamma)\},
\]
which is the trace space of ${\mathbf H}({\mathbf curl},\Omega)$, where $\Omega \subset {\field{R}}^3$
is a Lipschitz polyhedron such that $\Gamma = \partial\Omega$
(we refer to \cite{BuffaC_01_TFI,BuffaC_01_TII,BuffaCS_02_THL}
for the definition and properties of this and other trace spaces on $\Gamma$).

Let ${\mathbf X}'$ be the dual space of ${\mathbf X}$ (with duality pairing
extending the ${\mathbf L}^2(\Gamma)$-inner product for tangential vector fields).
Then, for a given source functional ${\mathbf f}\in {\mathbf X}'$
the variational formulation of the EFIE reads as:
{\em find a complex tangential field ${\mathbf u}\in{\mathbf X}$ such that}
\begin{equation} \label{bie_var}
    a({\mathbf u},{\mathbf v})
    :=
    \<\Psi_\kappa {\mathrm div}_{\G}\,{\mathbf u}, {\mathrm div}_{\G}\,{\mathbf v}\> - \kappa^2 \<\hbox{\mathversion{bold}$\Psi$}_\kappa{\mathbf u}, {\mathbf v}\>
    = \<{\mathbf f}, {\mathbf v}\> \quad\forall {\mathbf v} \in {\mathbf X}.
\end{equation}
Here, the brackets $\<\cdot,\cdot\>$ denote
dualities associated with $H^{1/2}(\Gamma)$ and ${\mathbf H}^{1/2}_{\|}(\Gamma)$.
To ensure the uniqueness of the solution to (\ref{bie_var}) we always assume
that $\kappa^2$ is not an electrical eigenvalue of the interior problem.

Let us recall the typical structure of the solution ${\mathbf u}$ to problem (\ref{bie_var}),
provided that the source functional ${\mathbf f}$ is sufficiently smooth
(we note that this regularity assumption is satisfied for the electromagnetic
scattering with plane incident wave). We use the results of
\cite[Appendix~A]{BespalovH_NpB}. These results were derived from the regularity
theory in \cite{CostabelD_00_SEF} for the boundary value problems
for Maxwell's equations in 3D by making use of trace arguments and some
technical calculations in the spirit of \cite{vonPetersdorffS_90_DEC} and \cite{vP}.

Let $V = \{v\}$ and $E = \{e\}$ denote the sets of vertices and edges of $\Gamma$, respectively.
For $v\in V$, let $E(v)$ denote the set of edges with $v$ as an end point.
Then the solution ${\mathbf u}$ of (\ref{bie_var}) can be written as
\begin{equation} \label{dec}
    {\mathbf u} = {\mathbf u}_{\rm reg} + {\mathbf u}_{\rm sing},
\end{equation}
where
\begin{equation} \label{reg}
    {\mathbf u}_{\rm reg} \in
    {\mathbf X}^k :=
    \{{\mathbf u} \in {\mathbf H}^k_{\;-}(\Gamma);\; {\mathrm div}_{\G}\,{\mathbf u} \in H^k_{\;-}(\Gamma)\}\ \ \hbox{with \ $k>0$}
\end{equation}
(the spaces ${\mathbf H}^k_{\;-}(\Gamma)$ and $H^k_{\;-}(\Gamma)$ are
defined in a piecewise fashion by localisation to each face of $\Gamma$,
and the space ${\mathbf X}^k$ is equipped with its graph norm $\|\cdot\|_{{\mathbf X}^k}$,
see \cite{BespalovH_10_hpA}),
\begin{equation} \label{sing}
    {\mathbf u}_{\rm sing} = \sum_{e\in E} {\mathbf u}^e + \sum_{v\in V} {\mathbf u}^v
                     + \sum_{v\in V}\sum_{e\in E(v)} {\mathbf u}^{ev},
\end{equation}
and ${\mathbf u}^e$, ${\mathbf u}^v$, and ${\mathbf u}^{ev}$ are the edge, vertex, and
edge-vertex singularities, respectively.

In order to write explicit expressions for the above singularities,
let us fix a vertex $v \in V$ and an edge $e \in E(v)$.
Then, on each face $\Gamma^{ev} \subset \Gamma$ such that $e \subset \partial\Gamma^{ev}$
we will use local polar and Cartesian coordinate systems $(r_v,\theta_v)$
and $(x_{e1},x_{e2})$, both with the origin at $v$, such that
$e = \{(x_{e1},x_{e2});\; x_{e2} = 0,\ x_{e1} > 0\}$ and
for a sufficiently small neighbourhood $B_\tau$ of $v$ there holds
$\Gamma^{ev} \cap B_\tau \subset \{(r_v,\theta_v);\; 0 < \theta_v < \omega_v\}$.
Here, $\omega_v$ denotes the interior angle (on $\Gamma^{ev}$)
between the edges meeting at $v$.
For simplicity of notation we write out here only the leading singularities
in ${\mathbf u}^e$, ${\mathbf u}^v$, and ${\mathbf u}^{ev}$ on the face $\Gamma^{ev}$,
thus omitting the corresponding terms of higher regularity
(see \cite[Appendix~A]{BespalovH_NpB} for complete expansions).

For the edge singularities ${\mathbf u}^{e}$ one has
\begin{equation} \label{ue}
    {\mathbf u}^e =
            {\mathbf curl}_{\Gamma^{ev}}\Big(x_{e2}^{\gamma_1^e}\,|\log x_{e2}|^{s_1^e}\,b_{1}^e(x_{e1})\,
            \chi_1^e(x_{e1})\,\chi_2^e(x_{e2})\Big)
            +\,
            x_{e2}^{\gamma_2^e}\,|\log x_{e2}|^{s_2^e}\,{\mathbf b}_{2}^e(x_{e1})\,
            \chi_1^e(x_{e1})\,\chi_2^e(x_{e2}),
\end{equation}
where ${\mathbf curl}_{\Gamma^{ev}} = (\partial/\partial x_{e2},\, -\partial/\partial x_{e1})$
is the tangential vector curl operator ${\mathbf curl}_{\G}$ restricted to the face $\Gamma^{ev}$
(cf. \cite{BuffaC_01_TFI,BuffaC_01_TII}),
$\gamma_{1}^e,\;\gamma_{2}^e > \frac 12$, and $s_1^e,\,s_2^e\ge 0$ are
integers. Here, $\chi_1^e$, $\chi_2^e$ are $C^\infty$ cut-off functions
with $\chi_1^e=1$ in a certain distance to the end points of $e$
and $\chi_1^e=0$ in a neighbourhood of these vertices.
Moreover, $\chi_2^e=1$ for $0\le x_{e2}\le\delta_e$
and $\chi_2^e=0$ for $x_{e2}\ge 2\delta_e$ with some $\delta_e \in (0,\frac 12)$.
The functions $b_{1}^e\chi_1^e \in H^{m_1}(e)$ and ${\mathbf b}_{2}^e\chi_1^e \in {\mathbf H}^{m_2}(e)$
for $m_1$ and $m_2$ as large as required.

The vertex singularities ${\mathbf u}^v$ have the form
\begin{equation} \label{uv}
    {\mathbf u}^v =
    {\mathbf curl}_{\Gamma^{ev}}\Big(r_v^{\lambda_1^{v}}\, |\log r_v|^{q_1^v}\,
    \chi^v(r_v)\, \chi^v_{1}(\theta_v)\Big)
    +\,r_v^{\lambda_2^{v}}\, |\log r_v|^{q_2^v}\,
    \chi^v(r_v)\, \hbox{\mathversion{bold}$\chi$}^v_{2}(\theta_v),
\end{equation}
where $\lambda_1^v,\;\lambda_2^v > - \frac 12$ are real numbers,
$q_1^v,\,q_2^v \ge 0$ are integers,
$\chi^v$ is a $C^\infty$ cut-off function with
$\chi^v=1$ for $0\le r_v\le\tau_v$ and $\chi^v=0$ for $r_v\ge 2\tau_v$
with some $\tau_v \in (0,\frac 12)$.
The functions $\chi_1^v,\,\hbox{\mathversion{bold}$\chi$}_2^v$ are such that
$\chi_1^v \in H^{t_1}(0,\omega_v)$,
$\hbox{\mathversion{bold}$\chi$}_2^v \in {\mathbf H}^{t_2}(0,\omega_v)$
for $t_1$, $t_2$ as large as required.

For the combined edge-vertex singularity ${\mathbf u}^{ev}$ one has
\[
    {\mathbf u}^{ev} = {\mathbf u}_1^{ev} + {\mathbf u}_2^{ev},
\]
where
\begin{eqnarray} \label{u1ev}
     {\mathbf u}_1^{ev}
     & = &
     {\mathbf curl}_{\Gamma^{ev}}\Big(
     x_{e1}^{\lambda_1^v-\gamma_1^e} x_{e2}^{\gamma_1^e}
     |\log x_{e1}|^{\beta_1}|\log x_{e2}|^{\beta_2}\,
     \chi^v(r_v) \chi^{ev}(\theta_v)
     \Big)
     \cr\cr
     & &
     \qquad\qquad +\,
     x_{e1}^{\lambda_2^v-\gamma_2^e} x_{e2}^{\gamma_2^e}
     |\log x_{e1}|^{\beta_3}|\log x_{e2}|^{\beta_4}\,
     \chi^v(r_v) \chi^{ev}(\theta_v)
     \left(
     \begin{array}{c}
     0 \\
     1
     \end{array}
     \right)
\end{eqnarray}
and
\begin{equation} \label{u2ev}
    {\mathbf u}_2^{ev} =
    {\mathbf curl}_{\Gamma^{ev}}\Big(x_{e2}^{\gamma_1^e}\,|\log x_{e2}|^{s_1^e}\,b_{3}^e(x_{e1},x_{e2})\,
    \chi_2^e(x_{e2})\Big)
    +\,
    x_{e2}^{\gamma_2^e}\,|\log x_{e2}|^{s_2^e}\,{\mathbf b}_{4}^e(x_{e1},x_{e2})\,
    \chi_2^e(x_{e2}).
\end{equation}
Here, $\lambda_i^v$, $\gamma_i^e$, $s_i^e$ ($i=1,2$), $\chi^v$, and $\chi_2^{e}$
are as above, $\beta_k \ge 0$ ($k = 1\ldots,4$) are integers,
$\beta_1 +\beta_2 = s_1^e +q_1^v$, $\beta_3 +\beta_4 = s_2^e +q_2^v$
with $q_1^v$, $q_2^v$ being as in (\ref{uv}),
$\chi^{ev}$ is a $C^\infty$ cut-off function with
$\chi^{ev}=1$ for $0\le\theta_v\le\beta_v$ and $\chi^{ev}=0$ for
$\frac 32\beta_v\le\theta_v\le\omega_v$ for some
$\beta_v \in (0,\min\{\omega_v/2,\pi/8\}]$.
The functions $b_{3}^e$ and ${\mathbf b}_{4}^e$, when extended by zero onto
$\mbox{\rm I\kern-.18em R}^{2+}:=\{(x_{e1},x_{e2});\; x_{e2}>0\}$, lie in $H^{m_1}(\mbox{\rm I\kern-.18em R}^{2+})$ and
${\mathbf H}^{m_2}(\mbox{\rm I\kern-.18em R}^{2+})$, respectively, with $m_1$, $m_2$ as large as required.
Finally, the supports of ${\mathbf u}_1^{ev}$ and ${\mathbf u}_2^{ev}$ are subsets of the sector
$\bar S_{ev} = \{(r_v,\theta_v);\;
                 0 \le r_v \le 2\tau_v,\, 0 \le \theta_v \le \frac 32 \beta_v\}$.

\begin{remark} \label{rem_sing1}
{\rm (i)}
The exponents $\gamma_i^e$ ($i=1,2$) of the edge and vertex-edge singularities
in {\rm (\ref{ue})}, {\rm (\ref{u1ev})}, {\rm (\ref{u2ev})} satisfy
$\gamma_{i}^e > \frac 12$. However, for our approximation analysis below
it suffices to require that $\gamma_{i}^e > 0$ ($i=1,2$).
Note that $\gamma_{i}^e > 0$ and $\lambda_{i}^v > -\frac 12$ ($i=1,2$)
are the minimum requirements to guarantee ${\mathbf u} \in {\mathbf X}$.

{\rm (ii)}
As mentioned above, the terms of higher regularity (i.e., with greater
singularity exponents) are omitted in {\rm (\ref{ue})--(\ref{u2ev})}.
These terms are necessary
to obtain the regular part ${\mathbf u}_{\rm reg} \in {\mathbf X}^k$ of decomposition {\rm (\ref{dec})}
as smooth as required. This can be done by considering sufficiently many
(omitted) singularity terms of each type.
\end{remark}

\begin{remark} \label{rem_sing2}
{\rm (i)}
By {\rm (\ref{ue})--(\ref{u2ev})} we conclude that any singular vector
field ${\mathbf u}^s$ in {\rm (\ref{sing})} ($s = e,\, v$, or $ev$) can be written as
\begin{equation} \label{us}
    {\mathbf u}^s = {\mathbf curl}_{\G}\, w^s + {\mathbf v}^s = {\mathbf curl}_{\G}\, w^s + (v^s_1,\,v^s_2)
\end{equation}
with corresponding (scalar) singular functions $w^s,\,v^s_1,\,v^s_2$
being defined on the whole surface $\Gamma$. Note that these
scalar functions are $H^{1/2}$-regular on $\Gamma$, i.e.,
\[
  w^s \in H^{1/2}(\Gamma),\qquad {\mathbf v}^s = (v^s_1,\,v^s_2) \in {\mathbf H}^{1/2}(\Gamma),\qquad
  s = e,\, v,\, ev.
\]

{\rm (ii)}
It is important to observe that the functions $w^s,\,v_1^s,\,v_2^s$
($s = e,\, v$, or $ev$) in {\rm (\ref{us})} are typical scalar singularities
inherent to solutions of the boundary integral equations
with hypersingular operator for the Laplacian on $\Gamma$
and with possibly singular right-hand side.
Continuous piecewise polynomial approximations of these scalar singularities
in fractional-order Sobolev spaces were analysed
in~{\rm \cite{BespalovH_08_hpB,Bespalov_09_NPA}}
and will be used to prove the main result of the present paper.
\end{remark}

\section{The $hp$-version of the BEM and the main result} \label{sec_BEM}
\setcounter{equation}{0}

For the approximate solution of (\ref{bie_var}) we apply
the $hp$-version of the BEM on quasi-uniform triangulations of $\Gamma$.
Our BEM is based on Galerkin discretisations with an appropriate
family of ${\mathbf H}({\mathrm div}_{\G},\Gamma)$-conforming surface elements
(here, ${\mathbf H}({\mathrm div}_{\G},\Gamma) := \{{\mathbf u} \in {\mathbf L}^{2}(\Gamma);\; {\mathrm div}_{\G}\,{\mathbf u} \in L^{2}(\Gamma)\}$).
In what follows, $h > 0$ and $p \ge 1$ will always specify the mesh parameter and
a polynomial degree, respectively.
For any $\Omega \subset {\field{R}}^n$ we will denote
$\rho_{\Omega} = \sup \{{\mathrm diam} (B);\; \hbox{$B$ is a ball in $\Omega$}\}$.

Let ${\cal T} = \{\Delta_h\}$ be a family of meshes
$\Delta_h = \{\Gamma_j;\; j=1,\ldots,J\}$. Each mesh is a partition of $\Gamma$
into triangular elements $\Gamma_j$ such that $\bar\Gamma = \cup_{j=1}^{J} \bar\Gamma_j$,
and the intersection of any two elements $\bar\Gamma_j,\,\bar\Gamma_k$ ($j \not= k$)
is either a common vertex, an entire side, or empty.

In the following we always identify a face of the polyhedron $\Gamma$
with a subdomain of ${\field{R}}^2$.
We denote $h_j={\mathrm diam} (\Gamma_j)$ for any $\Gamma_j\in \Delta_h$.
Furthermore, any element $\Gamma_j$ is the image
of the reference triangle
$K = \{(\xi_1,\xi_2);\; 0 < \xi_1 <1,\ 0 < \xi_2 < 1 - \xi_1\}$
under an affine mapping $T_j$, more precisely
\[
  \bar\Gamma_j = T_j(\bar K),\quad
  {\mathbf x} = T_j(\hbox{\mathversion{bold}$\xi$}) =  B_j\,\hbox{\mathversion{bold}$\xi$} + {\mathbf b}_j,
\]
where $B_j \in {\field{R}}^{2\times 2}$,
${\mathbf b}_j \in {\field{R}}^2$, ${\mathbf x} = (x_1,x_2) \in \bar\Gamma_j$, and
$\hbox{\mathversion{bold}$\xi$} = (\xi_1,\xi_2) \in \bar K$.
Then, the Jacobian matrix of $T_j$ is $B_j \in {\field{R}}^{2\times 2}$,
and its determinant $J_j := \hbox{det}(B_j)$
satisfies the relation $|J_j| \simeq h_j^2$.

We consider a family ${\cal T}$ of quasi-uniform shape-regular meshes $\Delta_h$
on $\Gamma$ in the sense that there exist positive constants
$\sigma_1,\,\sigma_2$ independent of $h = \max\limits_j h_j$ such
that for any $\Gamma_j\in\Delta_h$ and arbitrary $\Delta_h\in {\cal T}$ there holds
\begin{equation} \label{mesh}
    h_j \le \sigma_1\,\rho_{\Gamma_j},\qquad
    h \le \sigma_2\,h_j.
\end{equation}
Whereas the mapping $T_j$ introduced above is used to associate
scalar functions defined on the real element $\Gamma_j$ and on the reference
triangle $K$, the Piola transformation is used to transform
vector-valued functions between $K$ and $\Gamma_j$:
\begin{equation} \label{Piola}
    {\mathbf v} = {\cal M}_j(\hat {\mathbf v}) = \hbox{$\frac {1}{J_j}$} B_j \hat{\mathbf v} \circ T_j^{-1},
    \quad
    \hat{\mathbf v} = {\cal M}_j^{-1}({\mathbf v}) = J_j B_j^{-1} {\mathbf v} \circ T_j.
\end{equation}

Let us introduce the needed polynomial sets.
By ${\cal P}_p(K)$ we denote the set of polynomials of total degree $\le p$
on the reference triangle $K$.
We will use two families of ${\mathbf H}({\mathrm div}_{\G},\Gamma)$-conforming surface elements:
the Raviart-Thomas (RT) and Brezzi-Douglas-Marini (BDM) elements.
The corresponding spaces of degree $p \ge 1$ on the reference triangle $K$
will be denoted as follows (see, e.g., \cite{BrezziF_91_MHF, RobertsT_91_MHM}):
\begin{eqnarray*}
     \hbox{\mathversion{bold}$\cal P$}^{\rm RT}_p(K) & = &
     ({\cal P}_{p-1}(K))^2 \oplus \hbox{\mathversion{bold}$\xi$} {\cal P}_{p-1}(K);
     \\[3pt]
     \hbox{\mathversion{bold}$\cal P$}^{\rm BDM}_p(K) & = &
     ({\cal P}_{p}(K))^2.
\end{eqnarray*}
We will use the unified notation $\hbox{\mathversion{bold}$\cal P$}_p(K)$ which refers to either
the RT- or BDM-space on $K$ for $p \ge 1$.
Accordingly, all results in this paper
are formulated in a unified way and are valid for both
the RT- and BDM-based boundary element spaces defined on the triangulation of $\Gamma$.
Note, however, that in some cases we will need to provide arguments separately for
each type of these boundary elements.

Using transformations (\ref{Piola}), we set
\begin{equation} \label{Xp}
    {\mathbf X}_{hp} := \{{\mathbf v} \in {\mathbf H}({\mathrm div}_{\G},\Gamma);\;
                      {\cal M}_j^{-1}({\mathbf v}|_{\Gamma_j}) \in \hbox{\mathversion{bold}$\cal P$}_p(K),\ j=1,\ldots,J\}.
\end{equation}
Note that only one type of surface elements (i.e., either the RT- or BDM-elements)
is used in (\ref{Xp}) for all triangles $\Gamma_j$.
We will denote by $N = N(h,p)$ the dimension of the discrete space ${\mathbf X}_{hp}$.
One has $N \simeq h^{-2}$ for fixed $p$ and $N \simeq p^2$ for fixed $h$.

The $hp$-version of the Galerkin BEM for the EFIE reads as:
{\em Find ${\mathbf u}_{hp}\in {\mathbf X}_{hp}$ such that}
\begin{equation} \label{BEM}
    a({\mathbf u}_{hp},{\mathbf v}) = \<{\mathbf f}, {\mathbf v}\> \quad \forall {\mathbf v} \in {\mathbf X}_{hp}.
\end{equation}

Due to the infinite-dimensional kernel of the ${\mathrm div}_{\G}$-operator,
the bilinear form $a(\cdot,\cdot)$ in (\ref{bie_var}) is not ${\mathbf X}$-coercive,
and, hence, the unique solvability of (\ref{BEM}) cannot be proved by standard
arguments. However, the refined analysis in \cite{BespalovHH_Chp}
shows that the unique BEM-solution ${\mathbf u}_{hp} \in {\mathbf X}_{hp}$ does exist,
and it converges quasi-optimally to the exact solution ${\mathbf u} \in {\mathbf X}$ of the EFIE
as $N(h,p) \to \infty$. This result is formulated in the following proposition.

\begin{prop} \label{prop_solve}
{\rm \cite[Theorem~1.2]{BespalovHH_Chp}}
There exists $N_0 \ge 1$ such that for any ${\mathbf f} \in {\mathbf X}'$ and for arbitrary
mesh-degree combination satisfying $N(h,p) \ge N_0$ the discrete problem {\rm (\ref{BEM})}
is uniquely solvable and the $hp$-version of the Galerkin BEM
converges quasi-optimally, i.e.,
\begin{equation} \label{quasi-optimality}
    \|{\mathbf u} - {\mathbf u}_{hp}\|_{{\mathbf X}} \le
    C \inf\{\|{\mathbf u} - {\mathbf v}\|_{{\mathbf X}};\; {\mathbf v}\in {\mathbf X}_{hp}\}.
\end{equation}
Here, ${\mathbf u} \in {\mathbf X}$ is the solution of {\rm (\ref{bie_var})},
${\mathbf u}_{hp} \in {\mathbf X}_{hp}$ is the solution of {\rm (\ref{BEM})},
$\|\cdot\|_{{\mathbf X}}$ denotes the norm in ${\mathbf X}$, and
$C>0$ is a constant independent of $h$ and $p$.
\end{prop}

The next theorem is the main result of this paper.
It states convergence rates of the $hp$-BEM with quasi-uniform meshes
for the EFIE. These convergence rates (in both the mesh parameter $h$
and the polynomial degree $p$) are given explicitly in terms of the
singularity exponents of the vector fields in (\ref{ue})--(\ref{u2ev}).

\begin{theorem} \label{thm_main_hp}
Let ${\mathbf u} \in {\mathbf X}$ and ${\mathbf u}_{hp} \in {\mathbf X}_{hp}$ be the solutions of {\rm (\ref{bie_var})}
and {\rm (\ref{BEM})}, respectively. We assume that the source functional ${\mathbf f}$
in {\rm (\ref{bie_var})} is sufficiently smooth such that representation
{\rm (\ref{dec})--(\ref{u2ev})} holds for the solution ${\mathbf u}$ of {\rm (\ref{bie_var})}.
Let $v_0\in V$, $e_0\in E(v_0)$ be a vertex-edge pair such that
\[
  \min\{\lambda_1^{v_0}+1/2, \lambda_2^{v_0}+1/2, \gamma_1^{e_0}, \gamma_2^{e_0}\} =
  \min_{v\in V, e\in E(v)}\min\,\{\lambda_1^v+1/2, \lambda_2^v+1/2, \gamma_1^e, \gamma_2^e\}
\]
with $\lambda_i^v$ and $\gamma_i^e$ ($i = 1,2$) being as in {\rm (\ref{ue})--(\ref{u2ev})}.
Then for any $h > 0$ and for every
$p \ge \min\,\{\lambda_1^{v_0},\lambda_2^{v_0},\gamma_1^{e_0} - \frac 12,\gamma_2^{e_0} - \frac 12\}$
there holds
\begin{equation} \label{hp_rates_1}
    \|{\mathbf u} - {\mathbf u}_{hp}\|_{{\mathbf X}} \le
    C\,\bigg(\frac{h}{p^{\,2}}\bigg)^{\min\{\lambda_1^{v_0}+1/2,\lambda_2^{v_0}+1/2,\gamma_1^{e_0},\gamma_2^{e_0}\}}
    \Big(1 + \log\frac ph\Big)^{\beta + \nu},
\end{equation}
where
\begin{equation} \label{beta}
    \beta =
    \cases{
           \max\,\{q_1^{v_0}+s_1^{e_0}+\frac 12,\, q_2^{v_0}+s_2^{e_0}+\frac 12\}
                       & \hbox{if \ $\lambda_i^{v_0} = \gamma_i^{e_0}-\frac 12$ for $i=1,2$},\cr
    \noalign{\vskip3pt}
           \max\,\{q_1^{v_0}+s_1^{e_0}+\frac 12,\, q_2^{v_0}+s_2^{e_0}\}
                       & \hbox{if \ $\lambda_1^{v_0} = \gamma_1^{e_0}-\frac 12$,
                                  \ $\lambda_2^{v_0} \not= \gamma_2^{e_0}-\frac 12$},\cr
    \noalign{\vskip3pt}
           \max\,\{q_1^{v_0}+s_1^{e_0},\, q_2^{v_0}+s_2^{e_0}+\frac 12\}
                       & \hbox{if \ $\lambda_1^{v_0} \not= \gamma_1^{e_0}-\frac 12$,
                                  \ $\lambda_2^{v_0} = \gamma_2^{e_0}-\frac 12$},\cr
    \noalign{\vskip3pt}
           \max\,\{q_1^{v_0}+s_1^{e_0},\, q_2^{v_0}+s_2^{e_0}\}
                       & \hbox{otherwise}\cr
          }
\end{equation}
with the numbers $s_i^{e_0}$ and $q_i^{v_0}$ ($i=1,2$) given in
{\rm(\ref{ue})} and {\rm(\ref{uv})}, respectively, and
\begin{equation} \label{nu}
    \nu = \cases{
                 \frac 12
                 & \hbox{if $p = \min\,\{\lambda_1^{v_0},\lambda_2^{v_0},\gamma_1^{e_0}-1/2,\gamma_2^{e_0}-1/2\}$},\cr
                 \noalign{\vskip3pt}
                 0
                 & \hbox{otherwise}.\cr
                }
\end{equation}
If $1 \le p < \min\,\{\lambda_1^{v_0},\lambda_2^{v_0},\gamma_1^{e_0} - \frac 12,\gamma_2^{e_0} - \frac 12\}$,
then for any $h > 0$ there holds
\begin{equation} \label{hp_rates_2}
    \|{\mathbf u} - {\mathbf u}_{hp}\|_{{\mathbf X}} \le C\,h^{p+1/2}.
\end{equation}
\end{theorem}

\begin{proof}
Considering enough singularity terms in representation (\ref{sing})
the function ${\mathbf u}_{\rm reg}$ in (\ref{reg}) is as regular as needed.
Then, due to the quasi-optimal convergence (\ref{quasi-optimality})
of the $hp$-BEM with quasi-uniform meshes, the assertion follows immediately
from the general approximation result given in Theorem~\ref{thm_gen_approx} below.
\end{proof}

\begin{remark} \label{rem_parallelograms}
We have only considered meshes of triangular elements on $\Gamma$.
If the meshes contain also shape-regular parallelogram elements (i.e., affine images
of the reference square), then a priori error estimates of
Theorem~{\rm \ref{thm_main_hp}} remain valid only in the case of the RT-based BEM.
This is because all the auxiliary results needed for the proof are
valid in this case. This, however, is not true for the BDM-based BEM.
In particular, the arguments in the proofs of Proposition~{\rm \ref{prop_solve}} above
and Proposition~{\rm \ref{prop_approx_reg}} below rely essentially
on the fact that the involved polynomial spaces
form the exact ${\mathbf curl}\,{-}\,\div$ sequence (on the reference element),
a property the BDM-spaces fail to satisfy on the reference
square (see, e.g., {\rm \cite{Demkowicz_08_PES}}).
\end{remark}

\begin{remark} \label{rem_open_surface}
We have assumed that $\Gamma$ is a polyhedral (closed) surface.
However, all arguments in our proofs carry over only with minor modifications
to the case of a piecewise plane orientable open surface $\Gamma$. Note that in this case
there are no restrictions needed on the wave number $\kappa$ to ensure the uniqueness
of the solution to the EFIE; the strongest edge singularities in
{\rm (\ref{ue})} have the exponents $\gamma_i^e = \frac 12$ $(i=1,2)$;
the energy space for the EFIE is $\tilde{\mathbf H}^{-1/2}_0({\mathrm div}_{\G},\Gamma)$ and
the boundary element space ${\mathbf X}_{hp}$ consists of ${\mathbf H}({\mathrm div}_{\G},\Gamma)$-conforming
polynomial vector fields with normal components vanishing on $\partial\Gamma$
(see {\rm \cite{BespalovH_NpB,BespalovH_10_hpA}} for details).
\end{remark}

\section{General approximation result} \label{sec_gen_approx}
\setcounter{equation}{0}

In this section we prove the following general $hp$-approximation result for the
vector field ${\mathbf u}$ given by formulas (\ref{dec})--(\ref{u2ev})
with the singularity exponents $\gamma_{i}^e$ and $\lambda_{i}^v$ ($i=1,2$)
satisfying the minimum requirements to guarantee ${\mathbf u} \in {\mathbf X}$.

\begin{theorem} \label{thm_gen_approx}
Let the vector field ${\mathbf u}$ be given by {\rm (\ref{dec})--(\ref{u2ev})} on $\Gamma$
with $\gamma_1^e,\,\gamma_2^e>0$ and $\lambda_1^v,\,\lambda_2^v>-\frac 12$
for each edge $e$ and every vertex $v$.
Also, let $v_0\in V$, $e_0\in E(v_0)$ be such that
\[
  \min\{\lambda_1^{v_0}+1/2, \lambda_2^{v_0}+1/2, \gamma_1^{e_0}, \gamma_2^{e_0}\} =
  \min_{v\in V, e\in E(v)}\min\,\{\lambda_1^v+1/2, \lambda_2^v+1/2, \gamma_1^e, \gamma_2^e\}
\]
with $\lambda_i^v$ and $\gamma_i^e$ ($i = 1,2$) being as in {\rm (\ref{ue})--(\ref{u2ev})}.
Then for any $h > 0$ and for every
$p \ge \min\,\{\lambda_1^{v_0},\lambda_2^{v_0},\gamma_1^{e_0} - \frac 12,\gamma_2^{e_0} - \frac 12\}$,
there exists ${\mathbf u}^{hp} \in {\mathbf X}_{hp}$ such that
\begin{eqnarray}
     \|{\mathbf u} - {\mathbf u}^{hp}\|_{{\mathbf X}}
     & \hskip-4pt \le \hskip-4pt &
     C\,\max\left\{
     h^{\min\{k,p\}+1/2}\,p^{-(k+1/2)},\right.
     \nonumber
     \\
     &     &
       \qquad\qquad\left.
       \bigg(\frac{h}{p^{\,2}}\bigg)^{\min\{\lambda_1^{v_0}+1/2,\lambda_2^{v_0}+1/2,\gamma_1^{e_0},\gamma_2^{e_0}\}}
       \Big(1 + \log\frac ph\Big)^{\beta + \nu}
     \right\},
     \label{gen_approx_1}
\end{eqnarray}
where $\beta$ and $\nu$ are defined by {\rm (\ref{beta})} and {\rm (\ref{nu})}, respectively.

If $1 \le p < \min\,\{\lambda_1^{v_0},\lambda_2^{v_0},\gamma_1^{e_0} - \frac 12,\gamma_2^{e_0} - \frac 12\}$,
then for any $h > 0$ there exists ${\mathbf u}^{hp} \in {\mathbf X}_{hp}$ such that
\begin{equation} \label{gen_approx_2}
    \|{\mathbf u} - {\mathbf u}^{hp}\|_{{\mathbf X}} \le C\,h^{\min\{k,p\}+1/2}.
\end{equation}
\end{theorem}

In order to prove this theorem one needs to find discrete vector fields belonging to ${\mathbf X}_{hp}$
and approximating the smooth and singular parts of ${\mathbf u}$ such that the approximation errors
satisfy the upper bounds in (\ref{gen_approx_1}) and (\ref{gen_approx_2}).

We start with formulating the following $hp$-approximation result
for regular vector fields on $\Gamma$. This result will be used, in particular, to approximate
the vector field ${\mathbf u}_{\rm reg} \in {\mathbf X}^k$ in (\ref{dec}).

\begin{prop} \label{prop_approx_reg}
Let $P_{hp}:\, {\mathbf X} \rightarrow {\mathbf X}_{hp}$ be the orthogonal projection
with respect to the norm in ${\mathbf X}$. If ${\mathbf u} \in {\mathbf X}^k$ with
$k > 0$, then
\begin{equation} \label{approx_reg}
    \|{\mathbf u} - P_{hp}{\mathbf u}\|_{{\mathbf X}} \le
    C\,h^{\min\{k,p\}+1/2}\,p^{-(k+1/2)}\,\|{\mathbf u}\|_{{\mathbf X}^k}
\end{equation}
with a positive constant $C$ independent of $h$, $p$, and ${\mathbf u}$.
\end{prop}

The proof is given in \cite[Theorem~4.1]{BespalovH_10_hpA} for the case
of RT-spaces, and it carries over without
essential modifications to the case of BDM-spaces on {\em triangular}
elements (cf. Remark~\ref{rem_parallelograms}).

Now, we will study approximations of the singular part ${\mathbf u}_{\rm sing}$
in representation (\ref{dec}). By (\ref{sing})--(\ref{u2ev}) and due to the arguments
in Remark~\ref{rem_sing2}~(i), we conclude that ${\mathbf u}_{\rm sing}$ can be written as
\begin{equation} \label{dec_sing}
    {\mathbf u}_{\rm sing} = {\mathbf curl}_{\G}\, w + {\mathbf v} = {\mathbf curl}_{\G}\, w + (v_1,\,v_2),
\end{equation}
where $w \in H^{1/2}(\Gamma)$ and ${\mathbf v} \in {\mathbf H}^{1/2}(\Gamma)$.

Let us define the following discrete space (of continuous piecewise polynomials)
over the mesh $\Delta_h$:
\[
  S_{hp}(\Gamma) := \{v\in C^0(\Gamma);\; v|_{\Gamma_j} \circ T_j \in {\cal P}_p(K),\ j=1,\ldots,J\}.
\]
We will also need the following functions of $h$ and $p$:
\begin{equation} \label{f_j(hp)}
    f_j(h,p) :=
    h^{\alpha_j} p^{-2 \alpha_j} (1 + \log(p/h))^{\tilde\beta_j+\nu_j},
\end{equation}
where $j = 1,2$,
\begin{equation} \label{alpha_j}
    \alpha_j := \min\,\{\lambda_j^{v_0}+1/2,\gamma_j^{e_0}\},
\end{equation}
\begin{equation} \label{beta_j}
    \tilde\beta_j := \cases{
                            q_j^{v_0} + s_j^{e_0} +\frac 12
                            & \hbox{if $\lambda_j^{v_0} = \gamma_j^{e_0} - \frac 12$},\cr
                     \noalign{\vskip3pt}
                            q_j^{v_0} + s_j^{e_0}
                            & \hbox{otherwise},\cr
                           }
\end{equation}
\begin{equation} \label{nu_j}
    \nu_j := \cases{
                    \frac 12
                    & \hbox{if $p = \alpha_j - \frac 12$},\cr
                     \noalign{\vskip3pt}
                    0
                    & \hbox{otherwise},\cr
                   }
\end{equation}
and the numbers $\gamma_j^{e_0},\, \lambda_j^{v_0},\, s_j^{e_0},\, q_j^{v_0}$ ($j = 1,2$)
are given in (\ref{ue})--(\ref{u2ev}) for the vertex-edge pair $(v_0,e_0)$
introduced in the formulation of Theorem~\ref{thm_gen_approx}.

In the next two lemmas we formulate approximation results for the
vector fields ${\mathbf curl}_{\G}\, w$ and ${\mathbf v}$ on the right-hand side of (\ref{dec_sing}).
The proofs are given in Section~\ref{sec_proofs} below.

\begin{lemma} \label{lm_approx_w}
Let $w \in H^{1/2}(\Gamma)$ be the scalar singular function in representation {\rm (\ref{dec_sing})}.
Then for any $h > 0$ and $p \ge 1$, there exists $w^{hp} \in S_{hp}(\Gamma)$
such that ${\mathbf curl}_{\G}\, w^{hp} \in {\mathbf X}_{hp}$ and there holds
\begin{equation} \label{approx_w}
    \|{\mathbf curl}_{\G}\,w - {\mathbf curl}_{\G}\,w^{hp}\|_{{\mathbf X}} \le
    \cases{
           C\,f_1(h,p)
           & \hbox{if \ $p \ge \alpha_1 - \frac 12$},\cr
           \noalign{\vskip3pt}
           C\,h^{p+1/2}
           & \hbox{if \ $1 \le p < \alpha_1 - \frac 12$},\cr
          }
\end{equation}
where $f_1(h,p)$ and $\alpha_1$ are defined by {\rm (\ref{f_j(hp)})}
and {\rm (\ref{alpha_j})}, respectively.
\end{lemma}

\begin{lemma} \label{lm_approx_v}
Let ${\mathbf v} \in {\mathbf H}^{1/2}(\Gamma)$ be the singular vector field in representation {\rm (\ref{dec_sing})}.
Then for any $h > 0$ and $p \ge 1$, there exists ${\mathbf v}^{hp} \in {\mathbf X}_{hp}$
satisfying
\begin{equation} \label{approx_v}
    \|{\mathbf v} - {\mathbf v}^{hp}\|_{{\mathbf X}} \le
    \cases{
           C\,f_2(h,p)
           & \hbox{if \ $p \ge \alpha_2 - \frac 12$},\cr
           \noalign{\vskip3pt}
           C\,h^{p+1/2}
           & \hbox{if \ $1 \le p < \alpha_2 - \frac 12$},\cr
          }
\end{equation}
where $f_2(h,p)$ and $\alpha_2$ are defined by {\rm (\ref{f_j(hp)})}
and {\rm (\ref{alpha_j})}, respectively.
\end{lemma}

Now we are able to prove Theorem~\ref{thm_gen_approx}.
\medskip

\noindent{\bf Proof of Theorem~\ref{thm_gen_approx}.}
For the regular vector field ${\mathbf u}_{\rm reg} \in {\mathbf X}^k$ in (\ref{dec}) we use the
orthogonal projection $P_{hp}:\,{\mathbf X} \rightarrow {\mathbf X}_{hp}$
with respect to the norm in ${\mathbf X}$ to define
${\mathbf u}_{\rm reg}^{hp} := P_{hp} {\mathbf u}_{\rm reg} \in {\mathbf X}_{hp}$.
Then we have by Proposition~\ref{prop_approx_reg}
\begin{equation} \label{approx_reg_1}
    \|{\mathbf u}_{\rm reg} - {\mathbf u}_{\rm reg}^{hp}\|_{{\mathbf X}} \le
    C\,h^{\min\{k,p\}+1/2}\,p^{-(k+1/2)}.
\end{equation}
Since the singular part of decomposition (\ref{dec}) can be written as in (\ref{dec_sing}),
we use approximations $w^{hp} \in S_{hp}(\Gamma)$ and ${\mathbf v}^{hp} \in {\mathbf X}_{hp}$
from the above two lemmas to define
${\mathbf u}_{\rm sing}^{hp} := {\mathbf curl}_{\G}\, w^{hp} + {\mathbf v}^{hp} \in {\mathbf X}_{hp}$.
Then, applying the triangle inequality, we obtain by (\ref{approx_w}) and (\ref{approx_v})
\begin{equation} \label{approx_sing}
    \|{\mathbf u}_{\rm sing} - {\mathbf u}_{\rm sing}^{hp}\|_{{\mathbf X}} \le
    \cases{
           C\,\max \left\{f_1(h,p),\;f_2(h,p)\right\}
           & \hbox{if \ $p \ge \min\,\{\alpha_1,\alpha_2\} - \frac 12$},\cr
           \noalign{\vskip3pt}
           C\,h^{p+1/2}
           & \hbox{if \ $1 \le p < \min\,\{\alpha_1,\alpha_2\} - \frac 12$}.\cr
          }
\end{equation}
Now, we set ${\mathbf u}^{hp} := {\mathbf u}_{\rm reg}^{hp} + {\mathbf u}_{\rm sing}^{hp} \in {\mathbf X}_{hp}$.
Combining estimates (\ref{approx_reg_1}) and (\ref{approx_sing}),
applying the triangle inequality, and using expressions
(\ref{f_j(hp)}) and (\ref{alpha_j}) for the functions $f_j(h,p)$
and the parameters $\alpha_j$, respectively, we prove the desired estimates in
(\ref{gen_approx_1}) and (\ref{gen_approx_2}).\nopagebreak\hbox{ }\hfill$\Box$\bigskip

\section{Proofs of technical lemmas} \label{sec_proofs}
\setcounter{equation}{0}

In this section we prove Lemmas~\ref{lm_approx_w} and~\ref{lm_approx_v}.

\subsection{Proof of Lemma~\ref{lm_approx_w}} \label{sec_proof_w}

Recalling Remark~\ref{rem_sing2}~(ii), we use
the results of \cite{BespalovH_08_hpB,Bespalov_09_NPA}
to find the desired piecewise polynomial $w^{hp} \in S_{hp}(\Gamma)$
such that the norm $\|w - w^{hp}\|_{H^{1/2}(\Gamma)}$ is bounded as
in (\ref{approx_w}) (see Theorems~5.1,~5.2 in \cite{BespalovH_08_hpB}
and Theorem~4.1 in \cite{Bespalov_09_NPA}). Then, recalling the fact that the
the operator ${\mathbf curl}_{\G}: H^{1/2}(\Gamma) \rightarrow {\mathbf H}^{-1/2}_{\|}(\Gamma)$
is continuous (see \cite{BuffaC_01_TII}), we derive the estimate in (\ref{approx_w}):
\begin{eqnarray*}
     \|{\mathbf curl}_{\G}\,w - {\mathbf curl}_{\G}\,w^{hp}\|_{{\mathbf X}}
     & = &
     \|{\mathbf curl}_{\G} (w - w^{hp})\|_{{\mathbf H}^{-1/2}_{\|}(\Gamma)}
     \\[3pt]
     & \le &
     C \|w - w^{hp}\|_{H^{1/2}(\Gamma)} \le
    \cases{
           C\,f_1(h,p)
           & \hbox{if \ $p \ge \alpha_1 - \frac 12$},\cr
           \noalign{\vskip3pt}
           C\,h^{p+1/2}
           & \hbox{if \ $1 \le p < \alpha_1 - \frac 12$}.\cr
          }
\end{eqnarray*}
It remains to prove that ${\mathbf curl}_{\G}\, w^{hp} \in {\mathbf X}_{hp}$.
In fact, it is easy to check (see~\cite[p.~615]{BespalovH_NpB}) that
\[
     {\cal M}_j^{-1}({\mathbf curl}_{\G}\, w^{hp}|_{\Gamma_j}) =
     {\mathbf curl}\, (w^{hp}|_{\Gamma_j} \circ T_j),\qquad
     \hbox{where
           ${\mathbf curl} = \Big(
                          \hbox{$\frac{\partial}{\partial \xi_{2}}$},\,
                          -\hbox{$\frac{\partial}{\partial \xi_{1}}$}
                     \Big)$.
          }
\]
Hence, recalling that $w^{hp}|_{\Gamma_j} \circ T_j \in {\cal P}_p(K)$,
we conclude that
${\cal M}_j^{-1}({\mathbf curl}_{\G}\, w^{hp}|_{\Gamma_j}) \in
 ({\cal P}_{p-1}(K))^2 \subset \hbox{\mathversion{bold}$\cal P$}_p(K)$
(this is because $({\cal P}_{p-1}(K))^2$ is a subset of both $\hbox{\mathversion{bold}$\cal P$}^{\rm BDM}_p(K)$ and
$\hbox{\mathversion{bold}$\cal P$}^{\rm RT}_p(K)$).
Moreover, ${\mathbf curl}_{\G}\, w^{hp} \in {\mathbf H}({\mathrm div}_{\G},\Gamma)$,
because ${\mathrm div}_{\G}({\mathbf curl}_{\G}\, w^{hp}) \equiv 0$ on $\Gamma$.
Therefore, ${\mathbf curl}_{\G}\, w^{hp} \in {\mathbf X}_{hp}$, and the proof is finished.

\subsection{Proof of Lemma~\ref{lm_approx_v}} \label{sec_proof_v}

Let ${\mathbf v} = (v_1,v_2)$ be the second term in decomposition (\ref{dec_sing}) of the
singular vector field ${\mathbf u}_{\rm sing}$. Again, using the results
in \cite{BespalovH_08_hpB,Bespalov_09_NPA}, we find continuous piecewise polynomial
approximations to the scalar components $v_1,\ v_2$ of ${\mathbf v}$:
for any $h > 0$ and every $p \ge 1$,
there exist $v_1^{hp},\,v_2^{hp} \in S_{hp}(\Gamma)$ such that
for $i = 1,2$ there holds
\begin{equation} \label{v_approx1}
    \|v_i - v_{i}^{hp}\|_{H^{1/2}(\Gamma)} \le
    \cases{
           C\,f_2(h,p)
           & \hbox{if \ $p \ge \alpha_2 - \frac 12$},\cr
           \noalign{\vskip3pt}
           C\,h^{p+1/2}
           & \hbox{if \ $1 \le p < \alpha_2 - \frac 12$}\cr
          }
\end{equation}
with positive constants $C > 0$ independent of $h$ and $p$.

Let ${\mathbf v}^{hp} = (v_1^{hp},\,v_2^{hp})$. Then ${\mathbf v}^{hp} \in {\mathbf H}({\mathrm div}_{\G},\Gamma)$.
Furthermore, in the case of BDM-elements we observe that
for any $\Gamma_j$ there holds
\[
  {\cal M}_j^{-1}({\mathbf v}^{hp}|_{\Gamma_j}) =
  J_j B_j^{-1} ({\mathbf v}^{hp}|_{\Gamma_j}) \circ T_j \in
  ({\cal P}_{p}(K))^2 = \hbox{\mathversion{bold}$\cal P$}^{\rm BDM}_p(K).
\]
Therefore, ${\mathbf v}^{hp} \in {\mathbf X}_{hp}$ in this case.
Moreover, since ${\mathbf v} \in {\mathbf H}^{1/2}(\Gamma)$ and $v^{hp}_i \in S_{hp}(\Gamma)$,
we use estimate (\ref{v_approx1}) and the continuity of the operator
${\mathrm div}_{\G}: {\mathbf H}^{1/2}(\Gamma) \rightarrow H^{-1/2}(\Gamma)$ to obtain
\begin{eqnarray}
    \|{\mathbf v} - {\mathbf v}^{hp}\|_{{\mathbf X}}
    & = &
    \|{\mathbf v} - {\mathbf v}^{hp}\|_{{\mathbf H}^{-1/2}_{\|}(\Gamma)} +
    \|{\mathrm div}_{\G} ({\mathbf v} - {\mathbf v}^{hp})\|_{H^{-1/2}(\Gamma)}
    \nonumber
    \\
    & \le &
    C \|{\mathbf v} - {\mathbf v}^{hp}\|_{{\mathbf H}^{1/2}(\Gamma)} \le
    C \sum\limits_{i=1}^{2} \|v_i - v_i^{hp}\|_{H^{1/2}(\Gamma)}
    \nonumber
    \\
    & \le &
    \cases{
           C\,f_2(h,p)
           & \hbox{if \ $p \ge \alpha_2 - \frac 12$},\cr
           \noalign{\vskip3pt}
           C\,h^{p+1/2}
           & \hbox{if \ $1 \le p < \alpha_2 - \frac 12$}.\cr
          }
    \label{v_approx2}
\end{eqnarray}

Unfortunately, in the case of RT-elements this component-wise approximation
of ${\mathbf v}$ does not work since the dimension of the RT-space on the reference triangle
is smaller than the dimension of the BDM-space, and, in general,
${\cal M}_j^{-1}({\mathbf v}^{hp}|_{\Gamma_j}) \notin \hbox{\mathversion{bold}$\cal P$}^{\rm RT}_p(K)$, so that
${\mathbf v}^{hp} = (v_1^{hp},\,v_2^{hp}) \notin {\mathbf X}_{hp}$.
In this case we follow the procedure described in \cite{BespalovH_08_hpB} (for the scalar case).
More precisely, we use appropriate
$h$-scaled cut-off functions and represent the vector field ${\mathbf v}$ as the sum of
a singular vector field $\hbox{\mathversion{bold}$\varphi$}$ with small support in the vicinity of the edges
and a sufficiently smooth vector field $\hbox{\mathversion{bold}$\psi$}$.
In the rest of this subsection we will demonstrate
how this procedure works in the vector case. We will give a concise step-by-step
outline of the procedure referring frequently to \cite{BespalovH_08_hpB}
for particular error estimates and other technical details.

{\bf Step 1: decomposition of ${\mathbf v}$.}
Let $h_0 = (\sigma_1 \sigma_2)^{-1} h$ with $\sigma_1,\,\sigma_2$ from
(\ref{mesh}). Using the scaled cut-off functions
$\chi^e_2(x_{e2}/h_0)$ and $\chi^v(r_v/h_0)$ with $\chi^e_2$ and $\chi^v$
from (\ref{ue}) and (\ref{uv}), respectively, one can decompose ${\mathbf v}$
as follows (cf. \cite[eqs. (5.4), (5.21), (6.4)]{BespalovH_08_hpB})
\begin{equation} \label{v_approx3}
    {\mathbf v} = \hbox{\mathversion{bold}$\varphi$} + \hbox{\mathversion{bold}$\psi$},
\end{equation}
where ${\mathrm supp}\,\hbox{\mathversion{bold}$\varphi$} \subset \cup_{v \in V} \cup_{e \in E(v)} (\bar A_e \cup \bar A_v)$,
$\hbox{\mathversion{bold}$\varphi$} \in {\mathbf H}^{1/2}(\Gamma)$, $\hbox{\mathversion{bold}$\psi$} \in {\mathbf X}^m$ (see (\ref{reg}) for the notation), and
$\hbox{\mathversion{bold}$\psi$}$ vanishes in small ($h$-dependent) neighbourhoods
of each vertex and each edge of~$\Gamma$. Here, $A_e$ is the union of elements at one edge
$e$, i.e., $\bar A_{e} := \cup\{\bar\Gamma_j;\; \bar\Gamma_j \cap e \not= \mbox{\o}\}$
(note that the endpoints of $e$ are not included in $e$),
$A_v$ is the union of elements at a vertex $v$, i.e.,
$\bar A_{v} := \cup\{\bar\Gamma_j;\; v \in \bar\Gamma_j\}$,
$m$ is sufficiently large and depends
on the parameters $t_2$ and $m_2$ specified for the singularities ${\mathbf u}^v$
and ${\mathbf u}_2^{ev}$, respectively.

{\bf Step 2: approximation of $\hbox{\mathversion{bold}$\varphi$}$ for $p = 1$.}
If $p = 1$, then one can approximate $\hbox{\mathversion{bold}$\varphi$}$ by zero.
One has $\hbox{\mathversion{bold}$\varphi$}_{hp} \equiv {\mathbf 0} \in {\mathbf X}_{hp}$, and, recalling the continuity
of the operator ${\mathrm div}_{\G}$, we derive
\begin{equation} \label{v_approx4}
    \|\hbox{\mathversion{bold}$\varphi$} - \hbox{\mathversion{bold}$\varphi$}_{hp}\|_{{\mathbf X}} =
    \|\hbox{\mathversion{bold}$\varphi$}\|_{{\mathbf H}^{-1/2}_{\|}(\Gamma)} + \|{\mathrm div}_{\G}\,\hbox{\mathversion{bold}$\varphi$}\|_{H^{-1/2}(\Gamma)} \le
    C \|\hbox{\mathversion{bold}$\varphi$}\|_{{\mathbf H}^{1/2}(\Gamma)}.
\end{equation}
We will obtain $h$-estimates for the norms of $\hbox{\mathversion{bold}$\varphi$}$ in ${\mathbf H}^s(\Gamma)$
with $s=0$ and $s = \frac 12 + \varepsilon$ (for sufficiently small $\varepsilon > 0$)
by using the fact that $\hbox{\mathversion{bold}$\varphi$}$ has a small support.
First, we apply Lemma~3.1 of \cite{BespalovH_05_pBE} and Lemma~3.5 of \cite{BespalovH_08_hpB}
to localise the norm $\|\hbox{\mathversion{bold}$\varphi$}\|_{{\mathbf H}^{1/2+\varepsilon(\Gamma)}}$ to the faces $\Gamma^f \subset \Gamma$
and to the elements $\Gamma_j \subset \Gamma^f$, respectively. Then, we use scaling
on each element $\Gamma_j \subset {\mathrm supp}\,\hbox{\mathversion{bold}$\varphi$}$
(cf. \cite[Lemma~3.1, eqs. (5.6), (5.11), (5.23), (6.7)]{BespalovH_08_hpB}).
As a result we have
\begin{eqnarray*}
     \|\hbox{\mathversion{bold}$\varphi$}\|^2_{{\mathbf H}^s(\Gamma)}
     & \le &
     C \sum\limits_{f:\, \Gamma^f \subset \Gamma} \|\hbox{\mathversion{bold}$\varphi$}\|^2_{{\mathbf H}^s(\Gamma^f)}
     \\
     & \le &
     C \sum\limits_{f:\, \Gamma^f \subset \Gamma}\
       \sum\limits_{j:\, \Gamma_j \subset \Gamma^f}
       \bigg(h^{-2s} \|\hbox{\mathversion{bold}$\varphi$}\|^2_{{\mathbf L}^2(\Gamma_j)} + |\hbox{\mathversion{bold}$\varphi$}|^2_{{\mathbf H}^s(\Gamma_j)}\bigg)
     \\[3pt]
     & \le &
     C h^{2(\alpha_2 + 1/2 - s)} (1+\log(1/h))^{2\tilde\beta_2}\qquad
     \hbox{for $s \in \{0, 1/2 + \varepsilon\}$}.
\end{eqnarray*}
Here, $\alpha_2$ and $\tilde\beta_2$ are defined  by (\ref{alpha_j})
and (\ref{beta_j}), respectively.
Hence, using the interpolation between ${\mathbf H}^0(\Gamma)$ and ${\mathbf H}^{1/2+\varepsilon}(\Gamma)$,
we obtain by (\ref{v_approx4})
\begin{equation} \label{v_approx5}
    \|\hbox{\mathversion{bold}$\varphi$} - \hbox{\mathversion{bold}$\varphi$}_{hp}\|_{{\mathbf X}} \le
     C h^{\alpha_2} (1+\log(1/h))^{\tilde\beta_2}.
\end{equation}

{\bf Step 3: approximation of $\hbox{\mathversion{bold}$\varphi$}$ for $p \ge 2$.}
We approximate each component $\varphi_i$ of $\hbox{\mathversion{bold}$\varphi$}$ by a piecewise
polynomial $\varphi_i^{hp} \in S_{h,p-1}(\Gamma)$ ($i = 1,2$).
Here we can use the results of \cite{BespalovH_08_hpB,Bespalov_09_NPA} for
each type of singularity: there exist $\varphi_i^{hp} \in S_{h,p-1}(\Gamma)$ such that
for $i = 1,2$ there holds (cf. \cite[eqs. (5.12), (5.13), (5.24), (6.9)]{BespalovH_08_hpB})
\begin{equation} \label{v_approx6}
    \|\varphi_i - \varphi_{i}^{hp}\|_{H^{1/2}(\Gamma)} \le
    C\,h^{\alpha_2}\,(p-1)^{-2\alpha_2}\,(1+\log\hbox{$\frac{p-1}{h}$})^{\tilde\beta_2} \le
    C\,h^{\alpha_2}\,p^{-2\alpha_2}\,(1+\log\hbox{$\frac{p}{h}$})^{\tilde\beta_2},
\end{equation}
where $\alpha_2$ and $\tilde\beta_2$ are the same as in (\ref{v_approx5}).

Setting $\hbox{\mathversion{bold}$\varphi$}^{hp} = (\varphi_1^{hp},\varphi_2^{hp})$ we observe that
$\hbox{\mathversion{bold}$\varphi$}^{hp} \in {\mathbf H}({\mathrm div}_{\G},\Gamma)$, and for any element $\Gamma_j$ one has
\[
  {\cal M}_j^{-1}(\hbox{\mathversion{bold}$\varphi$}^{hp}|_{\Gamma_j}) =
  J_j B_j^{-1} (\hbox{\mathversion{bold}$\varphi$}^{hp}|_{\Gamma_j}) \circ T_j \in
  ({\cal P}_{p-1}(K))^2 \subset \hbox{\mathversion{bold}$\cal P$}^{\rm RT}_p(K).
\]
Hence, $\hbox{\mathversion{bold}$\varphi$}^{hp} \in {\mathbf X}_{hp}$. Moreover, since $\hbox{\mathversion{bold}$\varphi$} \in {\mathbf H}^{1/2}(\Gamma)$
and $\varphi_i^{hp} \in S_{h,p-1}(\Gamma)$ for $i=1,2$,
we estimate by analogy with (\ref{v_approx4}) and using (\ref{v_approx6}):
\begin{eqnarray} \label{v_approx7}
     \|\hbox{\mathversion{bold}$\varphi$} - \hbox{\mathversion{bold}$\varphi$}_{hp}\|_{{\mathbf X}}
     & \le &
     C \|\hbox{\mathversion{bold}$\varphi$} - \hbox{\mathversion{bold}$\varphi$}_{hp}\|_{{\mathbf H}^{1/2}(\Gamma)}
     \nonumber
     \\
     & \le &
     C \sum\limits_{i=1}^{2} \|\varphi_i - \varphi_i^{hp}\|_{H^{1/2}(\Gamma)} \le
     C\,h^{\alpha_2}\,p^{-2\alpha_2}\,(1+\log\hbox{$\frac{p}{h}$})^{\tilde\beta_2}.
\end{eqnarray}
Comparing (\ref{v_approx5}) and (\ref{v_approx7}) we observe that
one can use estimate (\ref{v_approx7}) also in the case $p = 1$.

{\bf Step 4: approximation of $\hbox{\mathversion{bold}$\psi$}$.}
Recalling that $\hbox{\mathversion{bold}$\psi$} \in {\mathbf X}^m$ we apply Proposition~\ref{prop_approx_reg}:
there exists $\hbox{\mathversion{bold}$\psi$}^{hp} \in {\mathbf X}_{hp}$ such that
\begin{equation} \label{v_approx8}
    \|\hbox{\mathversion{bold}$\psi$} - \hbox{\mathversion{bold}$\psi$}^{hp}\|_{{\mathbf X}} \le
    C\,h^{\min\{k,p\}+1/2}\,p^{-(k+1/2)}\,\|\hbox{\mathversion{bold}$\psi$}\|_{{\mathbf X}^k},\qquad
    0 \le k \le m.
\end{equation}
We also recall that $\hbox{\mathversion{bold}$\psi$}$ vanishes in the $h$-neighbourhood
of all edges and vertices of $\Gamma$. Therefore, having explicit expressions
of the singular vector field ${\mathbf v}$ (and thus, of $\hbox{\mathversion{bold}$\psi$}$) on each face $\Gamma^f \subset \Gamma$,
one can derive upper bounds for the norms $\|\hbox{\mathversion{bold}$\psi$}\|_{{\mathbf H}^k(\Gamma^f)}$ and
$\|\div_{\Gamma^f}\hbox{\mathversion{bold}$\psi$}\|_{H^k(\Gamma^f)}$ in terms of the mesh size $h$, the parameter $k$,
and the singularity exponents.
For the norm of $\hbox{\mathversion{bold}$\psi$}$ this can be done component-wise using the same calculations
as in \cite{BespalovH_08_hpB} (see, e.g., inequalities (5.15) and (6.10) therein):
\begin{equation} \label{v_approx9}
    \|\hbox{\mathversion{bold}$\psi$}\|_{{\mathbf H}^k(\Gamma^f)} \le
    C\,h^{\alpha_2 + 1/2 - k}\,(\log(1/h))^{\tilde\beta_2+\tilde\nu_2},\qquad
    \alpha_2 + \hbox{$\frac 12$} \le k \le m,
\end{equation}
where $\alpha_2$ and $\tilde\beta_2$ are defined by (\ref{alpha_j}) and (\ref{beta_j}),
respectively, $\tilde\nu_2 = \frac 12$ if $k = \alpha_2 + \frac 12$ and
$\tilde\nu_2 = 0$ otherwise.

To estimate the norm $\|\div_{\Gamma^f}\hbox{\mathversion{bold}$\psi$}\|_{H^k(\Gamma^f)}$ we observe that the operator
$\div_{\Gamma^f}$ reduces all singularity exponents by one (while preserving the structure
of the corresponding singularity). Then, using similar calculations as indicated above,
we obtain
\begin{equation} \label{v_approx10}
    \|\div_{\Gamma^f}\,\hbox{\mathversion{bold}$\psi$}\|_{H^k(\Gamma^f)} \le
    C\,h^{\alpha_2 - 1/2 - k}\,(\log(1/h))^{\tilde\beta_2+\bar\nu_2},\qquad
    \alpha_2 - \hbox{$\frac 12$} \le k \le m,
\end{equation}
where $\alpha_2,\,\tilde\beta_2$ are the same as in (\ref{v_approx9}), whereas
$\bar\nu_2 = \frac 12$ if $k = \alpha_2 - \frac 12$ and
$\bar\nu_2 = 0$ otherwise.

By (\ref{v_approx8})-(\ref{v_approx10}) we conclude that
\begin{equation} \label{v_approx11}
    \|\hbox{\mathversion{bold}$\psi$} - \hbox{\mathversion{bold}$\psi$}^{hp}\|_{{\mathbf X}} \le
    C\,h^{\min\{k,p\}+\alpha_2-k}\,p^{-(k+1/2)}\,(\log(1/h))^{\tilde\beta_2+\bar\nu_2},\qquad
    \alpha_2 - \hbox{$\frac 12$} \le k \le m
\end{equation}
with the same $\alpha_2$, $\tilde\beta_2$, and $\bar\nu_2$ as in (\ref{v_approx10}).

Let $p > 2\alpha_2 - \frac 12$. Since $m$ is large enough, we can select an integer $k$ satisfying
\[
  2\alpha_2 - \hbox{$\frac 12$} < k \le \min\,\{m,p\}.
\]
Then $\min\,\{k,p\} = k$, and $p^{-(k+1/2)} \le p^{-2\alpha_2}$.

If $\alpha_2 - \frac 12 < p \le 2\alpha_2 - \frac 12$ (i.e., $p$ is bounded),
we choose an integer $k \in (\alpha_2 - \frac 12,p]$,
and if $p = \alpha_2 - \frac 12$, then we take $k = p = \alpha_2 - \frac 12$.
In both these cases $\min\,\{k,p\} = k$, and $p^{-(k+1/2)} \le C\,p^{-2\alpha_2}$.

Thus, for any $p \ge \alpha_2 - \frac 12$, selecting $k$ as indicated above
we find by (\ref{v_approx11})
\begin{equation} \label{v_approx12}
    \|\hbox{\mathversion{bold}$\psi$} - \hbox{\mathversion{bold}$\psi$}^{hp}\|_{{\mathbf X}} \le
    C\,h^{\alpha_2}\,p^{-2\alpha_2}\,(\log(1/h))^{\tilde\beta_2+\nu_2}
\end{equation}
with $\alpha_2$, $\tilde\beta_2$, and $\nu_2$ being defined by
(\ref{alpha_j}), (\ref{beta_j}), and (\ref{nu_j}), respectively.

{\bf Step 5: approximation of ${\mathbf v} = \hbox{\mathversion{bold}$\varphi$} + \hbox{\mathversion{bold}$\psi$}$.}
Let us define ${\mathbf v}^{hp} := \hbox{\mathversion{bold}$\varphi$}^{hp} + \hbox{\mathversion{bold}$\psi$}^{hp} \in {\mathbf X}_{hp}$, where
$\hbox{\mathversion{bold}$\varphi$}^{hp}$ and $\hbox{\mathversion{bold}$\psi$}^{hp}$ are approximations constructed above (see Steps~2--4).
Then combining estimates (\ref{v_approx7}) and (\ref{v_approx12}) and
using the triangle inequality we prove (\ref{approx_v}) in the case
$p \ge \alpha_2 - \frac 12$.

It remains to consider the case $1 \le p < \alpha_2 - \frac 12$.
In this case one does not need decomposition (\ref{v_approx3}).
Observe that for each face $\Gamma^f$ one has $\div_{\Gamma^f}\,{\mathbf v} \in H^k(\Gamma^f)$
with $1 \le k < \alpha -\frac 12$. Therefore, ${\mathbf v} \in {\mathbf X}^k$
with $1 \le k < \alpha -\frac 12$, and applying Proposition~\ref{prop_approx_reg}
we find ${\mathbf v}^{hp} \in {\mathbf X}_{hp}$ satisfying
\[
  \|{\mathbf v} - {\mathbf v}^{hp}\|_{{\mathbf X}} \le
  C\,h^{\min\,\{k,p\}+1/2}\,\|{\mathbf v}\|_{{\mathbf X}^k}.
\]
Hence, selecting $k \in [p, \alpha - \frac 12)$ we arrive at the desired
upper bound in (\ref{approx_v}), and the proof is finished.