# Multi-Material Topology Optimization with Continuous Magnetization Direction
for Permanent Magnet Synchronous Reluctance Motors

Thomas Gauthey E-Mail:thomas.gauthey(at)geeps.centralesupelec.fr  Université
Paris-Saclay, CentraleSupélec, CNRS, Laboratoire de Génie Electrique et
Electronique de Paris, 91192, Gif-sur-Yvette, France.
Sorbonne Université, CNRS, Laboratoire de Génie Electrique et Electronique de
Paris, 75252, Paris, France.  Technische Universität Graz, Institut für
Angewandte Mathematik, 8010 Graz, Austria.  Peter Gangl E-Mail:
gangl(at)math.tugraz.at Technische Universität Graz, Institut für Angewandte
Mathematik, 8010 Graz, Austria.  Maya Hage Hassan E-Mail: maya.hage-
hassan(at)centralesupelec.fr  Université Paris-Saclay, CentraleSupélec, CNRS,
Laboratoire de Génie Electrique et Electronique de Paris, 91192, Gif-sur-
Yvette, France.
Sorbonne Université, CNRS, Laboratoire de Génie Electrique et Electronique de
Paris, 75252, Paris, France.

###### Abstract

Permanent magnet-assisted synchronous reluctance motors (PMSynRM) have a
significantly higher average torque than synchronous reluctance motors. Thus,
determining an optimal design results in a multi-material topology
optimization problem, where one seeks to distribute ferromagnetic material,
air and permanent magnets within the rotor in an optimal manner.This study
proposed a novel density-based distribution scheme, which allows for
continuous magnetization direction instead of a finite set of angles. Thus, an
interpolation scheme is established between properties pertaining to magnets
and non-linear materials. This allows for new designs to emerge without
introducing complex geometric parameterization or relying on the user’s biases
and intuitions. Toward reducing computation time, Nitsche-type mortaring is
applied, allowing for free rotation of the rotor mesh relative to the stator
mesh. The average torque is approximated using only four-point static
positions. This study investigates several interpolation schemes and presents
a new one inspired by the topological derivative. We propose to filter the
final design for the magnetization angle using K-mean clustering accounting
for technical feasibility constraints of magnets. Finally, the design of the
electrical motor is proposed to maximize torque value.

Keywords: Topology optimization, Permanent magnets machines, Design
optimization, Acceleration methods

##  1 Introduction

Synchronous reluctance machines (SynRM) are standard in households and
industrial applications, thanks to their cheap cost compared to permanent
magnet motors and advances in manufacturing techniques. Although the
deployment of these machines continues [15], PMSynRM offers an excellent
alternative for both structures, solving for SynRM, its poor power factor and,
for permanent magnet machines (PMM), its cost. The design of these machines
using parametric optimization often necessitates either complex analytical
models or the use of Finite Element Analysis (FEA) relying heavily on
experienced engineers and known good designs [23].
Density based optimizations allow for bypassing such cumbersome frameworks.
Although they were first developed for two materials application in continuous
mechanics [33], a rise in n-materials optimization in the field of
electromagnetics has allowed for new PMM and SynRM to emerge [18, 14]. In most
optimizations where permanent magnets are involved, magnetization direction
are fixed [29] or limited to a set of a couple values [20, 6, 7]. If
continuous directions are considered during the optimization process, the
final design is filtered to take into account only a couple of predefined
directions to meet manufacturing constraints [36, 17].
In this paper we propose a simultaneous density-based optimization scheme
consisting of three-material (air-iron-magnet) with a continuous magnetization
direction. The proposition is applied to design the rotor of a distributed
winding stator as described in [25, 11] to maximize the mean torque under
constraints. The final designs are filtered using an unbiased K-means
heuristic for accounting for feasibility constraints. Here, we propose also to
accelerate the torque calculation through a four-point method.

##  2 Problem description

We chose to investigate a SynRM described in [11, 25], of which the rotor
design had proven to be a challenging problem for topology optimization and
use it for our PMSynRM optimization problem.

###  2.1 Geometry description

The electrical machine geometry and current density distributions are given
respectively in Figure 1 and 2. The dimensions for the considered machine are
given in Table 1.

Figure 1: Machine geometry  Table 1: Geometric parameters  Parameter | Value
---|---
Slot number | 24
Axial length | 50.0 mm
Outer rotor radius | 18.5 mm
Inner stator radius | 26.5 mm
Outer stator radius | 47.5 mm
Air gap length | 8.0 mm

This machine differs from most conventional SynRM by its large air gap which
constrains the statoric winding distribution to only one pair of poles (cf.
Figure 2).

Figure 2: Statoric winding distribution and current parameters Table 2: Statoric winding parameters  Parameter | Value
---|---
Number of turn $N_{s}$ | 64
Winding type | Distributed
Connection type | Star
Resistance ($R_{S,20^{\circ}C}$) | 7.1 $\Omega$
Voltage $U_{eff}$ | 230 V
Peak intensity $I_{max}$ | 12 A
Number of pole pairs $n_{pp}$ | 1

We introduce the relationship between the electrical angle $\theta_{elec}$ and
the mechanical angle $\theta$

$\theta_{elec}=n_{pp}\theta,\ $ (2.1)

with $n_{pp}$ the number of pair of poles, here $n_{pp}=1$.

We define the three phases as follows:

$\begin{cases}I_{U}(\theta)=I_{max}cos\left(n_{pp}\theta+\varphi\right)\\\
I_{V}(\theta)=I_{max}cos\left(n_{pp}\theta+\varphi-\frac{2\pi}{3}\right)\\\
I_{W}(\theta)=I_{max}cos\left(n_{pp}\theta+\varphi-\frac{4\pi}{3}\right).\end{cases}$
(2.2)

Here, $\varphi$ is the phase angle. The computational domain $\Omega$ consists
of iron, air, permanent magnet and coils,

$\Omega=\Omega_{f}\cup\Omega_{air}\cup\Omega_{mag}\cup\Omega_{c}$ (2.3)

where we further subdivide the ferromagnetic and air subdomains into their
rotor and stator parts,

$\Omega_{f}=\Omega_{f,stat}\cup\Omega_{f,rot},\qquad\Omega_{air}=\Omega_{air,stat}\cup\Omega_{air,rot}.$
(2.4)

Moreover, we subdivide the coil subdomains according to the distribution shown
in Figure 2,

$\Omega_{c}=\Omega_{U^{+}}\cup\Omega_{U^{-}}\cup\Omega_{V^{+}}\cup\Omega_{V^{-}}\cup\Omega_{W^{+}}\cup\Omega_{W^{-}}.$
(2.5)

###  2.2 Partial differential equation

In the two-dimensional magnetostatic setting, the magnetic flux density
$\mathbf{B}=\mbox{curl}((0,0,u)^{\top})$ for rotor position
$\theta\in[0,2\pi]$ can be computed via the solution of the boundary value
problem

$\displaystyle\mbox{Find }u\in
H^{1}_{0}(\Omega):\int_{\Omega}\nu_{\theta}(x,|\nabla u|)\nabla u\cdot\nabla
v\,dx=$ (2.6)
$\displaystyle\int_{\Omega_{c}}j(\theta)\,v\,dx+\int_{\Omega_{mag}^{\theta}}R_{\theta}\begin{bmatrix}-M_{y}\\\
M_{x}\end{bmatrix}\cdot\nabla v\,dx$

for all $v\in H^{1}_{0}(\Omega)$, see e.g. [26].
Here, the magnetic reluctivity is a nonlinear function $\hat{\nu}$ of the flux
density $|\mathbf{B}|=|\nabla u|$ in the ferromagnetic subdomain and a
constant $\nu_{0}=10^{7}/(4\pi)$ elsewhere, i.e.,

$\nu_{\theta}(x,|\nabla u|)=\begin{cases}\hat{\nu}(|\nabla
u|)&x\in\Omega_{f}^{\theta}\\\
\nu_{0}&x\in\Omega_{air}^{\theta}\cup\Omega_{c}\cup\Omega_{mag}^{\theta}\end{cases}$
(2.7)

with the rotated domains

$\displaystyle\Omega_{f}^{\theta}=$ $\displaystyle\Omega_{f,stat}\cup
R_{\theta}\Omega_{f,rot}$ (2.8) $\displaystyle\Omega_{air}^{\theta}=$
$\displaystyle\Omega_{air,stat}\cup R_{\theta}\Omega_{air,rot}$ (2.9)
$\displaystyle\Omega_{mag}^{\theta}=$ $\displaystyle R_{\theta}\Omega_{mag}$
(2.10)

and $R_{\theta}$ a rotation matrix around angle $\theta$,

$R_{\theta}=\begin{bmatrix}\mbox{cos}\theta&-\mbox{sin}\theta\\\
\mbox{sin}\theta&\mbox{cos}\theta\end{bmatrix}.$ (2.11)

The first term on the right hand side of (2.6) represents the impressed
current density which is given by

$\begin{split}j(x,\theta)=&\;\chi_{\Omega_{U^{+}}}(x)j_{U}(\theta)+\chi_{\Omega_{V^{+}}}(x)j_{V}(\theta)+\chi_{\Omega_{W^{+}}}(x)j_{W}(\theta)\\\
&-\chi_{\Omega_{U^{-}}}(x)j_{U}(\theta)-\chi_{\Omega_{V^{-}}}(x)j_{V}(\theta)-\chi_{\Omega_{W^{-}}}(x)j_{W}(\theta),\end{split}$
(2.12)

where $\chi_{A}$ denotes the characteristic function of a set $A$,

$\chi_{A}(x)=\begin{cases}1&x\in A,\\\ 0&\mbox{else}.\end{cases}$

Here the current distribution is defined by :

$j_{p}(\theta)=\frac{1}{S_{slot}}N_{s}I_{p}(\theta),\quad p\in\\{U,V,W\\}$
(2.13)

with $S_{slot}$ the cross-sectional area of one coil, $N_{s}$ the number of
turns per coil and $I_{U}$, $I_{V}$, $I_{W}$ as defined in (2.2). The second
term on the right hand side of (2.6) represents the magnetization
$\mathbf{M}=(M_{x},M_{y})^{\top}$ coming from permanent magnets which will be
added in the course of the multi-material optimization procedure.

In the following, we will denote by $u_{\theta}$ the solution to the state
equation (2.6).
We present here after the properties of interest of the materials
(air,ferromagnetic,magnet) used in the machine.

