text:
\section{Introduction}

A particular challenge for electric machine simulation is the relative motion of stator and rotor and the computation of quantities of interest depending on the rotation angle, e.g., the torque as a measure for the magneto-mechanic energy conversion.
Various approaches have been proposed to tackle the coupling of subdomain problems on moving geometries, e.g., moving band,  mortar and Lagrange multiplier or discontinuous Galerkin methods \cite{Alotto01,Buffa01,Lange10,Tsukerman95}, all with their specific advantages and shortcomings.
In addition, different strategies have been discussed for the numerical computation of the torque, e.g., via virtual displacements or the Maxwell stress tensor; see \cite{HenrotteHameyer04} for an overview.

In this paper, we consider harmonic mortar finite element and isogeometric analysis methods proposed in \cite{Bontinck18,DeGersem04}, which are based on finite element or isogeometric analysis approximations of the magnetic problems in the stator and rotor subdomains, coupled by a Lagrange multiplier technique using trigonometric functions.
These methods have several interesting properties:
\begin{enumerate}[(a)]
\item They are based on the Galerkin approximation of a certain weak formulation of the problem using Lagrange multipliers. This will allow us to give an unambiguous definition of the torque resulting in an exact energy conservation principle also on the discrete level.
\item The computed torque, as a function of the rotation angle, can be shown to be smooth, i.e.,  infinitely differentiable, and therefore no spurious torque ripples are introduced by non-smoothness of the numerical approximation.
\item The use of trigonometric functions for the Lagrange multipliers addition-ally allows for an efficient assembling of the coupling matrices for multiple rotation angles, avoiding the repeated integration of interface terms.
\end{enumerate}
The resulting methods therefore are particularly well-suited for multiquery applications, e.g., the construction of performance maps, the setup of reduced order models, coupled frequency or time domain simulations. This will be illustrated by computation of cogging torque in the numerical tests.

The remainder of the manuscript is organized as follows: Section \ref{sec:problem} introduces the magnetostatic setting for electric machines in two dimensions. The strong and weak formulations of the problem are presented and the energy balance is used to derive the torque. The discretization of the problem is given in Section \ref{sec:galerkin} using the Galerkin approach. Section \ref{sec:harmonicmortar} discusses the harmonic mortar method and its implementation. Numerical results for the torque computation using the harmonic mortar method are demonstrated in Section \ref{sec:results}. Finally the paper is concluded with a discussion in Section \ref{sec:conclusions}.

\section{Problem statement}\label{sec:problem}

For ease of presentation, we consider a linear magnetostatic setting in two space dimensions.
The geometric setup represents the cross section of a cylindrical device and consists of a stator domain $\Omega_S$ completely enclosing the rotor domain $\Omega_R$ with common interface $\Gamma = \partial\Omega_S \cap \partial\Omega_R = \partial\Omega_R$; see Figure~\ref{fig:geometry}.
\begin{figure}
		\centering
\includegraphics{images/harmonicmortar-figure0}
		\caption{Sketch of the geometry. $x_{\theta}$ denotes a point in the fixed coordinate system (red) of the stator domain $\Omega_S$, whereas $\tilde{x}_{\theta}$ denotes a point in the rotated coordinate system (blue) of the rotor domian $\Omega_R$. The are related by $x_\theta=\tilde{x}_{\theta-\alpha}$ where $\alpha$ is the angle of rotation. }
		\label{fig:geometry}
	\end{figure}
Motivated by the application context, we assume that $\Gamma$ is a circle placed within the air gap and centered around the origin.
Using a standard vector potential formulation \cite{Alotto01,DeGersem04}, the magnetic fields in the two subdomains are described by
\begin{alignat}{5}
    -\div(\nu_S \nabla a_S) &= j_e \qquad && \text{in } \Omega_S \label{eq:sys1}\\
    -\div(\nu_R \nabla a_R) &= \div M^\perp \qquad && \text{in } \Omega_R \label{eq:sys2}
\end{alignat}
where $a_S$, $a_R$ are the $z$-component of the magnetic vector potential, $j_e$ is the impressed electric source density, and
$M^{\perp} = (m_{y},-m_{x})$ the rotated magnetization vector.
The equations are formulated in two different coordinate systems attached to the respective subdomains, such that the position of material continuities stays fixed under relative motion.
We further assume that
\begin{alignat}{5}
a_S &= 0 && \qquad \text{on } \Sigma \label{eq:sys3}
\end{alignat}
on the outer boundary $ \Sigma = \partial\Omega_S \setminus \Gamma$.
The continuity of the magnetic field and the magnetic vector potential across the interface are encoded in
\begin{alignat}{5}
a_S &= a_R \circ \rho_{-\alpha} \qquad && \text{on } \Gamma, \label{eq:sys4}\\
n_\Gamma \cdot (\nu_S \nabla a_S)  &= n_\Gamma \cdot (\nu_R \nabla a_R) \circ \rho_{-\alpha} \qquad && \text{on } \Gamma, \label{eq:sys5}
\end{alignat}
where $\rho_{\alpha}$ describes the rotation of a point by angle $\alpha$ around the origin.
Furthermore, $n_\Gamma$ denotes the unit normal vector on $\Gamma$ pointing from $\Omega_R$ to $\Omega_S$.
The second coupling condition allows us to introduce
\begin{alignat}{5}
\lambda :=  n_\Gamma \cdot (\nu_S \nabla a_S) \label{eq:lambda}
\end{alignat}
which amounts to the tangential trace of the magnetic field at the interface.

\medskip

\textbf{Dependence on $\alpha$.}
By \eqref{eq:sys4}--\eqref{eq:sys5} the solutions $(a_S,a_R,\lambda)$ of \eqref{eq:sys1}--\eqref{eq:lambda} implicitly depend on the rotation angle $\alpha$, even if the sources $j_e$ and $M$ are independent of $\alpha$, which we assume in the sequel.
We write $u(\cdot|\alpha)$ to highlight the dependence of a function $u$ on the angle $\alpha$ and denote by
\begin{align}
u'(\cdot|\alpha) := \tfrac{d}{d\alpha} u(\cdot|\alpha) := \lim_{\varepsilon \to 0} \tfrac{1}{\varepsilon}(u(\cdot|\alpha+\varepsilon) - u(\cdot|\alpha))
\end{align}
the derivative of such a parameter-dependent function $u$ with respect to the parameter $\alpha$.
By differentiation of \eqref{eq:sys4} with respect to $\alpha$, we can see that
\begin{alignat}{5}
a'_S(x_\theta|\alpha) = a'_R(x_{\theta-\alpha}|\alpha) +  \tfrac{d}{d\theta} a_R(x_{\theta-\alpha}|\alpha) \qquad x_\theta \in \Gamma,
\end{alignat}
where $x_{\theta-\alpha} = \rho_{-\alpha}(x_\theta)$ is the point in the rotated coordinate system of $\partial\Omega_R$ corresponding to a point $x_\theta=(r_\Gamma \cos \theta,r_\Gamma \sin \theta)$ on the fixed interface $\Gamma$; see Figure~\ref{fig:geometry}.
Furthermore, $\frac{d}{d\theta} u(x_\theta) = r_\Gamma n_\Gamma^\perp \cdot \nabla u(x_\theta)$ amounts to the tangential derivative of $u$ along $\Gamma$ and $r_\Gamma$ is the radius of $\Gamma$.

\medskip