Table 3: Material properties  Material | Reluctivity [$m.H^{-1}$] | Magnetization [$A.m^{-1}$]
---|---|---
Air | $\nu_{0}$ | 0
Copper | $\nu_{0}$ | 0
Ferromagnetic | $\hat{\nu}(|\vec{B}|)$ | 0
Magnet | $\nu_{0}$ | $M_{max}$

The maximum norm of the magnetization vector was chosen as $M_{max}=2.33\cdot
10^{5}A.m^{-2}$ to fit data from [28] on ferrite magnets. The reluctivity of
the magnets and of the copper coils is assimilated to the one of air to
simplify further material interpolation and avoid complex schemes like the
ones found in [37]. The non-linear behaviour of the ferromagnetic material is
modelled with a Marrocco’s BH curve approximation [24].

$\hat{\nu}(|\vec{B}|)=\begin{cases}\nu_{0}(\varepsilon+\frac{(c-\varepsilon)|\vec{B}|^{2\alpha}}{\tau+|\vec{B}|^{2\alpha}}`&\mbox{if
}|\vec{B}|\leq B_{max},\\\
\nu_{0}\left(1-\frac{M_{s}}{|\vec{B}|}\right)&\mbox{else if
}|\vec{B}|>B_{s},\\\
exp\left(\frac{\gamma\left(|\vec{B}|-\beta\right)}{|\vec{B}|}\right)&\mbox{otherwise},\end{cases}$
(2.14)

where $B_{s}=\beta+\frac{log\left(\frac{\nu_{0}}{\gamma}\right)}{\gamma}$ and
$M_{s}=B_{s}+\frac{1}{\gamma}$ and the coefficient of the Marrocco curve in
Figure3 are defined in the Table 4.

Table 4: Marrocco curve coefficient for the ferromagnetic material Parameter |  | Value
---|---|---
$\alpha$ |  | 6.84
$\beta$ |  | -1.30$\cdot 10^{-1}$
$\gamma$ |  | 4.86
$\varepsilon$ |  | 1.57 $\cdot 10^{-4}$
$\tau$ |  | 4.14 $\cdot 10^{3}$
c |  | 1.90$\cdot 10^{-2}$
$B_{max}$ |  | 1.80 (T)

(a) Magnetic permeability

(b) Marrocco BH-curve

Figure 3: Non-linear characteristic for the ferromagnetic material.

###  2.3 Torque computation method

For computing the torque, we chose a method based on Maxwell’s stress tensor,
Arkkio’s method [30]. While co-energy torque computation methods were proven
to be more precise and less prone to error, they are more costly in terms of
computation time and not fit for optimization. Using Arkkio’s method, the
torque can be computed as

$T=\frac{L_{z}\nu_{0}}{r_{s}-r_{r}}\int_{S}\sqrt{x^{2}+y^{2}}B_{r}B_{\phi}dS$
(2.15)

where $B_{r}$ and $B_{\phi}$ denote the radial and tangential magnetic
induction, respectively, $L_{z}$ denotes the length of the machine in
$z$-direction and $S$ denotes the surface between radii $r_{s}$ and $r_{r}$ in
the air gap (with $r_{s}>r_{r}$).

In the setting of two-dimensional magnetostatics, the torque for the rotor
position given by angle $\theta$ thus amounts to

$T_{\theta}=T(u_{\theta})=\frac{L_{z}\nu_{0}}{r_{s}-r_{r}}\int_{S}Q\nabla
u_{\theta}\cdot\nabla u_{\theta}\,dS$ (2.16) $\mbox{with
}Q(x,y)=\begin{bmatrix}\frac{xy}{\sqrt{x^{2}+y^{2}}}&&\frac{y^{2}-x^{2}}{2\sqrt{x^{2}+y^{2}}}\\\
\frac{y^{2}-x^{2}}{2\sqrt{x^{2}+y^{2}}}&&\frac{-xy}{\sqrt{x^{2}+y^{2}}}\end{bmatrix}\in\mathbb{R}^{2\times
2}.$ (2.17)

Determining the average torque by means of its instantaneous values can be
very expensive. It is shown in [5] that a good approximation to the average
torque can be obtained when evaluating the torque for only suitably chosen
rotor positions,

$\bar{T}=\frac{1}{4}\left(T_{0}+T_{\frac{\pi}{12}}+T_{\frac{\pi}{6}}+T_{\frac{\pi}{4}}\right).$
(2.18)

We compared the average torque obtained by evaluation at 500 equally
distributed rotor positions between $0$ and $2\pi$ with the value obtained by
the four-point formula (2.18).

When the torque value is not equal to zero, the error found was to be lower
than 0.1% as expected and described in literature [1]. This is solved beyond
the first iteration.

Table 5: Four static positions method error

Design | $\bar{T}[N.m]$ | $\bar{T}[N.m]$ | Error
---|---|---|---
| (500 points) | (4points) | [%]
Unbiased starting point* | 1.5790 $\cdot 10^{-6}$ | 5.7232 $\cdot 10^{-6}$ | 262.4
Final design Table 7 | 1.1123 | 1.1129 | 0.048
Final design Table 9 | 1.4513 | 1.4516 | 0.027

* in this design $\rho_{\nu}=0.5,\rho_{M_{x}}=0.5,\rho_{M_{y}}=0.5$ everywhere in the rotor

##  3 Optimization problem

In this section, we define our optimization problem and reformulate the
forward problem to fit the density-based topology optimization approach. Our
goal is to maximize the average torque computed via (2.18),

$(P_{1}):\begin{cases}\text{maximize
}\bar{T}=\frac{1}{4}\left(T(u_{0})+T(u_{\frac{\pi}{12})}+T(u_{\frac{\pi}{6}})+T(u_{\frac{\pi}{4}})\right)\\\
\mbox{s.t. $u_{\theta}$ is a solution of \eqref{eq:2dMag} for
}\theta\in\\{0,\frac{\pi}{12},\frac{\pi}{6},\frac{\pi}{4}\\}\end{cases}$ (3.1)

This is achieved by finding the optimal material distribution consisting of
ferromagnetic material, air and permanent magnets on the one hand, and the
optimal magnetization direction of the permanent magnets on the other hand.
Moreover, we will incorporate a bound on the maximum allowed permanent magnet
volume.

###  3.1 Density based topology optimization

Let us reformulate the forward problem (2.6), introducing the three density
variables respectively for the ferromagnetic material and the two components
of the permanent magnets magnetization, $\rho_{\nu},\rho_{M_{x}},\rho_{M_{y}}$
defined in
$\Omega^{\theta}_{rot}=R_{\theta}(\Omega_{f,rot}\cup\Omega_{air,rot}\cup\Omega_{mag}$).
Moreover we introduce the rotated design variables

$\displaystyle\rho_{\nu}^{\theta}(x,y)=$
$\displaystyle\rho_{\nu}(R_{\theta}((x,y)^{\top}))$
$\displaystyle\rho_{M_{x}}^{\theta}(x,y)=$
$\displaystyle\rho_{M_{x}}(R_{\theta}((x,y)^{\top}))$
$\displaystyle\rho_{M_{y}}^{\theta}(x,y)=$
$\displaystyle\rho_{M_{y}}(R_{\theta}((x,y)^{\top}))$

which represent the design given by $\rho_{\nu}$, $\rho_{M_{x}}$,
$\rho_{M_{y}}$ after rotation, and the vector of design variables
$\textbf{X}:=\begin{bmatrix}\rho_{\nu},&\rho_{M_{x}},&\rho_{M_{y}}\end{bmatrix}^{\top}$.

Given two interpolation functions

$f_{\nu}:[0,1]\rightarrow[0,1],\qquad f_{M}:[0,1]\rightarrow[0,1],$ (3.2)

we define the operator

$\displaystyle K_{\theta}:(\mathbf{X},u,v)\mapsto$
$\displaystyle\int_{\Omega}\nu(\rho_{\nu}^{\theta},|\nabla u|)\nabla
u\cdot\nabla v$ (3.3)
$\displaystyle-\int_{\Omega_{rot}}f_{\nu}(1-\rho_{\nu}^{\theta})\frac{M_{max}f_{M}(|\vec{M}^{\theta}|)}{|\vec{M}^{\theta}|}R_{\theta}\begin{bmatrix}-M_{y}^{\theta}\\\
M_{x}^{\theta}\end{bmatrix}\cdot\nabla v,$

with the reluctivity function

$\nu(\rho_{\nu}^{\theta},|\nabla u|)=\begin{cases}\hat{\nu}(|\nabla u|)&\text{
in }\Omega_{f,stat}\\\ \nu_{0}&\text{ in }\Omega_{c}\cup\Omega_{air,stat}\\\
\nu_{0}+f_{\nu}(\rho_{\nu}^{\theta})(\hat{\nu}(|\nabla u|)-\nu_{0})&\text{ in
}\Omega_{rot}\end{cases}$ (3.4)

and with the components of the magnetization vector
$\vec{M}^{\theta}=(M_{x}^{\theta},M_{y}^{\theta})$ given in dependence of the
two rotated density variables $\rho_{M_{x}}^{\theta}$,
$\rho_{M_{y}}^{\theta}$,

$(M_{x}^{\theta},M_{y}^{\theta})=\tilde{f}_{sd}(\rho_{M_{x}}^{\theta},\rho_{M_{y}}^{\theta})$
(3.5)

for a mapping $\tilde{f}_{sd}$ which will be discussed later on. Hence, the
state equation (2.6) can be reformulated into

$\displaystyle\text{Find }u\in H^{1}_{0}(\Omega):\,$ (3.6) $\displaystyle
K_{\theta}(u,v,\textbf{X})=\int_{\Omega_{c}}j(\theta)\,v\,dx,\text{ for all
}v\in H^{1}_{0}(\Omega).$

The optimization problem ($P_{1}$) can then be reformulated into

$(P_{2}):\begin{cases}\underset{\textbf{X}}{\text{maximize
}}\bar{T}=\frac{1}{4}\sum_{\theta\in\\{0,\frac{\pi}{12},\frac{\pi}{6},\frac{\pi}{4}\\}}T(u_{\theta})\\\
\text{s.t. $u_{\theta}$ is a solution of the
\eqref{eq:PDE_density}}\quad\theta\in\\{0,\frac{\pi}{12},\frac{\pi}{6},\frac{\pi}{4}\\}\end{cases}$
(3.7)

####  3.1.1 Material interpolation schemes

In density based topology optimization, the quality of the final solution is
dependant on the choice of interpolation functions (equation (3.2)). We
present here two existing schemes and a novel one based on properties of the
topological derivative.
The polynomial interpolation scheme

$f_{n}(\rho)=\rho^{n}\quad n>0,$ (3.8)

also referred to as SIMP (Solid Isotropic Material with Penalization), is the
most used material interpolation scheme for topology optimization and allows
for easy penalization of intermediate materials. However, it presents some
symmetry issues and favors low $\rho$ associated material in the final design.
In [31], the authors compared this scheme to other schemes and concluded that
the final design was not as good as many other proposed ones.

Figure 4: SIMP Polynomial interpolation scheme

To solve symmetry issues introduced by the classical polynomial interpolation,
D. Lukàš introduced a new scheme in [21]:

$f_{\lambda}(\rho)=\frac{1}{2}\left(1+\frac{1}{\arctan(\lambda)}\arctan(\lambda(2\rho-1))\right),\quad\lambda>0.$
(3.9)

In this equation the particular invariant point $\rho=0.5$ does not promote
intermediate materials, grey material depends on $\lambda$ values (cf. Figure
5).
High $\lambda$ values permit to penalize intermediate materials but can lead
to a poor convergence of the algorithm. A parameter study for $\lambda$ led us
to choose $\lambda=5$. This interpolation method is chosen for the norm of the
magnetization vector ($f_{M}$ in (LABEL:eq:PDE_density)).

Figure 5: D. Lukàš’s interpolation scheme

Finally, we propose a new interpolation scheme as given in Figure 6,which is
inspired by the topological derivative as done in [2], see also the the SIMP-
All method for linear elasticity [9]. Here, we seek to design a material
interpolation function whose derivative with respect to the density variable
$\rho$ coincides with the topological derivative of the problem at $\rho=0$
and $\rho=1$. When interpolating between two linear materials with reluctivity
values $\nu_{0}$ and $\nu_{1}$, the conditions for the material interpolation
function $f$ according to [2] would read

$\begin{cases}f(0)=0,\\\ f(1)=1,\\\
f^{\prime}(0)=\frac{2\nu_{0}}{\nu_{0}+\nu_{1}},\\\
f^{\prime}(1)=\frac{2\nu_{1}}{\nu_{0}+\nu_{1}}.\end{cases}$ (3.10)

Due to the involved formula of the topological derivative for nonlinear
magnetostatics [3], a mathematically rigorous extension of this method to the
nonlinear setting is not straightforward. However, inspired by the particular
behaviour of the Marrocco BH-curve where the magnetic reluctivity is almost
constant for low flux density values, see Figure 3, we simply use this idea
for that constant reluctivity value $\nu_{1}:=\nu_{0}\varepsilon\approx
124.94$. Using cubic Hermite interpolation for the conditions (3.10), we
obtain the polynomial

$f(\rho)=\frac{2\nu_{0}}{\nu_{0}+\nu_{1}}\rho-\frac{\nu_{0}-\nu_{1}}{\nu_{0}+\nu_{1}}\rho^{2}.$
(3.11)

Note that the term of order 3 happens to vanish.

Figure 6: Topological derivative inspired interpolation scheme

(a) Interpolated BH curve for $\rho=0.25$

(b) Interpolated BH curve for $\rho=0.75$

(c) Interpolated BH curve for $\rho=0.9$

Figure 7: Interpolated BH curves using material interpolation scheme (3.11).

####  3.1.2 Magnetization vector transform

We deal with two magnetization density variables $\rho_{M_{x}}$ and
$\rho_{M_{y}}$ in order to represent the magnetization direction
$(M_{x},M_{y})$. One way of relating these quantities to each other would be
to have $\rho_{M_{x}}$ represent the first and $\rho_{M_{y}}$ the second
coordinate, resulting in a representation in Cartesian coordinates, which was
also considered in [36]. In this case, however, some magnetization directions
exhibit higher maximum magnetization than others, e.g.
$\rho_{M_{x}}=\rho_{M_{y}}=1$ would correspond to
$|(M_{x},M_{y})^{\top}|=\sqrt{2}$ whereas for the magnetization direction
pointing to the right $\rho_{M_{x}}=1$, $\rho_{M_{y}}=0.5$ would yield a
maximum magnetization of $|(M_{x},M_{y})^{\top}|=1$, thus making the maximum
magnetization angle dependent. As an alternative, one could use polar
coordinates and represent the magnetization direction by just one periodic
density function. In this case, however, the ambiguity of angular values
causes problems in the gradient computation. We define a means of solving
these issues without resorting to polar coordinates. We decide on two density
variables $\rho_{M_{x}},\rho_{M_{y}}$ with values in $[0,1]$, but map them
onto a disk, thereby avoiding an angle dependent maximum magnetization value.

Figure 8: Square to disk transform for the Magnetization vector coordinates

There are several mappings that approximately realize such a square-to-disk
transformation. In order to preserve the angle,as much as possible, without
being too computationally heavy, the elliptic grid mapping.

$f_{sd}(x,y)=\begin{cases}x\sqrt{1-\frac{y^{2}}{2}}\\\
y\sqrt{1-\frac{x^{2}}{2}}\end{cases}\text{with }(x,y)\in[-1,1]^{2}.$ (3.12)

was chosen as a good compromise out of the methods described in [10]. While
the associated inverse transformation $f_{ds}=f_{sd}^{-1}$ given by

$f_{ds}(u,v)=\begin{cases}\frac{1}{2}\left(\sqrt{2+u^{2}-v^{2}+2\sqrt{2}u}-\sqrt{2+u^{2}-v^{2}-2\sqrt{2}u}\right)\\\
\frac{1}{2}\left(\sqrt{2-u^{2}+v^{2}+2\sqrt{2}v}-\sqrt{2-u^{2}+v^{2}-2\sqrt{2}v}\right)\end{cases}$
(3.13)

is computationally more costly, it is only used once per iteration and in
post-processing. The mapping between the magnetization density variables
$\rho_{M_{x}},\rho_{M_{y}}$ and the magnetization vector
$\vec{M}=(M_{x},M_{y})$ (3.5) is then given by

$\displaystyle\tilde{f}_{sd}(\rho_{M_{x}}^{\theta},\rho_{M_{y}}^{\theta})=f_{sd}(2(\rho_{M_{x}}^{\theta}-0.5),2(\rho_{M_{y}}^{\theta}-0.5))$
(3.14)

###  3.2 Incorporation of volume constraints

Constraints on iron and magnets volume are added to avoid having structures
with disproportionate volumes of material. Hence we add new constraints with
the operator

$I_{vol}:\rho\mapsto\frac{1}{V_{\Omega_{rot}}}\int_{\Omega_{rot}}\rho(x)\,dx$
(3.15)

representing the volume fraction inside the rotor domain $\Omega_{rot}$ of a
material given by a density function $\rho$. Here,
$V_{rot}=\int_{\Omega_{rot}}1\;dx$ denotes the total area of $\Omega_{rot}$.
Based on ($P_{2}$) we define a new constrained optimization problem

$(P_{3}):\begin{cases}\text{maximize
}\bar{T}=\frac{1}{4}\sum_{\theta\in\\{0,\frac{\pi}{12},\frac{\pi}{6},\frac{\pi}{4}\\}}T(u_{\theta})\\\
\text{s.t. $u_{\theta}$ is a solution of
$\eqref{eq:PDE_density},\;\theta\in\\{0,\frac{\pi}{12},\frac{\pi}{6},\frac{\pi}{4}\\}$}\\\
\quad\;\text{ $I_{vol}(\rho_{\nu})\leq f_{v,f}$}\\\ \quad\;\text{
$I_{vol}((1-\rho_{\nu})\,|\vec{M}|)\leq f_{v,mag}$}\end{cases}$ (3.16)

with given upper bounds on the allowed ferromagnetic and permanent magnet
material $f_{v,f}$, $f_{v,mag}\in[0,1]$, respectively. We reformulate the
inequality constraints of ($P_{3}$) using the augmented Lagrangian framework
as described in [27],

$(P_{4}):\begin{cases}\text{minimize
}L(\mathbf{X},\underline{u})=-\bar{T}(\underline{u})+\psi(h_{v,f}(\mathbf{X}),\gamma_{f},\mu)\\\
\quad\quad\quad\quad\quad\quad\quad\quad+\psi(h_{v,mag}(\mathbf{X}),\gamma_{mag},\mu)&\\\
\text{s.t. $u_{\theta}$ is a solution of
\eqref{eq:PDE_density}},\theta\in\\{0,\frac{\pi}{12},\frac{\pi}{6},\frac{\pi}{4}\\}\\\
\end{cases}$ (3.17)

With the state vector

$\displaystyle\underline{u}:=(u_{0},u_{\frac{\pi}{12}},u_{\frac{\pi}{6}},u_{\frac{\pi}{4}}).$
(3.18)

Here,

$\begin{cases}h_{v,f}(\mathbf{X})=f_{v,f}-I_{vol}(\rho_{\nu}),\\\
h_{v,mag}(\mathbf{X})=f_{v,mag}-I_{vol}((1-\rho_{\nu})\,|\vec{M}|),\\\
\end{cases}$ (3.19)

with the scalar function

$\psi(t,\sigma,\mu)=\begin{cases}-\sigma t+\frac{1}{2\mu}t^{2}&\text{if
}t-\mu\sigma\leq 0,\\\ -\frac{\mu}{2}\sigma^{2}&\mbox{otherwise}.\end{cases}$
(3.20)

The positive scalar multipliers $\gamma_{f},\gamma_{mag},\mu$ are updated
according to the LANCELOT-Method of Multipliers presented in [27].

###  3.3 Adjoint Method

To solve the optimization problem $(P_{4})$ as formulated in (3.17) by a
gradient descent algorithm, we introduce the Lagrangian for the PDE-
constrained problem (3.17)

$\mathcal{L}(\mathbf{X},\underline{u},\underline{w})=L(\mathbf{X},\underline{u})+\sum_{\theta\in\\{0,\frac{\pi}{12},\frac{\pi}{6},\frac{\pi}{4}\\}}K_{\theta}(\mathbf{X},u_{\theta},w_{\theta})-\int_{\Omega_{stat}}j(\theta)w_{\theta},$
(3.21)

where
$\underline{w}=(w_{0},w_{\frac{\pi}{12}},w_{\frac{\pi}{6}},w_{\frac{\pi}{4}})$
is a vector of Lagrange multipliers. The adjoint states $\lambda_{\theta}$
corresponding to problem (3.17) for different rotor positions $\theta$ are
defined by $\frac{\partial\mathcal{L}}{\partial
u_{\theta}}(\mathbf{X},u_{\theta},\lambda_{\theta})=0$, i.e.
$\lambda_{\theta}$ is the solution to

$\mbox{Find }\lambda_{\theta}\in H_{0}^{1}(\Omega):\frac{\partial
K_{\theta}}{\partial
u}(\mathbf{X},u_{\theta},\lambda_{\theta})(v)=-\frac{1}{4}\frac{\partial
T}{\partial u}(u_{\theta})(v)\quad\\\ \mbox{for all }v\in H_{0}^{1}(\Omega).$
(3.22)

Here, the left and right hand side, respectively, involve the operators

$\displaystyle\frac{d}{du}K_{\theta}(\mathbf{X},u,\lambda)(v)=$
$\displaystyle\frac{d}{du}\left(\int_{\Omega}\nu(\rho_{\nu}^{\theta},|\nabla
u|)\nabla u\cdot\nabla\lambda\right)(v)$ (3.23) $\displaystyle=$
$\displaystyle\int_{\Omega}\nu(\rho_{\nu}^{\theta},|\nabla u|)\nabla
v\cdot\nabla\lambda+\int_{\Omega}\frac{d}{du}\nu(\rho_{\nu}^{\theta},|\nabla
u|)(v)\nabla u\cdot\nabla\lambda$ $\displaystyle=$
$\displaystyle\int_{\Omega}\nu(\rho_{\nu}^{\theta},|\nabla u|)\nabla
v\cdot\nabla\lambda+\int_{\Omega}f_{\nu}(\rho_{\nu}^{\theta})\frac{\hat{\nu}^{\prime}(|\nabla
u|)}{|\nabla u|}(\nabla u\cdot\nabla v)(\nabla u\cdot\nabla\lambda)$
$\displaystyle\frac{d}{du}T(u)(v)=$
$\displaystyle\;2\frac{L_{z}\nu_{0}}{r_{s}-r_{r}}\int_{S}Q\nabla u\cdot\nabla
v\;dS.$

Similarly to (3.18), we introduce the adjoint vector

$\underline{\lambda}:=(\lambda_{0},\lambda_{\frac{\pi}{12}},\lambda_{\frac{\pi}{6}},\lambda_{\frac{\pi}{4}}).$
(3.24)

Hence, the sensitivity associated with (3.21) amounts to

$\begin{split}\frac{\partial\mathcal{L}}{\partial\mathbf{X}}(\mathbf{X},\underline{u},\underline{\lambda})=&\sum_{\theta\in\\{0,\frac{\pi}{12},\frac{\pi}{6},\frac{\pi}{4}\\}}\frac{\partial
K_{\theta}}{\partial\mathbf{X}}(\mathbf{X},u_{\theta},\lambda_{\theta})\\\
&+\frac{\partial\psi(h_{v,f}(\mathbf{X}),\gamma_{f},\mu)}{\partial\mathbf{X}}+\frac{\partial\psi(h_{v,mag}(\mathbf{X}),\gamma_{mag},\mu)}{\partial\mathbf{X}}\end{split}$
(3.25)

where we used that $\frac{\partial T}{\partial\mathbf{X}}=0$ and
$\frac{\partial j}{\partial\mathbf{X}}=0$.

###  3.4 Update method

Several methods exist to consider the movement in electrical motors, such as
the Moving Band (MB) technique described in [8]. Even when using high order
elements in the MB, this method remains less accurate than the mortar element
method [4]. A variant of the mortar method, the Nitsche method, is chosen to
take into account the rotation [16]. The operator $K_{\theta}$ will be
replaced by $K_{\theta}^{NM}$ and $\mathcal{L}$ by $\tilde{\mathcal{L}}$ to
fit the new formulation (detailed in Appendix.A).

We introduce a triangular mesh with a total of $N$ elements, with $N_{rot}$
many elements inside the rotor domain $\Omega_{rot}$. We use piecewise linear
and globally continuous finite element basis functions to solve the state and
adjoint equations and we represent the density variables
$\rho_{\nu},\rho_{M_{x}},\rho_{M_{y}}$ as piecewise constant functions on the
mesh corresponding to $\Omega_{rot}$. Thus, these density variables can be
represented by a vector of dimension $N_{rot}$ consisting of the values of the
discrete functions in each element. We will use the same notation
$\rho_{\nu},\rho_{M_{x}},\rho_{M_{y}}$ for the vectors representing the
discretized density variables.

To comply with the bounds of the density variables, we introduce the projected
gradient as defined in [27, p. 520]

$\mathcal{G}:=\begin{bmatrix}\mathbf{P}_{\rho_{\nu},[0,1]}\left(\frac{\partial\tilde{\mathcal{L}}}{\partial\rho_{\nu}}\right)&\mathbf{P}_{\rho_{M_{x}},[0,1]}\left(\frac{\partial\tilde{\mathcal{L}}}{\partial\rho_{M_{x}}}\right)&\mathbf{P}_{\rho_{M_{y}},[0,1]}\left(\frac{\partial\tilde{\mathcal{L}}}{\partial\rho_{M_{y}}}\right)\end{bmatrix}^{\top},$
(3.26)

with the projection operator
$\mathbf{P}_{\rho,[a,b]}:\mathbb{R}^{N_{rot}}\rightarrow\mathbb{R}^{N_{rot}}$
defined by

$\left(\mathbf{P}_{\rho,[a,b]}(X)\right)_{i}=\begin{cases}X_{i}&\mbox{if
}\rho_{i}\in]a,b[,\\\ min(0,X_{i})&\mbox{if }\rho_{i}=a,\\\
max(0,X_{i})&\mbox{if }\rho_{i}=b,\\\ \end{cases}$ (3.27)

for $i=1,...,N_{rot}$. This in itself allows for finer geometry to emerge by
amplifying the relative importance of the gradient where change in the density
function can be made.
From the sensitivity we derive the update equation at iteration n

$\mathbf{X}_{i}^{n+1}=Q_{[0,1]}\left(\mathbf{X}_{i}^{n}-s\frac{\mathcal{G}_{i}}{|\mathcal{G}_{i}|}\right),\quad\mbox{for
}i=1,...,N_{rot}$ (3.28)

with the projection operator
$Q_{[a,b]}:\mathbb{R}^{3N_{rot}}\rightarrow\mathbb{R}^{3N_{rot}}$ defined by

$\left(Q_{[a,b]}(\mathbf{v})\right)_{i}=max(a,min(b,v_{i}))$ (3.29)

for a vector
$\mathbf{v}=(v_{1},\dots,v_{N_{rot}})^{\top}\in\mathbb{R}^{3N_{rot}}$ and
$i=1,\dots,N_{rot}$, to enforce the bounds of the density variables. Here, $s$
denotes the step size which is chosen in such a way that a descent of the
augmented Lagrangian is obtained,

$L(\mathbf{X}^{n+1},\underline{u}^{n+1})<L(\mathbf{X}^{n},\underline{u}^{n})$
(3.30)

with $\underline{u}^{n}$ is the vector of states for the design represented by
$\mathbf{X}^{n}$.

###  3.5 Filtering and projection method

In density-based topology optimization, checkerboard patterns and small
isolated elements of one material are avoided using filtering methods at each
step of gradient descent. While this filtering has a regularizing effect on
the density variables, it may introduce more gray areas. Therefore, in the
next step, so-called projection methods are applied in order to get a more
defined shape. This combination allows achieving smoother and more defined
boundaries between the materials in the final design.

In our approach, in the first step we perform density filtering by solving the
PDE

$-r^{2}\nabla^{2}\rho+\rho=\rho_{ref}\quad r>0.$ (3.31)

with $\rho_{ref}$ the given density function, which is commonly referred to as
Helmholtz filtering [19]. This mesh-independency filter modifies the
sensitivity by averaging on the neighbor cells. The parameter $r$ is the
radius of influence of the filtering and is defined in our case as a factor
$\delta$ of the minimum representative mesh-element length $h$, i.e.
$r=\delta\;h$.

For the projection step, we choose the function proposed in [35],

$\displaystyle
f_{\rho_{cut},\beta}(\rho)=\frac{\tanh(\beta(\rho-\rho_{cut}))+\tanh(\beta\rho_{cut})}{\tanh(\beta(1-\rho_{cut}))+\tanh(\beta\rho_{cut})}$
(3.32) $\displaystyle\rho_{cut}\in[0,1];\beta>0$

with $\rho_{cut}=0.5$ such as not to favor one material. The other parameters
$\delta$ and $\beta$ are chosen to preserve the equilibrium between the two
parts of the filtering step.

Figure 9: Comparison of projection curves for a varying $\beta$ parameter.

We applied the filtering technique (3.31) and the projection technique (3.32)
for the material densities. The same PDE-based filter (3.31) is applied for
the densities $\rho_{M_{x}}$ and $\rho_{M_{y}}$, to favor uniform
magnetization direction in a magnet area while the projection (3.32) is only
applied on $|\vec{M}|$ to avoid scaling issues.

###  3.6 Direct penalization of intermediate materials

Some final designs can still present fuzzy boundaries and intermediate
material, especially if the optimization starting point is near a local
minimum. To help overcome this issue, we propose to penalize the intermediate
materials directly as done in the phase-field topology optimization method
[13] and add to the cost function the following term with a weight $\gamma>0$:

$I_{\gamma}(\rho)=\frac{4\gamma}{V_{\Omega_{rot}}}\int_{\Omega_{rot}}\rho(x)(1-\rho(x))dx.$
(3.33)

The penalization is only applied on iron density $\rho_{\nu}$ and the
magnetization norm $|\vec{M}|$.

###  3.7 Post-processing: K-mean heuristic

In our optimization problem we look for permanent magnetization directions
which may change continuously in space. In order to obtain designs which
comply with feasibility constraints, we here propose a post-processing step. A
K-mean heuristic [22] clustering method is applied. Here, we suggest adapting
it to create clusters of elements of similar magnetization direction.
Let us define a point $P=(p_{x},p_{y},p_{\beta_{M}})$ by its coordinates
$p_{x}$, $p_{y}$ in the 2D plane and its magnetization angle $p_{\beta_{M}}$.
We define a set of $k$ points $C_{1},\dots C_{k}$ where
$C_{j}=(c_{x,j},c_{y,j},-)$, which we will refer to as centroids, and which
are first randomly sampled in the 2D plane.
Let $P_{1},\dots,P_{N_{rot}}$ be the centroids of the triangles in the rotor
domain $\Omega_{rot}$. At each step of the algorithm, we associate each point
$P_{i}$ with the closest centroid $C_{j}\in S_{k}$ using a modified 3D
Euclidean distance $d_{\alpha}$. For two such points
$P_{i}=(p_{x,i},p_{y,i},p_{{\beta_{M}},i})$ and
$C_{j}=(c_{x,j},c_{y,j},c_{{\beta_{M}},j})$, this modified distance function
is defined as

$\displaystyle d_{\alpha}:$
$\displaystyle\mathbb{R}^{3}\times\mathbb{R}^{3}\rightarrow\mathbb{R}$ (3.34)
$\displaystyle(P_{i},C_{j})\mapsto\sqrt{\left(\frac{p_{x,i}-c_{x,j}}{N_{x}}\right)^{2}+\left(\frac{p_{y,i}-c_{y,j}}{N_{y}}\right)^{2}+\alpha\left(\frac{p_{{\beta_{M}},i}-c_{{\beta_{M}},j}}{2\pi}\right)^{2}}$

where $N_{x}$, $N_{y}$, $\alpha$ are three weighting constants which can be
used to tune the method.

By this procedure, we get $k$ point clusters $S_{1},\dots,S_{k}$ with
$S_{j}=\\{P_{i}:d_{\alpha}(P_{i},C_{j})<d_{\alpha}(P_{i},C_{l})$ for all
$l\neq j\\}$ for $j=1,\dots,k$. Thus, all points in cluster $S_{j}$ are closer
to point $C_{j}$ in terms of $d_{\alpha}$ than to any other centroid $C_{l}$.

The positions of the centroids $C_{1},\dots,C_{k}$ are then updated according
to the mean coordinates of the points associated with them,

$C_{j}\leftarrow(\mu_{x,j},\mu_{y,j},\mu_{\beta_{M,j}})\mbox{ where
}\mu_{x,j}=\frac{1}{\\#S_{j}}\sum_{P_{i}\in S_{j}}p_{x,i}$ (3.35)

and $\mu_{y,j},\mu_{\beta_{M,j}}$ defined analogously. Here, $\\#S_{j}$
denotes the cardinality of $S_{j}$. In the first iteration where
$c_{{\beta_{M}},j}$ is not defined, we use the modified 2D Euclidean distance
$d_{\alpha=0}$.

##  4 Application to the magnetostatics problem

All computations were conducted using the NGSolve [12, 32] framework with its
python interface. Three sets of results are presented in this section to
validate the proposed strategy. For the different optimization sets, the
starting points are given in Table 6.
The first design in Table 6 is used for the Iron-Air optimization. The results
for the material distribution are given in Table 7 for different volume
fractions $f_{\nu,f}$. These results are coherent with the literature of
synchro-reluctant actuators with distributed winding [34], with a phase angle
$\varphi=\frac{5\pi}{6}$ and one pair of poles. This phase angle corresponds
to the maximum torque of the machine. Results can also be compared with
results found in [25, 11] which are validated by means of JMAG .
For the multi-material topology optimization including magnets, the final
design obtained at $40\%$ of ferromagnetic material is used as a starting
point. Results for magnets distribution are presented in Table 8. As expected,
magnets are distributed on the air barriers domain. Several constraints on
magnets volume were chosen. One should note that the ferromagnetic
distribution is modified at the outer radius of the rotor. Without magnets,
the reluctance torque is equal to $1.0768$ [N.m]. Optimized torque and post-
processed values after k-mean clustering are given, the number of clusters is
fixed to $k=5$.
In order to validate the proposed strategy, a new phase angle is chosen. It is
equal to $\varphi=\frac{3\pi}{32}$. The objective is to find the optimal
materials distribution to maximize the torque at this phase angle. An unbiased
starting point is chosen, where homogeneous grey materials are set in the
rotor such as $\rho_{\nu}=0.5,\rho_{M_{x}}=0.5,\rho_{M_{y}}=0.5$. The optimal
rotor structure is given in Table 9. Again, we present optimized torques as
well as the torque values after post-processing, where we used $k=5$ clusters
for the K-mean clustering. Results on the torque value are comparable for the
previously obtained but at higher magnets volume, which is also coherent with
literature. On the other hand, the first optimal multi-material model has
larger reluctant torque due to optimal current supply in the q-axis.

Table 6: Starting designs |  |  Starting Design |
---|---|---|---
Rotor |  |  |
Table 7: Designs Iron-Air  |  |  Final filtered design |
---|---|---|---
Rotor |  |  |
$f_{\nu,f}$ |  10 % |  20 % |  40 %
$\bar{T}[N.m]$ |  0.20897 |  0.83351 |  1.1129
Table 8: Designs Magnet-Air-Iron  |  |  Final filtered designs |
---|---|---|---
Rotor |  |  |
$f_{\nu,mag}$ |  7.5 % |  15 % |  30 %
$\bar{T}[N.m]$ |  1.2830 |  1.4729 |  1.9000
$\bar{T}_{Kmeans}[N.m]$ |  1.2646 |  1.4469 |  1.8671
Table 9: Designs Iron-Air-Magnets |  |  Final Design |
---|---|---|---
Rotor |  |  |
$f_{v,mag}$ |  10% |  20% |  Not bounded
$\bar{T}\,[N.m]$ |  1.2710 |  1.4513 |  2.0484
$\bar{T}_{Kmeans}\,[N.m]$ |  1.1097 |  1.3863 |  1.6422

##  5 Conclusion and outlooks

In this study we proposed a novel multi-material interpolation method to
determine the optimal distribution of air, iron, and magnets for PMSynRM. The
interpolation takes into account magnetization amplitudes and direction, and a
post-processing clustering method is also suggested to homogenize magnets
direction for feasibility constraints. This study also investigated the use of
the four statics positions method to reduce global computation time. With
further work, an exhaustive comparison with other interpolation schemes will
be considered, and the extension to multiobjective optimization under
gradient-based methods.

##  Appendix A Nitsche-mortaring reformulation

Several methods exist to consider the movement in electrical motors, such as
the Moving Band (MB) technique described in [8], even when using high order
elements in the MB, this method remains less accurate than the mortar element
method [4]. A variant of the mortar method, the Nitsche method, is chosen to
take into account the rotation [16]. Let us define the adapted operator
$K_{\theta}^{NM}$ to fit the new formulation

$\begin{split}K_{\theta}^{NM}(\mathbf{X},u,\eta,v,\mu)=&\sum_{i\in\\{rot,stat\\}}\left(\int_{\Omega_{i}}\nu(\rho_{\nu}^{\theta},|\nabla
u_{i}|)\nabla u_{i}\cdot\nabla v_{i}\right)\\\
&-\int_{\Omega_{rot}}f_{\nu}(1-\rho_{\nu}^{\theta})\frac{M_{max}f_{M}(|\vec{M}|)}{|\vec{M}|}R_{\theta}\begin{bmatrix}-M_{y}^{\theta}\\\
M_{x}^{\theta}\end{bmatrix}\cdot\nabla v_{rot}\\\
&+\sum_{i\in\\{rot,stat\\}}\left(-\int_{\partial\Omega_{i}}\nu_{0}\frac{\partial
u_{i}}{\partial n}(v_{i}-\mu)-\int_{\partial\Omega_{i}}\nu_{0}\frac{\partial
v_{i}}{\partial n}(u_{i}-\eta)\right.)\\\
&\left.\phantom{+\sum_{\Omega\in\\{\Omega_{rot},\Omega_{stat}\\}}(}\;+\frac{\alpha
p^{2}}{h}\int_{\partial\Omega_{i}}(u_{i}-\eta)(v_{i}-\mu)\right),\end{split}$
(A.1)

where $(u,\eta,v,\mu)\in V\times W\times V\times W$ and $V,W$ are defined as

$\begin{cases}V=\\{v=(v_{rot},v_{stat})\in H^{1}(\Omega_{rot})\times
H^{1}(\Omega_{stat}),\,v=0\text{ on }\partial\Omega\\},\\\ W=\\{w\in
L^{2}(\partial\Omega_{rot}\cap\partial\Omega_{stat})\\}.\\\ \end{cases}$ (A.2)

Here, $\alpha>0$ is a stabilization parameter which we chose as $\alpha=160$,
$p=1$ denotes the polynomial degree of the finite element discretization and
$h$ the diameter of the largest element of the mesh. Moreover, recall the
implicit dependence of $\vec{M}=(M_{x}^{\theta},M_{y}^{\theta})^{\top}$ on the
density variables $\rho_{M_{x}^{\theta}},\rho_{M_{y}^{\theta}}$ (3.5), (3.14).

Hence the state equation (LABEL:eq:PDE_density) can be formulated as

$\mbox{Find }(u_{\theta},\eta_{\theta})\in V\times
W,\;K_{\theta}^{NM}(\mathbf{X}_{\rho},u_{\theta},\eta_{\theta},v,\mu)=\int_{\Omega_{stat}}j(\theta)\,vdx\quad\mbox{for
all }(v,\mu)\in V\times W.$ (A.3)

In a similar manner, the adjoint equation (3.22) can be redefined as

$\mbox{Find }(w_{\theta},\mu_{\theta})\in V\times W:\frac{\partial
K_{\theta}^{NM}}{\partial(u,\eta)}(\mathbf{X},u_{\theta},\eta_{\theta},w_{\theta},\mu_{\theta})(\hat{u},\hat{\eta})=-\frac{1}{4}\frac{\partial
T}{\partial u}(u_{\theta})(\hat{u})$ (A.4)

for all $(\hat{u},\hat{\eta})\in V\times W$ for each
$\theta\in\\{0,\frac{\pi}{12},\frac{\pi}{6},\frac{\pi}{4}\\}$.

The corresponding Lagrangian reads

$\tilde{\mathcal{L}}(\mathbf{X},\underline{u},\underline{\eta},\underline{v},\underline{\mu}):=L(\mathbf{X},\underline{u})+\sum_{\theta\in\\{0,\frac{\pi}{12},\frac{\pi}{6},\frac{\pi}{4}\\}}K_{\theta}^{NM}(\mathbf{X},u_{\theta},\eta_{\theta},v_{\theta},\mu_{\theta})-\int_{\Omega_{stat}}j(\theta)v_{\theta},$
(A.5)

and its gradient is given by

$\begin{split}\frac{\partial\mathcal{\tilde{L}}}{\partial\mathbf{X}}(\mathbf{X},\underline{u},\underline{\eta},\underline{v},\underline{\mu})=&\left(\sum_{\theta\in\\{0,\frac{\pi}{12},\frac{\pi}{6},\frac{\pi}{4}\\}}\frac{\partial
K^{NM}_{\theta}}{\partial\mathbf{X}}(\mathbf{X},u_{\theta},\eta_{\theta},v_{\theta},\mu_{\theta})-\frac{1}{4}\frac{\partial
T_{\theta}}{\partial\mathbf{X}}\right)\\\
&+\frac{\partial\psi(h_{v,f}(\mathbf{X}),\gamma_{f},\mu)}{\partial\mathbf{X}}+\frac{\partial\psi(h_{v,mag}(\mathbf{X}),\gamma_{mag},\mu)}{\partial\mathbf{X}}.\end{split}$
(A.6)

## References

  * [1] P. Akiki, M. Hage-Hassan, M. Bensetti, J.-C. Vannier, D. Prieto, and M. McClelland.  Axial ferrite-magnet-assisted synchronous reluctance motor.  In 2018 XIII International Conference on Electrical Machines (ICEM), pages 583–589, 2018.
  * [2] S. Amstutz, C. Dapogny, and À. Ferrer.  A consistent relaxation of optimal design problems for coupling shape and topological derivatives.  Numerische Mathematik, 140(1):35–94, 3 2018.
  * [3] S. Amstutz and P. Gangl.  Topological derivative for the nonlinear magnetostatic problem.  Electron. Trans. Numer. Anal., 51:169–218, 2019.
  * [4] O. Antunes, J. Bastos, N. Sadowski, A. Razek, L. Santandrea, F. Bouillault, and F. Rapetti.  Using hierarchic interpolation with mortar element method for electrical machines analysis.  IEEE Transactions on Magnetics, 41(5):1472–1475, 2005.
  * [5] N. Bianchi, L. Alberti, M. Popescu, and T. Miller.  MMF Harmonics Effect on the Embedded FE-Analytical Computation of PM Motors.  In Conference Record - IAS Annual Meeting (IEEE Industry Applications Society), pages 1544 – 1551, 2007.
  * [6] J. S. Choi and J. Yoo.  Optimal design method for magnetization directions of a permanent magnet array.  Journal of Magnetism and Magnetic Materials, 322(15):2145–2151, 2010.
  * [7] J. S. Choi, J. Yoo, S. Nishiwaki, and K. Izui.  Optimization of Magnetization Directions in a 3-D Magnetic Structure.  IEEE Transactions on Magnetics, 46(6):1603–1606, 2010.
  * [8] B. Davat, Z. Ren, and M. Lajoie-Mazenc.  The movement in field modeling.  IEEE Transactions on Magnetics, 21(6):2296–2298, 1985.
  * [9] A. Ferrer.  Simp-all: A generalized simp method based on the topological derivative concept.  International Journal for Numerical Methods in Engineering, 120(3):361–381, 2019.
  * [10] C. Fong.  Analytical methods for squaring the disc, 2019.
  * [11] P. Gangl, S. Köthe, C. Mellak, A. Cesarano, and A. Mütze.  Multi-objective free-form shape optimization of a synchronous reluctance machine.  arXiv:2010.10117 [cs, math], 10 2020.  arXiv: 2010.10117.
  * [12] P. Gangl, K. Sturm, M. Neunteufel, and J. Schöberl.  Fully and Semi-Automated Shape Differentiation in NGSolve.  arXiv:2004.06783 [math], 10 2020.  arXiv: 2004.06783.
  * [13] H. Garcke, C. Hecht, M. Hinze, and C. Kahle.  Numerical Approximation of Phase Field Based Shape and Topology Optimization for Fluids.  SIAM J. Sci. Comput., 37(4):A1846–A1871, 1 2015.
  * [14] F. Guo, M. Salameh, M. Krishnamurthy, and I. P. Brown.  Multimaterial Magneto-Structural Topology Optimization of Wound Field Synchronous Machine Rotors.  IEEE Transactions on Industry Applications, 56(4):3656–3667, 2020\.
  * [15] H. Heidari, A. Rassõlkin, A. Kallaste, T. Vaimann, E. Andriushchenko, A. Belahcen, and D. V. Lukichev.  A Review of Synchronous Reluctance Motor-Drive Advancements.  Sustainability, 13(2):729, 1 2021.
  * [16] K. Hollaus, D. Feldengut, J. Schoeberl, M. Wabro, and D. Omeragic.  Nitsche-type Mortaring for Maxwell’s Equations.  PIERS 2010 Cambridge - Progress in Electromagnetics Research Symposium, Proceedings, 2010.
  * [17] T. Ishikawa, P. Xie, and N. Kurita.  Topology Optimization of Rotor Structure in Permanent Magnet Synchronous Motors Considering Ease of Manufacturing.  IEEJ Journal IA, 4(4):469–475, 2015.
  * [18] J. Kim, K. H. Sun, W. Kim, and J. E. Kim.  Magnetic torque maximization in a camera shutter module by the topology optimization.  J Mech Sci Technol, 24(12):2511–2517, 12 2010.
  * [19] B. S. Lazarov and O. Sigmund.  Filters in topology optimization based on Helmholtz-type differential equations.  Int. J. Numer. Meth. Engng., 86(6):765–781, 5 2011.
  * [20] J. Lee, E. M. Dede, and T. Nomura.  Simultaneous Design Optimization of Permanent Magnet, Coils, and Ferromagnetic Material in Actuators.  IEEE Transactions on Magnetics, 47(12):4712–4716, 2011.
  * [21] D. Lukáš.  An Integration of Optimal Topology and Shape Design for Magnetostatics.  In A. M. Anile, G. Alì, and G. Mascali, editors, Scientific Computing in Electrical Engineering, pages 227–232, Berlin, Heidelberg, 2006. Springer Berlin Heidelberg.
  * [22] J. MacQueen and others.  Some methods for classification and analysis of multivariate observations.  In Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, volume 1, pages 281–297. Oakland, CA, USA, 1967\.  Issue: 14.
  * [23] S. S. Maroufian and P. Pillay.  PM assisted synchronous reluctance machine design using AlNiCo magnets.  In 2017 IEEE International Electric Machines and Drives Conference (IEMDC), pages 1–6, Miami, FL, USA, 5 2017. IEEE.
  * [24] A. Marrocco.  Analyse numérique de problèmes d’électrotechnique.  Ann. Sc. Math. Québec, 1(2):271–296, 1977.
  * [25] C. Mellak, K. Krischan, and A. Muetze.  Synchronous Reluctance Machines as Drives for Rotary Anode X-Ray Tubes-A Feasibility Study.  In 2018 XIII International Conference on Electrical Machines (ICEM), pages 2613–2618, Alexandroupoli, 9 2018. IEEE.
  * [26] P. Monk and others.  Finite element methods for Maxwell’s equations.  Oxford University Press, 2003.
  * [27] J. Nocedal and S. J. Wright, editors.  Penalty, Barrier, and Augmented Lagrangian Methods, pages 488–525.  Springer New York, New York, NY, 1999.
  * [28] A. S. Nunes, L. Daniel, M. Hage-Hassan, and M. Domenjoud.  Modeling of the magnetic behavior of permanent magnets including ageing effects.  Journal of Magnetism and Magnetic Materials, 512:166930, 10 2020\.
  * [29] M. Risticevic, D. Iles, and A. Moeckel.  Design of an interior permanent magnet synchronous motor supported by the topology optimization algorithm.  In 2016 International Symposium on Power Electronics, Electrical Drives, Automation and Motion (SPEEDAM), pages 221–225, Capri, Italy, 6 2016. IEEE.
  * [30] N. Sadowski, Y. Lefevre, M. Lajoie-Mazenc, and J. Cros.  Finite element torque calculation in electrical machines while considering the movement.  IEEE Trans. Magn., 28(2):1410–1413, 3 1992.
  * [31] S. Sanogo and F. Messine.  Topology optimization in electromagnetism using SIMP method: Issues of material interpolation schemes.  COMPEL, 37(6):2138–2157, 11 2018.
  * [32] J. Schöberl.  C++11 Implementation of Finite Elements in NGSolve.  Technical report, Institute for Analysis and Scientific Computing, 9 2014\.
  * [33] O. Sigmund.  A 99 line topology optimization code written in Matlab.  Structural and Multidisciplinary Optimization, 21, 2001.
  * [34] A. Vagati, G. Franceschini, I. Marongiu, and G. Troglia.  Design criteria of high performance synchronous reluctance motors.  In Conference Record of the 1992 IEEE Industry Applications Society Annual Meeting, pages 66–73 vol.1, 1992.
  * [35] F. Wang, B. S. Lazarov, and O. Sigmund.  On projection methods, convergence and robust formulations in topology optimization.  Structural and Multidisciplinary Optimization, 43(6):767–784, 6 2011.
  * [36] S. Wang, D. Youn, H. Moon, and J. Kang.  Topology optimization of electromagnetic systems considering magnetization direction.  Magnetics, IEEE Transactions on, 41:1808 – 1811, 2005.
  * [37] W. Zuo and K. Saitou.  Multi-material topology optimization using ordered SIMP interpolation.  Struct Multidisc Optim, 55(2):477–491, 2 2017.


# Circle action of the punctured mapping class group and cross homomorphism

Lei Chen  Lei Chen
Department of Mathematics
University of Maryland
4176 Campus Drive
College Park, MD 20742, USA
<EMAIL_ADDRESS>
###### Abstract.

In the following short note, we give a new geometric interpretation of the
generator of the infinite cyclic group
$H^{1}(\text{Mod}(S_{g,1});H^{1}(S_{g};\mathbb{Z}))$ (this computation is
proved by Morita). There are several construction of this class given by
Earle, Morita, Trapp and Furuta. The construction we give here uses the action
of $\text{Mod}(S_{g,1})$ on the circle and its rotation numbers. We suspect
that our construction is the same as the construction by Furuta and Trapp
using winding numbers and provide half of the proof.

##  1\. The construction and result

Let $S_{g}$ be a genus $g$ surface and let $p$ be a point on $S_{g}$. When
$g>1$, the universal cover of $S_{g}$ is the hyperbolic plane
$\mathbb{H}^{2}$. Pick a lift $\tilde{p}$ of $p$ in $\mathbb{H}^{2}$, any
diffeomorphism $\phi$ of $S_{g}$ fixing $p$ can be lifted to a diffeomorphism
$\tilde{\phi}$ of $\mathbb{H}^{2}$ fixing $\tilde{p}$. Furthermore, the lift
$\tilde{\phi}$ extends to the boundary $\partial\mathbb{H}^{2}\cong S^{1}$.
The boundary action $\partial(\tilde{\phi})$ does not depend on the isotopy
type of $\phi$. The above construction describes the following Gromov boundary
action

$G:\operatorname{Mod}(S_{g,1})\to\operatorname{Homeo}^{+}(S^{1}).$

For a proof of the above description, see [FM12, Chapter 5.5]. Note that
Mann–Wolff [MW20] prove that any action of $\operatorname{Mod}(S_{g,1})$ on
the circle is either trivial or semiconjugate to $G$.

Let $\widetilde{\operatorname{Homeo}^{+}(S^{1})}$ be the homeomorphism of
$\mathbb{R}$ that commutes with translation by one map on $\mathbb{R}$. Then
we have the following short exact sequence.

(1)
$1\to\mathbb{Z}\to\widetilde{\operatorname{Homeo}^{+}(S^{1})}\xrightarrow{q}\operatorname{Homeo}^{+}(S^{1})\to
1$

Let
$\operatorname{trans}:\widetilde{\operatorname{Homeo}^{+}(S^{1})}\to\mathbb{R}$
be the translation number, which is defined as the following limit, which
always exists

$\operatorname{trans}(f)=\lim_{n\to\infty}\frac{f^{n}(0)}{n}.$

Let $P:\pi_{1}(S_{g})\to\operatorname{Mod}(S_{g,1})$ be the point-pushing
homomorphism. Let $e:F_{2g}\to\pi_{1}(S_{g})$ be the natural homomorphism
where $F_{2g}$ is generated by $a_{1},b_{1},...,a_{g},b_{g}$ and the only
relation of $\pi_{1}(S_{g})$ is given by $c:=[a_{1},b_{1}]...[a_{g},b_{g}]=1$.
Even though homomorphism $G\circ P$ does not have a lift under the map q to
$\widetilde{\operatorname{Homeo}^{+}(S^{1})}$ in the exact sequence (1), the
free group $G\circ P\circ e$ can be lifted to
$\widetilde{\operatorname{Homeo}^{+}(S^{1})}$. Let

$\tilde{G}:F_{2g}\to\widetilde{\operatorname{Homeo}^{+}(S^{1})}$

be a lift of $G\circ P\circ e$ (which is not unique).

For a group $H$, an outer-automorphism of $H$ is an automorphism of $H$ up to
a conjugation of an element in $H$. By Dehn–Nielsen–Baer Theorem of
$\operatorname{Mod}(S_{g,1})$ (see, e.g., [FM12, Chapter 8]), we know that
$\operatorname{Mod}(S_{g,1})\cong\operatorname{Out}^{*}(F_{2g})$, where
$\operatorname{Out}^{*}(F_{2g})$ denotes the outer-automorphism group of
$F_{2g}$ consisting of elements $f:F_{2g}\to F_{2g}$ such that $f(c)$ is
conjugate to $c$. By composing with a conjugation from $F_{2g}$, any element
in $\operatorname{Out}^{*}(F_{2g})$ has a representative $f:F_{2g}\to F_{2g}$
such that $f(c)=c$ (such $f$ is still not unique but differed by an element in
the centralizer of $c$, which is a power of $c$).

We now define a map

$R:\operatorname{Mod}(S_{g,1})\to\text{Map}(F_{2g},\mathbb{Z}).$

Let $\phi\in\operatorname{Mod}(S_{g,1})$ and $f:F_{2g}\to F_{2g}$ be a
representative of $\phi$ such that $f(c)=c$. Then the definition of $R$ is
given by the following formula

$R(\phi)(\gamma)=\operatorname{trans}(\tilde{G}(f(\gamma)))-\operatorname{trans}(\tilde{G}(\gamma))$

For a group $H$ and an $H$-module $M$, we call a map $\rho:H\to M$ a cross
homomorphism if

$\rho(gh)=\rho(g)+g(\rho(h))$

The main result of this paper is the following.

######  Theorem 1.1.

The map $R$ is well-defined, has image in $\text{Hom}(F_{2g},\mathbb{Z})\cong
H^{1}(\pi_{1}(S_{g});\mathbb{Z})$ and $R$ is a cross-homomorphism.
Furthermore, as an element in the cohomology class, $[R]$ is a generator of
$H^{1}(\operatorname{Mod}(S_{g,1});H_{1}(\pi_{1}(S_{g});\mathbb{Z})$, which is
equal to the cross ratio defined by Morita [Mor89].

Notice that Earle [Ear78], Morita [Mor89], Trapp [Tra92] and Furuta [Mor97]
(described by Morita) give various constructions of the generator of the group

$H^{1}(\operatorname{Mod}(S_{g,1});H_{1}(\pi_{1}(S_{g});\mathbb{Z}))\cong\mathbb{Z}.$

Trapp and Futura use the unit tangent bundle and winding numbers. Earle uses
Moduli spaces and Morita gives a combinatorial construction. All of the above
definition seem very different. It is an interesting question to consider how
those construction interact with each other. Kuno [Kun09] gave an
understanding of the difference between Earle construction and Morita
construction. This cross homomorphism is also related to Johnson homomorphism
as discussed in Morita [Mor93].

The second result is to relate our construction with that of Trapp [Tra92] and
Furuta [Mor97] (described by Morita). Let $D\subset S_{g}$ be an open disk.
Let $X$ be a nowhere zero vector field on $S_{g}-D$. Then there is a winding
number map

$\omega_{X}:\pi_{1}(S_{g}-D)=F_{2g}\to\mathbb{Z}$

given by how $X$ rotates along a smooth representative of a curve in
$\pi_{1}(S_{g}-D)$.

For any map $h:F_{2g}\to\mathbb{Z}$, we can define its defect

$D(h):F_{2g}\times F_{2g}\to\mathbb{Z}$

as $D(h)(a,b)=h(ab)-h(a)-h(b)$, which measures how far $h$ is from a
homomorphism. We ask the following question.

######  Problem 1.2.

Do we have $D(\omega_{X})=D(\operatorname{trans}\circ\tilde{G})$? Moreover,
does there exist a nowhere zero vector field $X$ on $S_{g}-D$ such that
$\omega_{X}=\operatorname{trans}\circ\tilde{G}$?

The positive answer of the first question implies the positive answer of the
second question. We also prove the following.

######  Theorem 1.3.

For $\alpha,\beta$ which generate a free group of order $2$ corresponds to a
cover $U$ of $S_{g,1}$ that is a once-punctured torus, we have that

$D(\omega_{X})(\alpha,\beta)=D(\operatorname{trans}\circ\tilde{G})(\alpha,\beta)=0.$

Acknowledgement. We thank Bena Tshishiku for introducing me the paper of
Morita.

##  2\. The proof of Theorem 1.1

Step 1: $R$ is well-defined, has image in $\text{Hom}(F_{2g},\mathbb{Z})\cong
H^{1}(\pi_{1}(S_{g});\mathbb{Z})$.

Let $\phi\in\operatorname{Mod}(S_{g,1})$. To prove that $R(\phi)$ is well-
defined, we need to prove that it does not depend on the choice of
$f:F_{2g}\to F_{2g}$ representing $\phi\in\operatorname{Out}^{*}(F_{2g})$ such
that $f(c)=c$. This comes from the fact that a different choice is given by
composing $f$ with a conjugation by a power of $c$ where $\tilde{G}(c)$ is a
translation by integers. Conjugation by a translation of integers has no
effect on the translation number since translation by integers commute with
$\widetilde{\operatorname{Homeo}^{+}(S^{1})}$.

We now show that $R$ has image in $\text{Hom}(F_{2g},\mathbb{Z})\cong
H^{1}(\pi_{1}(S_{g});\mathbb{Z})$.

######  Claim 2.1.

Let $\phi\in\operatorname{Mod}(S_{g,1})$ and $\alpha,\beta\in F_{2g}$, then

$R(\phi)(\alpha\beta)=R(\phi)(\alpha)+R(\phi)(\beta)$

###### Proof.

Let $f:F_{2g}\to F_{2g}$ be a representative of
$\phi\in\operatorname{Out}^{*}(F_{2g})$ such that $f(c)=c$.

We have the following canonical Euler cocycle

$\tau:\operatorname{Homeo}^{+}(S^{1})\times\operatorname{Homeo}^{+}(S^{1})\to\mathbb{Z}$

given by

$\tau(f,g):=\operatorname{trans}(\tilde{f}\tilde{g})-\operatorname{trans}(\tilde{f})-\operatorname{trans}(\tilde{g}),$

where $\tilde{f}$ and $\tilde{g}$ are lifts of $f,g$ in
$\widetilde{\operatorname{Homeo}}^{+}(S^{1}))$. The value of $\tau$ is in
$\\{-1,0,1\\}$.

By Matsumoto [Mat87, Theorem 3.3], $\tau$ is the same as a geometric cocyle

$\theta:\pi_{1}(S_{g})\times\pi_{1}(S_{g})\to\mathbb{Z}$

when restricted to the Fuchsian representation
$\pi_{1}(S_{g})\to\operatorname{Homeo}^{+}(S^{1})$. The cocycle $\theta$ is
given by the following rule: Suppose that the subgroup $F$ generated by
$\alpha,\beta$ is a free group on two generator, consider the covering $U\to
S_{g}$ associated with $F$, $U$ is either a punctured torus or a pair of
pants. If $U$ is a pair of pants and that $\alpha$ and $\beta$ are represented
by two boundary curves which match (resp. oppose) the orientation of $U$, we
define $\theta(\alpha,\beta)=-1(resp.-1)$. In all the other cases, define
$\theta(\alpha,\beta)=0$.

The formula $\theta(\alpha,\beta)$ is equivariant under the action of the
mapping class group $\operatorname{Mod}(S_{g,1})$. Thus from the
identification of $\theta$ and $\tau$, we obtain that

$\tau(\phi(\alpha),\phi(\beta))=\tau(\alpha,\beta).$

This implies that

(2)
$\begin{split}&\operatorname{trans}(\tilde{G}(\alpha)\tilde{G}(\beta))-\operatorname{trans}(\tilde{G}(\alpha))-\operatorname{trans}(\tilde{G}(\beta))\\\
=&\operatorname{trans}(\tilde{G}(f(\alpha))\tilde{G}(f(\beta)))-\operatorname{trans}(\tilde{G}(f(\alpha)))-\operatorname{trans}(\tilde{G}(f(\beta)))\end{split}$

(3)
$\begin{split}&\operatorname{trans}(\tilde{G}(f(\alpha))\tilde{G}(f(\beta)))-\operatorname{trans}(\tilde{G}(\alpha)\tilde{G}(\beta))\\\
=&\operatorname{trans}(\tilde{G}(f(\alpha)))-\operatorname{trans}(\tilde{G}(\alpha))+\operatorname{trans}(\tilde{G}(f(\beta)))-\operatorname{trans}(\tilde{G}(\beta))\end{split}$

Then

$R(\phi)(\alpha\beta)=R(\phi)(\alpha)+R(\phi)(\beta).\qed$

Step 2: $R$ is a cross homomorphism. We need to show that

$R(\phi\eta)=R(\phi)+\phi^{*}(R(\eta))$

Let $f,h:F_{2g}\to F_{2g}$ be representative of $\phi,\eta$ such that
$f(c)=h(c)=c$. Then $h\circ f$ is a representative of $\phi\eta$ such that
$h\circ f(c)=c$. Then

(4) $\begin{split}&R(\phi\eta)(\alpha)\\\ =&\operatorname{trans}(h\circ
f(\alpha))-\operatorname{trans}(\alpha)\\\
=&\operatorname{trans}(h(f(\alpha)))-\operatorname{trans}(f(\alpha))+\operatorname{trans}(f(\alpha))-\operatorname{trans}(\alpha)\\\
=&R(\eta)(f(\alpha))+R(\phi)(\alpha)\\\
=&\phi^{*}R(\eta)(\alpha)+R(\phi)(\alpha)\end{split}$

This proves what we need.

Step 3: $R$ is the generator of $\operatorname{Mod}(S_{g,1})\to
H^{1}(S_{g};\mathbb{Z})$. For this part, we only need to check the restriction
on the point-pushing subgroup
$P:\pi_{1}(S_{g})\to\operatorname{Mod}(S_{g,1})$. This gives a homomorphism

$R|_{P}:\pi_{1}(S_{g})\to H^{1}(S_{g};\mathbb{Z}),$

since the point-pushing subgroup acts trivially on $H^{1}(S_{g};\mathbb{Z})$.
The homomorphism $R|_{P}$ factors through the abelianization
$H_{1}(S_{g};\mathbb{Z})$ of $\pi_{1}(S_{g})$. Let
$i:\pi_{1}(S_{g})\times\pi_{1}(S_{g})\to\mathbb{Z}$ be the algebraic
intersection number, we have the following.

######  Claim 2.2.

We have the following formula

$R|_{P}(a)(b)=(2-2g)i(a,b).$

###### Proof.

We will do this by a calculation. Let $a_{1},b_{1},...,a_{g},b_{g}$ be the
standard generating set of $\pi_{1}(S_{g})$. The point-pushing $P(a_{1})$ has
the following representative $f:F_{2g}\to F_{2g}$ such that

(5) $a_{1}\to a_{1},b_{1}\to a_{1}^{-1}cb_{1}a_{1},a_{j}\to
a_{1}^{-1}ca_{j}c^{-1}a_{1},b_{j}\to a_{1}^{-1}cb_{j}c^{-1}a_{1}.$

We obtain the above representative by first find some representative using the
exact mapping class, then we compose with a conjugation to make sure that
$f(c)=c$.

Now we can easily obtain that $R(P(a_{1}))(a_{1})=0$,
$R(P(a_{1}))(b_{1})=2-2g$ and $R(P(a_{1}))(a_{j})=0,R(P(a_{1}))(b_{j})=0$ for
$1<j\leq g$. Thus we obtain that

$R(P(a_{1}))(b)=(2-2g)i(a_{1},b).$

Let $c$ be an element in $H_{1}(S_{g};\mathbb{Z})$ represented as a non-
separating simple closed curve. Then there is an element
$\phi\in\operatorname{Mod}(S_{g,1})$ such that $\phi(a_{1})=c$. Then we have
that

$R(\phi^{-1}P(a_{1})\phi)=R(P(c))$

By cross homomorphism of $R$, we know that

$R(\phi^{-1}P(a_{1})\phi)=R(\phi^{-1}P(a_{1}))+(\phi^{-1}P(a_{1}))^{*}R(\phi)=R(\phi^{-1})+(\phi^{*})^{-1}(P(a_{1}))+(\phi^{*})^{-1}(R(\phi))$

Since $R(\phi^{-1}\phi)=R(id)=0$, we know that
$R(\phi^{-1})+(\phi^{*})^{-1}(R(\phi^{-1}))=0$.

Thus we obtain that $R(P(c))=(\phi^{*})^{-1}(R(P(a_{1})))$, then we have

$\begin{split}&R(P(c))(x)=(\phi^{*})^{-1}(R(P(a_{1})))(x)=R(P(a_{1}))((\phi_{*})^{-1}(x))\\\
=&i(a_{1},(\phi_{*})^{-1}(x))=i(\phi_{*}(a_{1}),x)=i(c,x)\qed\end{split}$

Then by Morita [Mor89, Section 4], we know that $R$ as a cohomology class in
$H^{1}(\operatorname{Mod}(S_{g,1}),\mathbb{Z})$ is the same as Morita’s class.

##  3\. Open questions

We first identify the defects of two winding numbers. The following claim
implies that if the first question in Problem 1.2 has a positive answer, so
does the second question.

######  Claim 3.1.

Let $X,Y$ be two different nowhere zero vector fields on $S_{g}-D$. Then
$D(\omega_{X})=D(\omega_{Y})$. The defect
$D(\operatorname{trans}\circ\tilde{G})$ also does not depend on the choice of
the lift $\tilde{G}$.

###### Proof.

We have that $\omega_{X}-\omega_{Y}:F_{2g}\to\mathbb{Z}$ is a homomorphism,
which is proved by Trapp [Tra92, Section 1.3]. Thus we have that

$\omega_{X}(ab)-\omega_{Y}(ab)=\omega_{X}(a)-\omega_{Y}(a)+\omega_{X}(b)-\omega_{Y}(b).$

A rearrangement gives us the result. The proof for the defect
$D(\operatorname{trans}\circ\tilde{G})$ is similar. ∎

We now prove Theorem 1.3.

###### Proof.

We now prove that if $\alpha,\beta$ generate a free group of order $2$ such
that the cover determined by $\langle\alpha,\beta\rangle$ is a once-punctured
torus, then

$D(\omega_{X})(\alpha,\beta)=0.$

We can choose $X$ to be a non-where zero vector field on the whole torus. Then
since every loop can be straightened on the torus, we know that the winding
number $\omega_{X}(\gamma)=0$ for any $\gamma$. The fact that
$D(\operatorname{trans}\circ\tilde{G})(\alpha,\beta)=0$ is given by Matsumoto
[Mat87, Theorem 3.3]. ∎

In Morita’s construction, he constructed a map $f:F_{2g}\to\mathbb{Z}$ such
that

$D(f)(a,b)=f(ab)-f(a)-f(b)=\langle a,b\rangle,$

where $\langle a,b\rangle$ denotes the algebraic intersection number. Then he
define the cross homomorphism the same way as the following

$C_{f}(\phi)(a)=f(\phi(a))-f(a).$

Whether or not $C_{f}$ is a cross homomorphism or not only depends on $D(f)$
that $D(f)$ has to be invariant under $\operatorname{Mod}(S_{g,1})$ and
satisfies an association equation. We wonder about the following.

######  Question 3.2.

Other than $f$, $\omega_{X}$ and $\operatorname{trans}\circ\tilde{G}$
(possibly the same as $\omega_{X}$), do we have other functions
$f:F_{2g}\to\mathbb{Z}$ such that $C_{f}$ is a cross homomorphism?

We also would like to consider the difference between $C_{f}$ and $R$ as in
Kuno [Kun09].

######  Question 3.3.

What’s the coboundary of $C_{f}-R$?

## References

  * [Ear78] Clifford J. Earle.  Families of Riemann surfaces and Jacobi varieties.  Ann. of Math. (2), 107(2):255–286, 1978.
  * [FM12] Benson Farb and Dan Margalit.  A primer on mapping class groups, volume 49 of Princeton Mathematical Series.  Princeton University Press, Princeton, NJ, 2012.
  * [Kun09] Yusuke Kuno.  A combinatorial formula for Earle’s twisted 1-cocycle on the mapping class group $\mathcal{M}_{g,\ast}$.  Math. Proc. Cambridge Philos. Soc., 146(1):109–118, 2009.
  * [Mat87] Shigenori Matsumoto.  Some remarks on foliated $S^{1}$ bundles.  Invent. Math., 90(2):343–358, 1987.
  * [Mor89] Shigeyuki Morita.  Families of Jacobian manifolds and characteristic classes of surface bundles. I.  Ann. Inst. Fourier (Grenoble), 39(3):777–810, 1989.
  * [Mor93] Shigeyuki Morita.  The extension of Johnson’s homomorphism from the Torelli group to the mapping class group.  Invent. Math., 111(1):197–224, 1993.
  * [Mor97] Shigeyuki Morita.  Casson invariant, signature defect of framed manifolds and the secondary characteristic classes of surface bundles.  J. Differential Geom., 47(3):560–599, 1997.
  * [MW20] Kathryn Mann and Maxime Wolff.  Rigidity of mapping class group actions on $S^{1}$.  Geom. Topol., 24(3):1211–1223, 2020.
  * [Tra92] Rolland Trapp.  A linear representation of the mapping class group ${\mathcal{M}}$ and the theory of winding numbers.  Topology Appl., 43(1):47–64, 1992.