\textbf{Weak formulation.}
Using standard arguments, sufficiently regular solutions of \eqref{eq:sys1}--\eqref{eq:sys5} can be seen to satisfy the variational identities
\begin{alignat}{5}
(\nu_S \nabla a_S, \nabla v_S)_{\Omega_S} + \langle \lambda, v_S\rangle_\Gamma &= (j_e, v_S)_{\Omega_S} \label{eq:var1}\\
(\nu_R \nabla a_R, \nabla v_R)_{\Omega_R} - \langle \lambda, v_R \circ \rho_{-\alpha}\rangle_\Gamma &= (-M^\perp, \nabla v_R)_{\Omega_R} \label{eq:var2}\\
\langle a_S - a_R \circ \rho_{-\alpha}, \mu\rangle_\Gamma &= 0  \label{eq:var3}
\end{alignat}
for all smooth test functions $v_S$, $v_R$, and $\mu$ defined on $\Omega_S$, $\Omega_R$, and $\Gamma$, respectively, with $v_S=0$ on $\Sigma$.
Using the considerations outlined above, the coupling condition \eqref{eq:var3} can also be rephrased as
\begin{alignat}{5}
\langle a'_S - a'_R \circ \rho_{-\alpha}, \mu\rangle_\Gamma &= \langle \tfrac{d}{d\theta} a_R \circ \rho_{-\alpha}, \mu\rangle_\Gamma \label{eq:var3a}
\end{alignat}
which will be used in the following considerations.

\medskip

\textbf{Energy balance and torque computation.}
As a next step, we introduce the magnetic energy of the system.
For ease of notation, we assume a unit length $L$  of the device in axial direction. Then
\begin{align} \label{eq:energy}
    E(a_S,a_R) = \int_{\Omega_S} \tfrac{\nu_S}{2} |\nabla a_S|^2 dx + \int_{\Omega_R} \tfrac{\nu_R}{2} |\nabla a_R|^2 + M^\perp \cdot \nabla a_R \, dx.
\end{align}
Let us note that, via the fields $a_S$, $a_R$, the energy also implicitly depends on the rotation angle $\alpha$. By elementary calculations, we may then compute
\begin{align*}
\tfrac{d}{d\alpha} E(a_S,a_R)
&= (\nu_S \nabla a_S,  \nabla a'_S)_{\Omega_S} + (\nu_R \nabla a_R, \nabla a'_R )_{\Omega_R} + (M^\perp, \nabla a'_R)_{\Omega_R} \\
&= (j_e, a'_S)_{\Omega_S} - \langle \lambda,  a'_S\rangle_\Gamma
+ \langle \lambda, a'_R \circ \rho_{-\alpha}\rangle_\Gamma,
\end{align*}
where we used the variational identities \eqref{eq:var1}--\eqref{eq:var2} with $v_S=a'_S$ and $v_R = a'_R$ in the second step. Using the differential form~\eqref{eq:var3a} of the coupling condition, we arrive at the energy balance
\begin{align*}
\tfrac{d}{d\alpha} E(a_S,a_R)
&= (j_e,a'_S)_{\Omega_S}- \langle \lambda, \tfrac{d}{d\theta} a_R \circ \rho_{-\alpha} \rangle_\Gamma,
\end{align*}
where the first term denotes the electric work required to maintain the electric current $j_e$ in the stator coils (neglecting Ohmic losses), and the second term is the mechanic work required for an infinitesimal rotation of the system.
These \emph{energy-based} considerations immediately give rise to the following definition of the torque
\begin{align} \label{eq:torque}
T(\alpha) := \langle \lambda, \frac{d}{d\theta} a_R \circ \rho_{-\alpha} \rangle_\Gamma.
\end{align}
Note that $a_R$ and $\lambda$ and hence $T$ are functions of the rotation angle~$\alpha$.
In general, the energy and torque have to be scaled by the axial length~$L$.

\begin{remark}
The above formula can be derived in various ways, see \cite{Bossavit90,Coulomb83,HenrotteHameyer04}. Since the stator and rotor were both assumed to be rigid bodies, the torque can be computed here directly without resorting to the Maxwell stress tensor.
Using rotational symmetry as well as integration by parts along the interface $\Gamma$, we can express the torque alternatively as
\begin{align}
T(\alpha) = -\langle \tfrac{d}{d\theta} \lambda \circ \rho_{\alpha}, a_R\rangle_\Gamma,
\end{align}
which may be a more convenient representation depending on the particular setting. \end{remark}

\section{Galerkin approximation}\label{sec:galerkin}

For discretization of the problem \eqref{eq:sys1}--\eqref{eq:sys5}, we consider Galerkin approximations of the variational equations \eqref{eq:var1}--\eqref{eq:var3} using appropriate finite-dimensional sub-spaces $V_{h,S} \subset H^1_{\Sigma}(\Omega_S)$, $V_{h,R} \subset H^1(\Omega_R)$, and $M_N \subset H^{-1/2}(\Gamma)$.
Let us note that the boundary conditions \eqref{eq:sys3} are incorporated explicitly in the definition of $H^1_\Sigma(\Omega_S)=\{v \in H^1(\Omega_S) : v=0 \text{ on } \Sigma\}$ here.
The discrete problem for a fixed angle $\alpha$ then reads as follows.
\begin{problem} \label{prob:galerkin}
Find $a_{h,S} \in V_{h,S}$, $a_{h,R} \in V_{h,R}$ and $\lambda_N \in M_N$ such that
\begin{alignat}{5}
(\nu_S \nabla a_{h,S}, \nabla v_{h,S})_{\Omega_S} + \langle \lambda_N, v_{h,S}\rangle_\Gamma &= (j_e, v_{h,S})_{\Omega_S} \label{eq:var1h}\\
(\nu_R \nabla a_{h,R}, \nabla v_{h,R})_{\Omega_R} - \langle \lambda_N, v_{h,R} \circ \rho_{-\alpha}\rangle_\Gamma &= -(M^\perp, \nabla v_{h,R})_{\Omega_R} \label{eq:var2h}\\
\langle \alpha_{h,S} - \alpha_{h,R} \circ \rho_{-\alpha}, \mu_N\rangle_\Gamma &= 0  \label{eq:var3h}
\end{alignat}
for all discrete test functions $v_{h,S} \in V_{h,S}$, $v_{h,R} \in V_{h,R}$, and $\mu_N \in M_N$.
\end{problem}

\textbf{Well-posedness.}
For the following considerations, we assume that $V_{h,S}$, $V_{h,R}$, and $M_N$ are finite-dimensional and that $\nu_S$, $\nu_R$ are uniformly positive and bounded from above. Moreover,the source currents $j_e$ and magnetization vector $M$, as well as the domains $\Omega_S$, $\Omega_R$ are assumed to be sufficiently regular. We can then establish the well-posedness of the discretization scheme as follows.

\begin{lemma} \label{lem:well-posed}
Assume that
\begin{align} \label{eq:infsup}
\sup_{v_{h,S} \in V_{h,S}} \frac{\langle v_{h,S}, \mu_N\rangle_\Gamma}{\|v_{h,S}\|_{H^1(\Omega_S)}} \ge \beta \|\mu_N\|_{H^{-1/2}(\Gamma)} \qquad \text{for all } \mu_N \in M_N.
\end{align}
Then the discrete variational problem \eqref{eq:var1h}--\eqref{eq:var3h} is uniquely solvable and
\begin{align*}
    &\|a_S - a_{h,S}\|_{H^1(\Omega_S)} + \|a_R - a_{h,R}\|_{H^1(\Omega_R)} + \|\lambda - \lambda_N\|_{H^{-1/2}(\Gamma)} \\
    &\le C(\beta) \left( \|a_S - v_{h,S}\|_{H^1(\Omega_S)} + \|a_R - v_{h,R}\|_{H^1(\Omega_R)} + \|\lambda - \mu_N\|_{H^{-1/2}(\Gamma)}\right)
\end{align*}
for all $v_{h,S} \in V_{h,S}$, $v_{h,R} \in V_{h,R}$ and $\mu_N \in M_N$.
\end{lemma}
The assertion is a direct consequence of Brezzi's saddle-point theory \cite{Brezzi74}; also see \cite{Buffa01,Egger_2021aa} for application in the current context.
Under the stability condition \eqref{eq:infsup}, the Galerkin approximations obtained by Problem~\ref{prob:galerkin} thus always lead to quasi-optimal error estimates with respect to the chosen approximation spaces. Two particular examples of Galerkin approximations will be discussed below.

\medskip

\textbf{Discrete torque computation and energy balance.}
Mimicking the formula \eqref{eq:torque} for the torque obtained on the continuous level, it seems natural to define the discretized approximation of the torque as
\begin{align} \label{eq:torqueh}
    T_{h,N}(\alpha) :=
    \langle \lambda_h, \tfrac{d}{d\theta} a_{h,R} \circ \rho_{-\alpha}\rangle_\Gamma.
\end{align}
The torque computation for a given angle $\alpha$ thus only requires the solution of one magnetostatic interface problem \eqref{eq:var1h}--\eqref{eq:var3h} for this specific angle.
With the same reasoning as on the continuous level, we obtain the following result.

\begin{lemma} \label{lem:energy_discrete}
Let $(a_{h,S},a_{h,R},\lambda_N)$ denote a solution of \eqref{eq:var1h}--\eqref{eq:var3h}. Then
\begin{align}
    \tfrac{d}{d\alpha} E(\alpha_{h,S},\alpha_{h,R}) = (j_e,a'_{h,S})_{\Omega_S} - T_{h,N}
\end{align}
with discrete torque $T_{h,N}$ defined by the formula \eqref{eq:torqueh}.
\end{lemma}
This assertion illustrates that the electro-magneto-mechanic energy balance holds exactly also on the discrete level.

\medskip

\textbf{Alternative representation and smooth dependence on $\alpha$.}
If the space $M_N$ of discrete Lagrange multipliers consists of sufficiently smooth functions, we may express the discrete torque equivalently as
\begin{align} \label{eq:torqueh2}
    T_{h,N} =
-\langle \tfrac{d}{d\theta} \lambda_N \circ \rho_{\alpha}, a_{h,R}\rangle_\Gamma.
\end{align}
By a recursive argument, one can then show the following result.
\begin{lemma} \label{lem:smooth}
Let $M_N$ consist of $k-1$ times continuously differentiable functions.
Then the solution $(a_S,a_R,\lambda_N)$
is $k$ times differentiable with respect to $\alpha$.
\end{lemma}
A proof of this assertion follows from the algebraic structure of the discretized problem; see the end of the next section.

\section{Harmonic mortar discretizations.} \label{sec:harmonicmortar}

As a particular choice of a Galerkin approximation, we now consider the following choice of spaces: $V_{h,S}$ and $V_{h,R}$ are standard $H^1$-conforming finite element spaces \cite{Braess} over appropriate triangulations of the stator and rotor domain $\Omega_S$ and $\Omega_R$. Note that no conformity of the mesh across the interface is required. For approximation of the Lagrange multipliers, we consider the space
\begin{align} \label{eq:M}
M_N = \left\{ \mu_N(x_\theta) = \frac{c_0}{2} + \sum\nolimits_{n=1}^N c_n \cos(n \theta) + d_n \sin(n \theta) : c_n, d_n \in\mathbb{R}\right\}
\end{align}
of trigonometric polynomials of degree $\le N$.

\begin{remark}
As shown in \cite{Egger_2021aa}, the inf-sup condition \eqref{eq:infsup} holds true, if $N$ is chosen appropriately, in particular not too large. By the results of the previous section, the harmonic mortar finite element method then is stable and yields quasi-optimal error estimates. The same conclusions hold for a method based on isogeometric analysis discretization in the stator and rotor domains \cite{Bontinck18,Cottrell_2009aa}.
Following the considerations of the previous section, both methods allow for a torque computation consistent with a discrete energy conservation principle.
\end{remark}

\textbf{Implementation of the harmonic mortar finite element method.}
In the following, we discuss some details concerning the implementation of the harmonic mortar finite element method outlined above.
As we will see, the use of trigonometric polynomials for the discrete Lagrange multipliers brings advantages also for the numerical realization.

After choosing appropriate basis functions for the finite element subspaces $V_{h,S}$, $V_{h,N}$ and the trigonometric polynomials in $M_N$,
the discrete variational problem \eqref{eq:var1h}--\eqref{eq:var3h} can be turned into a linear system
\begin{alignat}{5}
    K_S a_S(\alpha) & &    &+B_S^\top \lambda(\alpha) &&= j_e \label{eq:lin1}\\
    K_R a_R(\alpha) & &    &-B_R(\alpha)^\top \lambda(\alpha) \ &&= j_M \label{eq:lin2}\\
    B_S a_S(\alpha) &-& B_R(\alpha) a_R(\alpha)& &&= 0. \label{eq:lin3}
\end{alignat}
The source vectors $j_e$ and $j_M$ and the matrices $K_S$, $K_R$ and $B_S$ are independent of the angle $\alpha$, whereas $B_R$ and hence also the solution vectors $a_S$, $a_R$ and $\lambda$ depend on the rotation angle.
For the harmonic mortar, the entries of the coupling matrix $B_R(\alpha)$ can be computed by the following expressions
\begin{align*}
[B_R(\alpha)]_{2n,j} &= \int_0^{2\pi} \cos(n \theta) \, \phi_j(x_{\theta-\alpha}) r_\Gamma d\theta, \qquad n \ge 0,  \\
[B_R(\alpha)]_{2n-1,j} &=  \int_0^{2\pi} \sin(n \theta) \, \phi_j(x_{\theta-\alpha}) r_\Gamma d\theta, \qquad n \ge 1,
\end{align*}
where $\phi_j$ are the basis functions for the Galerkin approximation on $\Omega_R$.

\begin{remark}
By trigonometric summation formulas, we get
\begin{align*}
    [B_R(\alpha)]_{2n,j} &= \cos(n\alpha)[B_R(0)]_{2n,j} -\sin(n\alpha)[B_R(0)]_{2n-1,j}  \\
    [B_R(\alpha)]_{2n-1,j} &= \sin(n\alpha)[B_R(0)]_{2n,j} + \cos(n\alpha)[B_R(0)]_{2n-1,j},
\end{align*}
with smooth coefficients depending on $\alpha$, hence the integrals in the coupling matrix $B_R(\alpha)$ in fact only have to be evaluated for a single angle $\alpha=0$.
In matrix notation, we may then write
\begin{align} \label{eq:BRa}
    B_R(\alpha) = R(\alpha) B_R(0)
\end{align}
where $R(\alpha)$ is a block diagonal matrix with $2\times 2$ blocks.
Further, note that the integrals in the definition of $B_R(0)$ consist of products of trigonometric functions with polynomials, which can be computed analytically \cite{Neuman_1981aa}.
\end{remark}

\textbf{Algebraic energy balance and torque computation.}
On the algebraic level, the energy is given by $E=\frac{1}{2} \|a_S\|^2_{K_S} + \frac{1}{2} \|a_R\|_{K_R}^2 - j_M^\top a_R$ where $\|x\|_M^2 = x^\top M x$ is the square of the norm associated with a symmetric positive semi-definite matrix $M$. Then
\begin{align*}
    \tfrac{d}{d\alpha} E
    &= (K_S a_S, a_S') + (K_R a_R, a_R') - (j_M,a_R') \\
    &= -(B_S^\top \lambda, a_S') + (j_e, a_S') + (B_R^\top \lambda, a_R') \\
    &= (j_e, a_S') - (B_S a_S' - B_R a_R', \lambda)
    =(j_e, a_S') - (B_R' a_R, \lambda).
\end{align*}
Hence the discrete torque can be simply computed by the algebraic formula
\begin{align*}
    T_{h,N}(\alpha) = \lambda^\top(\alpha) B_R'(\alpha) a_R(\alpha) =\langle \lambda_h, \tfrac{d}{d\theta}a_{h,R}\circ\rho_{-\alpha}\rangle .
\end{align*}
Note that $B_R'(\alpha)$ can be expressed as $B_R'(\alpha)=R'(\alpha)B_R(0)$ by using \eqref{eq:BRa}. Moreover, the matrix $R'(\alpha)$ is block diagonal with $2 \times 2$ blocks whose entries can be computed analytically; see above.

\medskip

\textbf{Derivatives with respect to $\alpha$ and proof of Lemma~\ref{lem:smooth}.}
By formal differentiation of the above equations, we can see that
\begin{alignat}{5}
    K_S a'_S(\alpha) - B_S^\top \lambda'(\alpha) &= 0 \label{eq:lin1d}\\
    K_R a'_R(\alpha) - B_R(\alpha)^\top \lambda'(\alpha) &= B_R'(\alpha)^\top \lambda(\alpha) \label{eq:lin2d}\\
    B_S a'_S(\alpha) - B_R(\alpha) a'_R(\alpha) &=  B_R'(\alpha) a_R(\alpha). \label{eq:lin3d}
\end{alignat}
Together with \eqref{eq:lin1d}--\eqref{eq:lin3d}, one can see that the solution of \eqref{eq:lin1}--\eqref{eq:lin3} is differentiable with respect to  the angle $\alpha$.
By recursion, we may further conclude that the torque $T_{h,N}$ of the harmonic mortar method is a smooth function of the angle $\alpha$. The same argument provides a proof of Lemma~\ref{lem:smooth}.

\section{Numerical results} \label{sec:results}

In the following, we demonstrate the application of the harmonic mortar finite element and isogeometric analysis method for computing the torque in a six-pole permanent magnet synchronous machine shown in Fig.~\ref{fig:pmsm}.
\begin{figure}[ht!]
 \centering
 \includegraphics[width=0.38\linewidth]{images/PMSM}
 \includegraphics[width=0.61\linewidth]{images/full_machine_isolines_small}
 \caption{\textit{Left:} Model of a six-pole permanent magnet synchronous machine. The rotor on the inside consists of iron (green) and permanent magnets (blue). The stator on the outside consists of iron (green) and the slots for the copper windings (light green). \textit{Right:} Magnetic vector potential $a$ and magnetic field lines in the machine.}
 \label{fig:pmsm}
\end{figure}
The geometry consists of a rotor with three pole pairs with embedded permanent magnets and a stator with $36$ slots filled with copper windings.
The length of the actual machine is $L=\SI{0.1}{m}$, so the energy and torque have to be scaled by $L$ as outlined above.
Some information about the geometry and material parameters is summarized in Table~\ref{table:machine_parameters};
for details see \cite[Chapter V.A]{Bontinck_2018af}.
\begin{table}[ht!]
\centering
\caption{Geometry and material parameters of the machine.}
\label{table:machine_parameters}
\begin{tabular}{ll}
\hline
\multicolumn{2}{c}{geometry}\\
\hline
inner radius rotor & \SI{16}{mm}\\
outer radius rotor & \SI{44}{mm}\\
inner radius stator & \SI{45}{mm}\\
outer radius stator & \SI{67.5}{mm}\\
\hline
\multicolumn{2}{c}{materials}\\
\hline
relative permeability iron & 500\\
relative permeability copper & 1\\
relative permeability permanent magnets & 1.05\\
remanence field permanent magnet & \SI{0.94}{T}\\
\end{tabular}
\end{table}

\medskip

\textbf{Simulation setup.}
In our numerical tests, we compute the cogging torque via formula \eqref{eq:torqueh2} in the absence of excitation currents $j_e=0$.
For the discretization of the magnetostatic problems in the stator and rotor domains, we used open source package \textsc{GeoPDEs} \cite{Vazquez_2016aa}. The computations are carried out with B-splines of degree $p=1$ and $p=2$ with maximal continuity as basis functions for the solution, which corresponds to a standard finite element method (FEM) and a multipatch spline discretization with patchwise $C^1$-continuous splines of second order (IGA). To accommodate for material discontinuities, only $C^0$-continuity is enforced across patch borders; see \cite{Bontinck18} for details.
At the external boundaries of the stator and rotor, the magnetic vector potential $a$ is set to zero.
For the computations of different degree $p$ the mesh refinement has been adapted such that the number of degrees of freedom is comparable for $p=1$ and $p=2$. The number of degrees of freedom for $p=1$ is $N_{\mathrm{dof},R}=26532$ in the rotor domain and $N_{\mathrm{dof},S}=39528$ in the stator domain. For $p=2$ the number of degrees of freedom is $N_{\mathrm{dof},R}=25140$ in the rotor domain and $N_{\mathrm{dof},S}=37764$ in the stator domain.
The reference solution is computed on a  refined mesh with $p=2$ leading to $N_{\mathrm{dof},R}=93360$ and $N_{\mathrm{dof},S}=128664$ degrees of freedom in the rotor and stator domain, and with trigonometric degree $N=200$.
The saddlepoint problems \eqref{eq:lin1}--\eqref{eq:lin3} are solved efficiently by a Schur complement method;  details are given in Section~\ref{sec:appendix}.

\medskip

\textbf{Symmetry of the solution.}
A snapshot of the magnetic vector potential and its isolines, which correspond to the magnetic flux lines, is depicted in Fig.~\ref{fig:pmsm}.
Note that by symmetry of the rotor and stator geometry, the solution components $(a_S,a_R,\lambda)$ have a $60^\circ$ anti-symmetry. In particular, for the Fourier expansion of the computed Lagrange multiplier
\begin{align} \label{eq:lagrange}
    \lambda_N(\theta) = \frac{c_0}{2} + \sum_{n=1}^N c_n \cos (n \theta) + d_n \sin(n \theta),
\end{align}
only the coefficients $c_n$, $d_n$ with index $n\in \mathcal{N}=\{3,9,15,...\}$
are non-zero.
By careful construction of the meshes, the geometric symmetries are preserved on the discrete level and the same behavior is observed for the numerical solution, as can be seen from Table~\ref{tab:lagrange}.
\begin{table}[ht!]
\centering
\caption{Relevant and irrelevant Fourier modes of the Lagrange multipliers of an IGA simulation result for $\alpha=\SI{7}{\degree}$, $C_n=\sqrt{c_n^2+d_n^2}$ is the total amplitude.\label{tab:lagrange}}
\begin{tabular}{c|c|c|c}
$C_3$ & $C_9$ & $C_{15}$ & $\sum\nolimits_{n\not\in\mathcal{N}}C_n$ \\
\hline
$\SI{9.2621e+03}{}$ & $\SI{1.6651e+03}{}$ &  $\SI{5.5848e+03}{}$ & $\SI{4.9720e-06}{}$
\end{tabular}
\end{table}
Using this prior knowledge, most of the Lagrange multipliers can be eliminated already during assembling leading to more efficient computations, in particular in connection with the Schur complement technique outlined in Section~\ref{sec:appendix}.

\medskip

 \textbf{Torque computation.}
In Fig.~\ref{fig:torque}, we display the torque computed via formula \eqref{eq:torque} with FEM ($p=1$) and IGA ($p=2$) for different orders $N$ of harmonic basis functions.
\begin{figure}[ht!]
\centering
\includegraphics{images/harmonicmortar-figure1}
\caption{Cogging torque of the machine. Computed using trigonometric polynomials of degree $N$ as basis for the Lagrange multiplier space. Computed with FEM and IGA using basis functions of degree $p=1,2$.}
 \label{fig:torque}
 \end{figure}
By the $10^\circ$ symmetry of the stator slots and the corresponding symmetry of the mesh, the computed torque can be shown to exhibit a periodic behavior with a period length of $10^\circ$.
Moreover, the torque can be seen to be anti-symmetric with respect to half of the period. We further observe that the torque converges towards that of the reference solution when increasing the number $N$ of harmonic basis functions. As explained by the theoretical results in \cite{Egger_2021aa}, the computation becomes unstable when increasing the dimension of the Lagrange multiplier spaces too much; here for $N=250$.

Due to its periodicity, the torque $T_{h,N}$ can be approximated by a Fourier series
\begin{align} \label{eq:torque_fourier}
T_{h,N}^M(\alpha) = \frac{\hat{c}_0}{2} + \sum\limits_{m=1}^{M} \hat{c}_m \cos(m \alpha) + \hat{d}_m \sin(m \alpha)
\end{align}
with coefficients $\hat{c}_m, \hat{d}_m \in\mathbb{R}$.
By the symmetry properties discussed above, one can see that only $\hat{d}_m$ for $m\in \mathcal{M}=\{36,72,108,...\}$ is non-zero, while $\hat{d}_m=0$ of all $m \not\in \mathcal{M}$ and $\hat{c}_m=0$ for all $m$.
The corresponding results of our computations are shown in Table~\ref{table:torque_modes}.
\begin{table}[ht!]
\centering
\caption{Relevant and irrelevant Fourier modes of the torque of an IGA simulation result for a full rotation of the machine.}
\label{table:torque_modes}
\begin{tabular}{c|c|c|c|c}
$|\hat{d}_{36}|$ & $|\hat{d}_{72}|$ & $|\hat{d}_{108}|$ & $\sum\limits_{m}|\hat{c}_m|$ & $\sum\limits_{m\not\in\mathcal{M}}|\hat{d}_m|$ \\
\hline
$\SI{2.293e-01}{}$ & $\SI{1.784e-01}{}$ &  $\SI{9.17e-02}{}$ & $\SI{6.2079e-11}{}$ & $\SI{51.7843e-11}{}$
\end{tabular}
\end{table}

\section{Discussion} \label{sec:conclusions}

We have discussed the harmonic mortar method using trigonometric polynomials as Lagrange multipliers for the efficient treatment of rotation in electric machine computation. The energy balance has been used to derive an explicit formulation for the computation of the torque in the continuous and discretized settings. It has been shown that the electro-magneto-mechanic energy balance holds exactly also on the discrete level.
Furthermore, if the functions in the Lagrange multiplier space are smooth then this is inherited by the torque.
The derived formulation has been applied to a typical permanent magnet synchronous machine model using isogeometric analysis and a lowest order finite element method with an exact geometry mapping. The cogging torque has been computed with the proposed method and the influence of the number of Lagrange multiplier basis functions on the solution has been investigated.
\section*{Acknowledgements}
This work is supported by the Graduate School CE within the Centre for Computational Engineering at Technische Universität Darmstadt, by the German Research Foundation via the project SCHO 1562/6-1, and by the Defense Advanced Research Projects Agency (DARPA), under contract HR0011-17-2-0028. The views, opinions and/or findings expressed are those of the author and should not be interpreted as representing the official views or policies of the Department of Defense or the U.S. Government.


text:
\section{Polyominoes.}
\emph{Polyominoes}, first studied systematically by Golomb \cite{golomb1954checker}, are shapes that can be made by gluing together finitely many unit squares, edge to edge. A polyomino with $n$ squares is sometimes called an \emph{$n$-omino}.

Tiling problems involving polyominoes are well studied---see, for example, \cite{MR1291821} or \cite{boltjanskij1991geometric}. In tiling problems, one almost always restricts to simply-connected polyominoes, i.e.,\ polyominoes without holes. Our main interest here is maximizing the number of holes. Figure \ref{Fonehole} illustrates an $8$-omino, a $20$-omino, and a $60$-omino with $1$, $5$, and $21$ holes respectively.

To be precise about the topology, we consider the tiles to be closed unit squares in the plane. Polyominoes are finite unions of these closed squares, so they are compact. The holes of a polyomino are defined to be the bounded connected components of its complement in the plane.

\input{onehole.txt}

\subsection{Counting Holes.}
We consider two $n$-ominoes to be equivalent if they agree after a translation. We denote the set of all $n$-ominoes by $\mathcal{A}_{n}$.
Given a polyomino $A$, we denote by $h(A)$ the number of holes in $A$. For $n\ge 1$ define,
\begin{equation}
f(n):= \max_{A\in \mathcal{A}_n}h(A).
\end{equation}

Other notions of polyomino equivalence have been studied, especially in the context of polyomino enumeration. For example, one might consider two polyominoes to be the same if they agree after rotation or reflection. However, the number of holes is invariant under rigid motions, so for our purposes, all that matters is that $\mathcal{A}_n$ is a finite set so the maximum $f(n)$ is well defined. This maximum is the same under any of these notions of equivalence.

Similarly, define $g(m)$ to be the minimum number $N$ such that there exists an $N$-omino with $m$ holes. We check below that $g$ is a right inverse of $f$.

The function $g$ is listed at The On-Line Encyclopedia Of Integer Sequence as sequence A118797.
Tom\'as Oliveira e Silva enumerated free polyominoes according  to area and number of holes, up to $n=28$, which at the time of this writing seems to be the state of the art. As a corollary of his calculations, we know $f(n)$ for $n \le 28$ and $g(m)$ for $m \le 8$. \cite{bworld}. See Table 1.

\begin{table}
\centering
\begin{tabular}{l | l | l | l | l | l | l | l | l}
m & 1 & 2  & 3  & 4  & 5  & 6  & 7  & 8  \\ \hline
g(m) & 7 & 11 & 14 & 17 & 19 & 23 & 25 & 28 \\
\end{tabular}
\caption{$g(m)$ for $m \le 8$.}
\end{table}

\subsection{Statement of main results.}

Here and throughout,
\begin{equation}\label{Defn}
n_k=\frac{1}{3} \left( 2^{2k+1}+3 \cdot 2^{k+1}+4 \right),
\end{equation}
and
\begin{equation}\label{Defh}
h_k=\frac{1}{3} \left(2^{2k}-1 \right).
\end{equation}

\begin{theorem}\label{Ttremebundoexactly}
For $k \ge 1$ we have
$f(n_k)=h_k$. Moreover,
$f(n_k - 1)=h_k$ and $f(n_k -2) =h_k-1$.
\end{theorem}

\begin{corollary}\label{aboutexactg}
$g(h_k) = n_k-1$ for all $k \ge 1$.\\
\end{corollary}

We also give bounds for large $n$.

\begin{theorem}\label{Ttremebundo}
Let $f(n)$ denote the maximum number of holes that a polyomino with $n$ squares can have. Given $C_{1}>\sqrt{5 / 2}$ and $C_{2}<\sqrt{3 / 2}$, there exists an $n_0=n_0(C_1,C_2)$ such that
\begin{equation*}
\frac{1}{2}n-C_{1}\sqrt{n} \leq f(n) \leq \frac{1}{2}n-C_{2}\sqrt{n},
\end{equation*}
for $n>n_0$.
\end{theorem}

Theorem \ref{Ttremebundoexactly} implies that the constant $\sqrt{3/2}$ in Theorem \ref{Ttremebundo} can not be improved.

\section{Preliminary Results.}\label{Spreliminaryr}

It is clearly always possible to attach a square tile to an $n$-omino and obtain an $(n+1)$-omino with at least the same number of holes. This implies that $f(n)\leq f(n+1)$ for evey $n \ge 1$, so $f$ is monotonically increasing. The following lemma tells us that $f$ never increases by more than one.

\begin{lemma}\label{Ljumpsoff}
For every $n \ge 1$, $$ f(n+1)-f(n)\leq 1.$$
\end{lemma}
\begin{proof}
We will show that if $A$ is any $(n+1)$-omino, then there exists an $n$-omino $B$ such that $h(B) \ge h(A)-1$.

Let $A$ be an $(n+1)$-omino, and $k$ be the number of tiles in the bottom row of $A$, and denote by $l$ the leftmost tile in this row

If $k=1$, then $l$ is only connected to one other tile, so we can delete $l$ without disconnecting $A$ or destroying any holes.

Now suppose the statement holds whenever $k =1, 2, \dots,  m$. Then let $k=m+1$ and denote by $l_1$, $l_2$, and $l_3$ the tile sites that share boundary with $l$ that are to the up and right to $l$.
Each of the tile sites $l_1$, $l_2$, and $l_3$ could either be occupied by tiles in $A$ or not. The six possibilities are depicted in Figure \ref{Flcomb2}. Any other combination would result in $A$ having disconnected interior.

\input{lcomb_2.txt}

If $A$ is such that the local tile structure around $l$ coincides with $C_1$, $C_2$, $C_3$, $C_5$, or $C_6$, then it is already possible to delete $l$ from $A$ to generate an $n$-omino $B$ such that $h(B) \geq h(A)-1$.

However, if $C_4$ is the local structure around $l$, then it is possible that deleting $l$ disconnects $A$. In this case, we delete $l$ and then add a new tile at the empty tile site $l_2$. This yields a new polyomino $A^{\prime}$ with the same number of tiles. If the addition of this new tile causes the coverage of a hole, then $h(A^\prime)= h(A)-1$, and $C_1$ or $C_3$ must then be the new local structure around the leftmost tile of the bottom row of $A^\prime$. This then allows us to terminate the process by deleting the bottom leftmost tile from $A^\prime$ without destroying more holes. If we have not destroyed any holes, then we have an $(n+1)$-omino $A^\prime$ with $h(A)=h(A^\prime)$ and with $k=m$ tiles in the bottom row. Hence, we can apply the induction hypothesis and the desired result follows.
\end{proof}

\begin{lemma}\label{Lfandg}
For every $m \ge 1$, we have that $f(g(m))=m$. Also, $g(m) = n$ if and only if $f(n)=m$ and $f(n-1) = m-1$.
\end{lemma}
\begin{proof}
By the definitions of $f$ and $g$, we have immediately that $f \left( g(m) \right) \geq m $. Suppose by way of contradiction that for some $m$ we have $f(g(m)) \ge m+1$. By Lemma \ref{Ljumpsoff}, $f(g(m)-1)\geq m$. This implies that there exists a polyomino with $g(m)-1$ tiles and at least $m$ holes, but this contradicts the definition of $g$. We then conclude that $f ( g(m) ) = m $ for every $m$. It follows immediately that $g(m) = n$ if and only if $f(n)=m$ and $f(n-1) = m-1$.
\end{proof}

\section{An Upper Bound.}
\subsection{Perimeter.}
Define the perimeter $p(A)$ of a polyomino $A \in \mathcal{A}_n$ as the number of edges in the topological boundary of $A$. That is, $p(A)$ counts the number of edges with a square in $A$ on one side, and a square not in $A$ on the other side. For example, the 8-omino in Figure \ref{Fonehole} has a perimeter equal to 16.

For $n \ge 1$ we denote $p_{\min}(n)$ and $p_{\max}(n)$ as the minimum perimeter and the maximum perimeter possible for a polyomino with an area of $n$.

In 1976, F. Harary and H. Harborth \cite{harary1976extremal} proved that the minimum perimeter possible in an  $n$-omino is given by,
\begin{equation}\label{perimeter1976}
p_{min}(n)=2\ceil{2\sqrt{n}}.\\ \\
\end{equation}

\input{spiral.txt}

 Let $A \in \mathcal{A}_{n}$, then the number of edges that are on the boundary of two squares of $A$ will be denoted by $b(A)$,  in which case the edges are contained in the interior of the polyomino. Observe that all the edges of the squares of a polyomino either belong to the perimeter or to the interior of the polyomino. This means that,
 \begin{equation}\label{perimeterequal}
 4n= p(A)+ 2 b(A).
 \end{equation}

For example, if $A$ is any 7-omino depicted in Figure \ref{wormlike}, then $b(A)=n-1=6$ and $p(A)=4n-2(n-1)=16$.

Let $b_{min}(n)$ be the minimum number of edges shared by two squares that an $n$-omino can have. It is possible to associate a dual graph with any polyomino by considering each square as a vertex and by connecting any two of these vertices if they share an edge.

\input{wormlike.txt}

 Let $A$ be an $n$-omino. Because any polyomino has to have a connected interior, then the associated dual graph of $A$ is connected and there exists a spanning tree of this graph. If the dual graph has $n$ vertices, then there are at least $n-1$ edges in the spanning tree. This implies that there are at least $n-1$ different common edges in $A$. That is, $b(A)\geq n-1$, and this is true for any $n$-omino. This then gives us that $b_{\min}(n)\geq n-1$. Observing that the $n$-omino $C$ with only one column has $b(C)=n-1$, we conclude that $b_{min}(n)=n-1$.

As a consequence of this and equality (\ref{perimeterequal}), for every $n\geq 1$,
\begin{equation}
p_{max}(n)\leq 4n-2b_{min}(n)=4n-2(n-1)=2n+2,
\end{equation}
and equality is only achieved by polyominoes with the number of common edges equal to $b_{min}(n)$.

We need to distinguish when an edge that is on the perimeter of a polyomino $A \in \mathcal{A}_{n}$ is an edge that forms part of a hole in the polyomino. Define such edges as being part of the hole perimeter. We represent by $p_{h}(A)$ the number of edges on the hole perimeter of $A$. We define the outer perimeter of $A$, denoted by $p_{o}(A)$, as the difference between the perimeter $p(A)$ and the hole perimeter,
\[p_{o}(A)=p(A)-p_{h}(A).\]

If a polyomino $A$ is simply connected,
 then $p(A)=p_{o}(A)$. In general, $p(A)=p_{o}(A)+p_{h}(A)$ by definition.

Polyominoes with holes might achieve the maximum perimeter. However, the next lemma checks the intuitive fact that the minimum perimeter cannot be achieved by polyominoes with holes.
\begin{lemma}\label{nop}
If $A \in \mathcal{A}_{n}$ and $A$ has at least one hole, then $p_{min}(n)<p(A)$.
\end{lemma}
\begin{proof}
By (\ref{perimeter1976}), an $n$-omino with $k$ holes, with $k>1$,  has an outer perimeter at least equal to $p_{min}(k+n)$. A polyomino with $k$ holes has a hole perimeter greater or equal to $4k$. Then, we have,
\begin{equation}\label{holesnominp}
p(A)=p_{o}(A)+p_{h}(A)\geq p_{min}(n+k)+4k.
\end{equation}
Because the function $h(x)=2\ceil{2\sqrt{x}}$ is a non-decreasing function, we can conclude from (\ref{holesnominp}) that,
\[p(A)\geq p_{min}(n+k)+4k > p_{min}(n)+4k.\]
\end{proof}

We denote this minimum outer perimeter over all polyominoes with $n$ tiles and $m$ holes (whenever such polyominoes exist) by $p_{\min}^{out}(n,m)$.
Note that we always have,
\begin{eqnarray}\label{pminpminout}
p_{\min}(n+m) \leq p_{\min}^{out}(n,m),
\end{eqnarray}
by definition.

\subsection{Main upper bound.}\label{SSmainub}

We prove in the next lemma a general result which allows us to generate an upper bound of $f$ from any lower bound. We will apply this lemma again in Sections \ref{FractalSequence} and \ref{Sgeneralb}.

\begin{lemma}\label{LTupperbound}
Let $n$ be any natural number. If $f(n)$ denotes the maximum number of holes that an element of $\mathcal{A}_n$ can have, then,
\begin{equation}\label{ub}
f(n) \leq \frac{4n-2b_{min}(n)-p_{min}(n+f(n))}{4}.
\end{equation}
\end{lemma}

\begin{proof}
Let $A$ be an element in $\mathcal{A}_{n}$ and let $h(A)$ denote the number of holes in $A$. Then,
\begin{equation}\label{upbound1}
p_{h}(A) = 4n-2b(A)-p_{o}(A)\leq  4n-2b_{min}(n)-p_{o}(A).
\end{equation}

By (\ref{pminpminout}) we have,
\begin{equation}\label{upbound2}
p_{o}(A)\geq p_{min}^{out}(n,h(A))\geq p_{min}(n+h(A)).
\end{equation}
From inequalities (\ref{upbound1}) and (\ref{upbound2}) we get,
\begin{equation}\label{upbound3}
p_{h}(A) \leq 4n-2b_{min}(n)-p_{min}(n+h(A)).
\end{equation}
Then, if $A$ is a polyomino that has $f(n)$ holes, we get,
\begin{equation*}
f(n)\leq \frac{p_{h}(A)}{4}\leq \frac{4n-2b_{min}(n)-p_{min}(n+f(n))}{4},
\end{equation*}
which establishes inequality (\ref{ub}).
\end{proof}

As a corollary for Lemma \ref{LTupperbound}, from any lower bound of $f$, we get the following upper bound for $f$,
\begin{corollary}\label{Cublb}
If $lb_f(n)\leq f(n)$, then
\begin{equation}\label{ublb}
f(n)\leq \frac{1}{2}n-\frac{1}{2}\ceil*{2\sqrt{n+lb_f(n)}}+\frac{1}{2},
\end{equation}
\end{corollary}
\begin{proof}
Let $lb_f(n)\leq f(n)$ for a natural number $n$, then $p_{min}(n+lb_f(n))\leq p_{min}(n+f(n))$, by monotonicity. This inequality allow us to obtain (\ref{ublb}) by substituting $p_{min}(n+f(n))$ with $p_{min}(n+lb_f(n))$ in  (\ref{ub}).
\end{proof}

\section{Polyominoes that attain the maximum number of holes.}\label{FractalSequence}
\subsection{Construction of the sequence.}
We are going to describe below how to construct a sequence $\{S_{k}\}_{k=1}^\infty$ of polyominoes with $n_k$ squares and $h(S_k) = h_k$ holes. Remember that $n_k$ and $h_k$ were defined in (\ref{Defn}) and (\ref{Defh}) respectively.

The first three elements, $S_1$, $S_2$, and $S_3$ of the sequence were shown in Figure \ref{Fonehole}.

To generate the rest of the sequence for $n \geq 2$, we follow the next general recursion process:

\begin{itemize}
\item First, place a rotation point in the center of the top right tile of $S_{n-1}$.
\item Then, rotate $S_{n-1}$ with respect to this point ninety degrees four times creating four, overlapping copies.
\item Finally, remove the tile containing the rotation point.
\end{itemize}

\input{Ifractalconstruction.txt}

\subsection{Properties of the sequence $S_k$ and proof of Theorem \ref{Ttremebundoexactly}.}
From this construction, we observe that for $k \ge 1$ we have the recursion $h(S_{k}) = 4h(S_{k-1})+1$. The factor of four is due to the four reflected copies of $S_{k-1}$ generated in the process of constructing $S_k$. The one hole added is generated by the square removed after the rotation process. Then, because $h(S_1)=1$, we get $h(S_k)= h_k$ for all $k\geq 1$.

Let $s_k$ be the number of tiles in $S_k$ for $k\geq 1$. The sequence $s_k$ satisfies the recursion, $$s_{k+1}=4s_k-4(2^k+1),$$ because the polyominoes $S_k$ have side lengths of $2^k +1$ tiles and, in the rotation process, $4(2^k+1)$ tiles overlap. Additionally, the sequence $n_k$ satisfies the relationship, $$n_{k+1}=4n_k-4(2^k+1).$$ Then, because both $n_k$ and $s_k$ satisfy the same recursion relationship and are equal in the first element ($s_1=8$, $n_1=8$), we can conclude that $n_k$ and $s_k$ are the same sequences.

We have proved the following lemma.
\begin{lemma}\label{Lsequencek}
There exists a sequence of polyominoes $\{S_k\}_{k=1}^\infty$, such that $S_k$ has $n_k$ tiles and $h_k$ holes.
\end{lemma}

\begin{proof}[Proof of Theorem \ref{Ttremebundoexactly}]
First we show that $f(n_k)=h_k$. From Lemma \ref{Lsequencek} we know,
\begin{equation}\label{eq}
h_k \leq f(n_k).
\end{equation} Substituting this lower bound in  (\ref{ublb}) we have,

\begin{equation}\label{pluginmaxexactly1}
f(n_k) \leq \frac{1}{2}n_k-\frac{1}{2}\ceil*{2\sqrt{n_k + h_k} } + \frac{1}{2}.
\end{equation}
From the easily verified identity,
\[h_k +\frac{1}{2} = \frac{1}{2}n_k-\frac{1}{2}\ceil*{2\sqrt{n_k+h_k}}+\frac{1}{2},
\] and inequalities (\ref{pluginmaxexactly1}) and (\ref{eq}), we get, \begin{equation}
h_k  \leq f(n_k) \leq h_k + \frac{1}{2}.
\end{equation}
This implies that,
\begin{equation}\label{Efnkhk}
f(n_k) = h_k,
\end{equation}
because $f(n_k)$ and $h_k$ are integers.

Now, we prove that $f(n_k-1) = h_k$. By removing the upper leftmost square from each $S_k$, it is possible to generate a sequence of polyominoes $\{A_k\}_{k=1}^\infty$ with $n_k-1$ tiles each $A_k$, such that $h(A_k) = h_k$. This implies that $h_k\leq f(n_k-1).$
Then, because $f(n_k-1) \leq f(n_k)$ and $f(n_k)=h_k$, we can conclude that,
 \[f(n_k-1)=h_k. \]

Finally, we prove that $f(n_k-2) = h_k-1$.
Because $f(n_k-1) = h_k$ and $f$ is a non decreasing function, it is forced that $f(n_k-2) \leq  h_k$. If we assume by way of contradiction that $f(n_k-2) = h_k$, then, using (\ref{ublb}) it must be the case that,
\begin{equation}
h_k=f(n_k-2)\leq \frac{1}{2}(n_k-2)-\frac{1}{2}\ceil*{2\sqrt{(n_k-2)+h_k}}+\frac{1}{2}.
\end{equation}
Substituting (\ref{Defn}) and (\ref{Defh}) in this inequality leads to a contradiction. Then $f(n_k-2)< h_k$; and, from this inequality and Lemma \ref{Ljumpsoff}, we can conclude that,
\begin{equation}
f(n_k-2) = h_k - 1.
\end{equation}

\end{proof}

\section{ General bounds}\label{Sgeneralb}
In this section we prove Theorem \ref{Ttremebundo}. We first prove the lower bound. We construct a sequence of polyominoes $\{R_k\}_{k=1}^\infty$ of polyominoes with  $m_k=40k^2+20k$ tiles and $t_k=20k^2$ holes, as follows.

We first place $10k^2$ copies of the pattern $S$ (in Figure \ref{pattern}) into a rectangle $6k$ high and $10k$ long. We add a top row of  $10k$ tiles, and a leftmost column of $6k-1$  tiles. Finally, we attach $2k$ vertically aligned dominoes (for a total of $4k$ tiles) to the rightmost column. The polyomino $R_2$ is depicted in Figure \ref{Fm2}, and the $40k^2$ tiles just described are colored with the lightest gray in this figure. The initial, repeated pattern $S$ is bordered in black within $R_2$.

Next ``we fill in the gaps'' between these constructions to a family of polyominoes $\{R_{k,l} \}$ with $m_{k,l} = 40k^2 + 20k + l$ tiles and $t_{k,l} = 20k^2 + \lfloor l / 2 \rfloor$ holes, defined whenever
$$0 \le l \le \sum_{i=1}^{2k-1} 2i = 2k(2k-1).$$
The polyomino $R_{k,l}$ is constructed from $R_{k}$ by adding tiles along the right side, continuing the domino pattern, as in Figure \ref{m1}. Every two tiles added creates one hole.

Note that $m_{k+1} - m_k  = 80k+ 60$, and $$80k + 60 \le 2k(2k-1)$$ for $k \ge 42$.
We define a sequence of polyominoes $\{ R'_n \}$ for all $n \geq m_{42}$, as follows.
Let $k$ be the largest number such that $m_k \le n$, let $l = n - m_k$, and then define $R'_n = R_{k,l}$.

Now we check that if $C_1 > \sqrt{5/2}$, then $f(n) \ge n/2 - C_1 \sqrt{n}$ for all large enough $n$.

The polyomino $R'_n= R_{k,l}$ has $m_{k,l}$ tiles and $t_{k,l}$ holes. Since $k$ and $l$ are nonnegative, we have
\begin{equation}
(5/2) \left(40k^2+ 20k+ l \right) = 100k^2 + 50k + (5/2) l \ge 100 k^2 + 20 k + 1.
\end{equation}
Taking square roots of both sides gives
\begin{equation}
\sqrt{5/2} \sqrt{40k^2 +20k+l } \ge 10 k+1.\\
\end{equation}

Since $C_1 > \sqrt{5/2}$, this gives
\begin{equation} \label{eqn:40k2}
C_1 \sqrt{40k^2 +20k+l } \ge 10 k+1.\\
\end{equation}

Then
\begin{equation} \label{eqn:20k2}
20k^2 + \left\lfloor \frac{l}{2} \right\rfloor \ge \frac{1}{2} \left( 40 k^2 + 20k + l \right)  - C_1  \sqrt{40k^2 +20k+l }.
\end{equation}

Inequality (\ref{eqn:20k2}) was obtained from (\ref{eqn:40k2}) by adding $20k^2$ to both sides, adding  the inequality $\left\lfloor \frac{l}{2} \right\rfloor  \ge \frac{l}{2} - 1$, and rearranging terms.

By considering the sequence $R'_n$, we have
\begin{equation} \label{eqn:mainlower}
f(n) \ge \frac{n}{2} - C_1 \sqrt{n}
\end{equation} for all large enough $n$, as desired.

Combining inequality (\ref{eqn:mainlower}) with Corollary \ref{Cublb}, we get an upper bound that allows us to complete the proof of Theorem  \ref{Ttremebundo}. Indeed, we have
\begin{align}
f(n) & \le \frac{1}{2} n -\frac{1}{2} \left\lceil 2 \sqrt{n + \frac{n}{2} - C_1 \sqrt{n}} \right\rceil + \frac{1}{2} \\
& \le  \frac{1}{2} n - \sqrt{\frac{3n}{2} - C_1 \sqrt{n} } + 1 \\
& \le \frac{1}{2} n - C_2 \sqrt{n}, \label{eqn:mainupper}
\end{align}
for $C_2 < \sqrt{3/2}$ and large enough $n$.

\input{patternS.txt}

\input{R2.txt}
\input{R3.txt}

\section{Concluding remarks}

Theorem \ref{Ttremebundoexactly} makes us guess that the upper bound in Theorem \ref{Ttremebundo} is correct, and that the lower bound might be improved to the following.

\begin{conjecture} For every $C_{1}>\sqrt{3 / 2}$ there exists an $n_0=n_0(C_1)$ such that
$$f(n) \ge \frac{1}{2} n - C_1 \sqrt{n}.$$
for all $n \ge n_0$. \end{conjecture}
It might even be possible to find an exact formula for $f$. The sequence of polyominoes $\{ S_k \}$ shows that
$$f(n) =\frac{1}{2} n - \sqrt{ \frac{3}{2} n +\frac{1}{4}} + \frac{1}{2},$$
infinitely often.

The recursive construction suggests that the main sequence of polyominoes $\{ S_k \}$ is approaching some limiting fractal shape. Elliot Paquette pointed out to us that one way to make sense of this idea is to consider the ``inner boundary'' of a polyomino $S_n$ to be an immersed circle in $\mathbb{R}^2$. Appropriately rescaling and reparameterizing, these circles seem to converge to a space-filling curve in $[0,1]^2$.

\input{circle.txt}

We noticed also that one can derive an aperiodic tiling of the plane from our main sequence, as follows. The `planar dual' of a polyomino $S_k$ is a planar graph with one vertex for every tile and one bounded face for every hole---see Figure \ref{Tess}. One can find a limiting infinite polyomino by centering every $S_k$ at the origin, and taking the union of all of them. The planar dual of this infinite polyomino is an aperiodic tiling by squares, pentagons, and hexagons.
\input{tessellation.txt}

\noindent

\section*{Acknowledgments.}
We thank CIMAT (Centro de Investigaci\'on en Matem\'aticas) for hosting us in Summer 2016 and ICERM (Institute for Computational and Experimental Research in Mathematics) for hosting us in Autumn 2016. This work was also supported in part by NSF DMS-1352386.

\bibliographystyle{plain}

redpajama_set_name: RedPajamaArXiv