# Prime Counting Function identity

Suyash Garg  Department of Mathematics
BITS-Pilani Hyderabad Campus
Hyderabad, India
<EMAIL_ADDRESS>
###### Abstract

In this paper, a new formula for $\pi^{(2)}(N)$ is formulated, it is a
function that counts the number of semi-primes not exceeding a given number N.
A semi-prime is a natural number that is the product of exactly two prime
numbers, the two primes in the product may equal each other. Semi-prime
numbers are also a case of almost primes. Since a formula for this is already
known, a new identity that uses the prime counting function is created by
equating the two functions.

###### Index Terms:

Semi-primes, Semi-prime counting function, Prime counting function, k-almost
primes

##  I Introduction

A semi-prime, also called a 2-almost prime, biprime[1][2], or $pq$-number, is
a composite number that is the product of two (possibly equal) primes. The
first few are 4,6,9,10,14,15,21,22…(OEIS A001358). The first few semi-primes
whose factors are distinct are 6, 10, 14, 15, 21, 22, 26…(OEIS A006881). A
formula for the number of semi-primes less than or equal to $N$ is given by
[3]

$\pi^{(2)}(N)=\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})-k+1],$ (1)

here, $p_{k}$ is the $k^{th}$ prime, and $\pi(N)$ is the prime counting
function, which is a function that gives the number of prime numbers not
exceeding the given number N[4].
In this paper, a new way to calculate the number of semi-primes not exceeding
$N$, is explored.

##  II Semi-Prime Counting function

Lemma 1. The number of numbers of the form $2p$ not exceeding a given number
$N$, is $\pi(N/2)$, where $p$ is any prime number and $\pi(x)$ is the prime
counting function.

Proof. The number of primes such that $2p\leq N$ is to be calculated. Just by
rearranging, it can be shown that $p\leq N/2$. This means p can be any prime
smaller than $N/2$, which equals $\pi(N/2)$.

Corollary 1. The number of numbers of the form $pq$ not exceeding a given
number $N$, is $\pi(N/q)$, where $p$ is a given constant prime number, $q$ can
be any prime number, and $\pi(x)$ is the prime counting.

Proof. If $p$ is equal to $2$ then, It would be the same as Lemma 1, So, for
proving the corollary, it would be sufficient to replace $2$ with a given
prime $p$ in Lemma 1.

Lemma 2. The number of pairs of the form $(p,q)$ where $p$ and $q$ can be any
prime such that $pq\leq N$ for a given $N$ is

$\sum_{k=1}^{\pi(N/2)}[\pi(N/p_{k})]$ (2)

Proof. In Corollary 1, each prime $p_{k}$ less than N can be used as the
constant prime $p$,
for $p=2$, number of pairs of the form $(2,q)=\pi(N/2)$,
for $p=3$, number of pairs of the form $(3,q)=\pi(N/3)$,
for $p=p_{k}$, number of pairs of the form $(p_{k},q)=\pi(N/p_{k})$,
Now, if all of these were to be added, the summation acquired would be,

$\sum_{k=1}^{\pi(N)}[\pi(N/p_{k})]$

Since, all the primes till $N$ were considered for $p$, and $q$ can be any
prime, therefore the above formula suffices as a formula for the number of
pairs of the form $(p,q)$ such that $pq\leq N$.
Also Since, for all primes $p>N/2$, $N/p<2$ , which means $\pi(N/p)=0$,
therefore the upper limit of the summation can be reduced to $\pi(N/2)$. So,
the new formula will be

$\sum_{k=1}^{\pi(N/2)}[\pi(N/p_{k})]$

Note. In Lemma 2, for all the pairs $(p,q)$, where $p$ and $q$ are distinct
primes, if order doesn’t matter, then each pair is counted twice.

Lemma 3. The number of primes such that $p^{2}\leq N$, where $p$ can be any
prime is $\pi(\sqrt{N})$.

Proof. Since $p>0$ and $N>0$,we can get the inequality $p\leq\sqrt{N}$.The
number of primes such that $p^{2}\leq N$ is equal to the number of primes such
that $p\leq\sqrt{N}$, which is equal to $\pi(\sqrt{N})$.

Theorem 1.

$\pi^{(2)}(N)=((\sum_{k=1}^{\pi((N/2))}[\pi(N/p_{k})])+\pi(\sqrt{N}))/2,$ (3)

is the function that gives the number of semi-primes not exceeding a given
number $N$, i.e. a semi-prime counting function.

Proof. Lemma 2 takes into account all the number of pairs $(p,q)$, such that
$p$ and $q$ can be any primes, and $(p,q)$ and $(q,p)$ are considered as $2$
distinct pairs. So, It is trivial that, all the pairs where $p$ and $q$ are
distinct will be counted twice and pairs of the form $(p,p)$ will be counted
once. The number of pairs of the form $(p,p)$, is equivalent to the number of
primes such that $p^{2}\leq N$, and Lemma 3 helps get the answer to this
question. Now, if this is added to (2), the new formula would be

$\sum_{k=1}^{\pi(N/2)}[\pi(N/p_{k})]+\pi(\sqrt{N})$ (4)

In (4) each and every pair is counted twice. so if (4) is halved, the number
of pairs of the form $(p,q)$ such that $p$ and $q$ can be any primes, without
considering order, can be obtained. This formula will be equivalent to the
formula for semi-primes less than or equal to $N$, as, if we multiply the
numbers $p$ and $q$ in all the given pairs, we get a list of all the semi-
prime less than or equal to $N$, and since we have counted each such pair only
once, there would be no repetitions in this list. So, counting the number of
elements in this list is equivalent to counting the number of pairs as
described above. Therefore, the final formula for the semi-prime counting
function is,

$((\sum_{k=1}^{\pi((N/2))}[\pi(N/p_{k})])+\pi(\sqrt{N}))/2$

##  III Example

Now, there are two formulas for semi-prime counting function.

Example 1.
Let’s take $n=25$.
Semi-primes not exceeding $25$ are
$4,6,9,10,14,15,21,22,25$
so,

$\pi^{(2)}(25)=9.$ (5)

Using (1)

$\pi^{(2)}(25)=\sum_{k=1}^{\pi(5)}[\pi(25/p_{k})-k+1]$
$\pi(5)=3,\pi(8.33)=4,\pi(12.5)=5$ $\pi^{(2)}(25)=(5-1+1)+(4-2+1)+(3-3+1)$
$\pi^{(2)}(25)=5+3+1=9$ (6)

Using (3)

$\pi^{(2)}(25)=((\sum_{k=1}^{\pi((12.5))}[\pi(25/p_{k})])+\pi(5))/2$
$\pi(12.5)=5,\pi(8.33)=4,\pi(5)=3,\pi(3.57)=2,\pi(2.27)=1$
$\pi^{(2)}(25)=((5+4+3+2+1)+3)/2$ $\pi^{(2)}(25)=(15+3)/2=18/2=9$ (7)
$(5)=(6)=(7)$

Therefore, both formulas are correct.

## Identity

Theorem 2.

$\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})]-\sum_{k=\pi(\sqrt{N})+1}^{\pi(N/2)}[\pi(N/p_{k})]=(\pi(\sqrt{n}))^{2}$
(8)

Where, $\pi(x)$ is the prime counting function, $p_{k}$ is the $k^{th}$ prime.

Proof. As proved in theorem 1 and illustrated with Example 1, both of the
formulas give the correct count of the number of semi-primes not exceeding a
given number $N$, it is logical to equate the two formulas, i.e.

$(1)=(3)$
$\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})-k+1]=((\sum_{k=1}^{\pi((N/2))}[\pi(N/p_{k})])+\pi(\sqrt{N}))/2$

Simplifying LHS,

$=\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})]-\sum_{k=1}^{\pi(\sqrt{N})}[k-1]$
$=\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})]-(\pi(\sqrt{N})(\pi(\sqrt{N})+1))/2+\pi(\sqrt{N})$
$=\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})]-((\pi(\sqrt{N}))^{2}-\pi(\sqrt{N})+2\pi(\sqrt{N}))/2$
$=(2(\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})])-(\pi(\sqrt{N}))^{2}+\pi(\sqrt{N}))/2$

Equating LHS and RHS

$(2(\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})])-(\pi(\sqrt{N}))^{2}+\pi(\sqrt{N}))/2=$
$((\sum_{k=1}^{\pi((N/2))}[\pi(N/p_{k})])+\pi(\sqrt{N}))/2$

Multiplying both sides by 2,

$2(\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})])-(\pi(\sqrt{N}))^{2}+\pi(\sqrt{N})=$
$(\sum_{k=1}^{\pi((N/2))}[\pi(N/p_{k})])+\pi(\sqrt{N})$

Subtracting both sides by $(\pi(\sqrt{N})-(\pi(\sqrt{N}))^{2})$

$2(\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})])=(\sum_{k=1}^{\pi((N/2))}[\pi(N/p_{k})])+(\pi(\sqrt{N}))^{2}$

Subtracting both sides by $(\sum_{k=1}^{\pi((N/2))}[\pi(N/p_{k})])$

$2(\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})])-(\sum_{k=1}^{\pi((N/2))}[\pi(N/p_{k})])=(\pi(\sqrt{N}))^{2}$

which is equivalent to,

$2(\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})])-\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})]-\sum_{k=\pi(\sqrt{N})+1}^{\pi(N/2)}[\pi(N/p_{k})]$
$=(\pi(\sqrt{N}))^{2}$
$\implies\sum_{k=1}^{\pi(\sqrt{N})}[\pi(N/p_{k})]-\sum_{k=\pi(\sqrt{N})+1}^{\pi(N/2)}[\pi(N/p_{k})]=(\pi(\sqrt{N}))^{2}$

Example 2.
For $N=25$, RHS=

$(\pi(5))^{2}$ $=3^{2}=9$

LHS=

$\sum_{k=1}^{\pi(5)}[\pi(25/p_{k})]-\sum_{k=\pi(5)+1}^{\pi(12.5)}[\pi(25/p_{k})]$
$\pi(12.5)=5,\pi(8.33)=4,\pi(5)=3,\pi(3.57)=2,\pi(2.27)=1$ $=(5+4+3)-(2+1)$
$=12-3=9$

Therefore, LHS=RHS

## Conclusion

In this research, a new formula for the semi-prime counting function was
found. This has been shown to be a fairly straightforward method for
determining the number of semi-prime numbers that do not exceed a given
integer $N$. This yields an identity made up of many prime counting functions
when compared to the previously known formula.

## Forward looking

Similar, to how two formulas for semi-prime counting functions were found, it
might be possible to find more than one formula for k-almost prime numbers. A
k-almost prime is a composite number that is the product of exactly $k$ (few
of which might be equal) primes[2]. They then could be equated to find more
such identities, which might give insights into the prime counting function
itself.

## References

  * [1] Sloane, N. J. A. (ed.). ”Sequence A001358”. The On-Line Encyclopedia of Integer Sequences. OEIS Foundation.
  * [2] Conway, J. H.; Dietrich, H.; O’Brien, E. A. ”Counting Groups: Gnus, Moas, and Other Exotica.” Math. Intell. 30, 6-18, 2008.
  * [3] Ishmukhametov, Sh. T.; Sharifullina, F. F. (2014). ”On distribution of semiprime numbers”. Russian Mathematics. 58 (8): 43–48. doi:10.3103/S1066369X14080052. MR 3306238
  * [4] Weisstein, Eric W. ”Prime Counting Function”. MathWorld.
  * [5] Sloane, N. J. A. (ed.). ”Sequence A001358”. The On-Line Encyclopedia of Integer Sequences. OEIS Foundation.
  * [6] ”Sequence A006881”. The On-Line Encyclopedia of Integer Sequences. OEIS Foundation.


# Global monopoles in the two-Higgs-doublet-model

Richard A. Battye, Steven J. Cotterill and Dominic G. Viatic  Jodrell Bank
Centre for Astrophysics, Department of Physics and Astronomy, University of
Manchester, Oxford Road, Manchester M13 9PL

###### Abstract

We discuss monopoles formed due to the spontaneous breakdown of a global
$SO(3)_{\rm HF}$ symmetry within the global two-Higgs doublet model. We
explain that the Higgs sector dynamics can be described in terms of two
vectors one of which is null, $R^{A}=(R^{0},R^{a},R^{4},R^{5})$ for $a=1,2,3$,
with 5 independent components describing the Higgs family symmetry and
another, $n^{a}$, with 3 independent components related to the “would-be”
Goldstone bosons. When formed from random initial conditions we find that
monopoles are formed with a charged vacuum in the centre which couples the two
fields together. We find a spherical symmetric solution which is an
approximately uniform, unit winding of the sphere in both the $R^{a}$ and
$n^{a}$ vectors. These global monopoles are closely related to the Nambu
monopole. The additional complexity and structure contained in these monopoles
does not appear to prevent the scaling of their density.

##  1 Introduction

Monopoles are topological defects which arise from the breaking of a spherical
symmetry [1] and can be formed during the breaking of both global or
local/gauge symmetries. Gauge monopoles are commonly predicted in GUTs [2] as
they typically occur in phase transitions where the “little group” contains a
U(1) symmetry. As with domain walls, gauge monopoles can present issues for
late-time cosmology (see, for example, ref. [3]). The energy density due to
monopoles is expected to scale like matter in an FRW universe, and the ratio
of the two will remain constant. If the initial density of monopoles is larger
than that of matter, it will remain so throughout cosmic history preventing a
matter dominated epoch.

A commonly suggested solution to this issue is a period of cosmic inflation in
the early Universe which could dilute the initial energy density of monopoles
to a level that would be compatible with the fact that we have not observed
them (see, for example, [4]). However, the global monopoles considered in this
article are predicted to emerge at the electroweak scale and are, therefore,
expected to be post-inflation monopoles.

The Higgs mechanism for electroweak symmetry breaking was verified by the
measurement of a Higgs boson [5, 6] of mass 125 GeV at the LHC [7]. The
properties of this scalar particle so far match those predicted for the
Standard Model (SM) Higgs scalar [8, 9]. Nonetheless, current experimental
measurements do not prohibit the existence of more scalar particles. One
minimal and theoretically well-motivated extension which can be made to the SM
is to introduce a second complex Higgs doublet into the theory. This is the
so-called two-Higgs-doublet-model (2HDM) [10, 11, 12].

The 2HDM allows for the emergence of a variety of topological defects, such as
domain walls, vortices and global monopoles, from the breaking of accidental
symmetries which the model can possess under certain parameter choices [13,
14, 15, 16, 17, 18, 19, 20, 21, 22]. Phase transitions associated with the
spontaneous breaking of these symmetries can leave relic topological defects
which can serve as probes of high energy physics in the early Universe [23,
24, 25]. In this letter we will focus our attention on the
$\text{SO}(3)_{\text{HF}}$-symmetric 2HDM which predicts global monopoles
[13]. The key question we want to answer is whether the phenomena found in
ref. [17], whereby a massive photon (what is often termed a “charged vacuum”)
is predicted at the centre of the domain wall when the defects are formed from
random initial conditions, is a generic one for other topological defects in
the 2HDM.

We will ignore the SM gauge degrees of freedom in our simulations since these
make the simulations technically more difficult. The introduction of gauge
fields will introduce new length-scales into the problem which would need to
be carefully dealt with in the simulations and it is our strong belief, based
on experience of the evolution of topological defects in other contexts, that
is unlikely to have a significant impact on the qualitative picture of the
topological defects produced.

##  2 Two-Higgs-doublet-model with $\text{SO}(3)_{\text{HF}}$ symmetry

The Lagrangian density for the model is

$\mathcal{L}=\sum_{i=1}^{2}(\partial^{\mu}\Phi_{i})^{\dagger}(\partial_{\mu}\Phi_{i})-V(\Phi_{1},\Phi_{2})\,,$
(1)

where $\Phi_{1}$ and $\Phi_{2}$ are two complex doublet fields, and the
potential is given by

$\begin{split}V(\Phi_{1},\Phi_{2})=&-\mu_{1}^{2}\left(\Phi_{1}^{\dagger}\Phi_{1}+\Phi_{2}^{\dagger}\Phi_{2}\right)+\lambda_{1}\left(\Phi_{1}^{\dagger}\Phi_{1}+\Phi_{2}^{\dagger}\Phi_{2}\right)^{2}\\\
&+\lambda_{4}\left[(\Phi_{1}^{\dagger}\Phi_{2})(\Phi_{2}^{\dagger}\Phi_{1})-(\Phi_{1}^{\dagger}\Phi_{1})(\Phi_{2}^{\dagger}\Phi_{2})\right]\,.\end{split}$
(2)

This potential possesses a $\text{SO}(3)_{\text{HF}}\cong
SU(2)_{\text{HF}}/Z_{2}$ symmetry and has 3 real parameters;
$\mu_{1}^{2},\lambda_{1}$ and $\lambda_{4}$111The parameters $\mu_{1}$,
$\lambda_{1}$ and $\lambda_{4}$ have been chosen to fit with standard notation
eg. [13]., where $\lambda_{4}$ prevents the potential from being symmetric
under a larger symmetry group, but is also crucial for preventing neutral
vacuum violation [13]. Under a $\text{SO}(3)_{\text{HF}}$ transformation the
complex scalar Higgs doublets, $\Phi_{1}$ and $\Phi_{2}$, transform as
$\Phi_{1}\rightarrow
e^{-i\alpha}c_{\gamma}\Phi_{1}+e^{-i\beta}s_{\gamma}\Phi_{2}$,
$\Phi_{2}\rightarrow-e^{i\beta}s_{\gamma}\Phi_{1}+e^{i\alpha}c_{\gamma}\Phi_{2}$,
where $0\leq\alpha,\beta,\gamma<\pi$ with $s_{\gamma}=\sin\gamma$ and
$c_{\gamma}=\cos\gamma$, that is, an $SO(3)$ Higgs Family (HF) rotation
between the two Higgs doublets which is broken to $U(1)$ and hence the vacuum
manifold contains an extra $S^{2}$. This is an additional symmetry on top of
those associated with the SM, see ref. [13] for a detailed discussion.

The most general parametrization of the 2HDM vacuum is [12]

${\bar{\Phi}}_{1}=\frac{1}{\sqrt{2}}\left(\begin{matrix}0\\\
v_{1}\end{matrix}\right),\quad{\bar{\Phi}}_{2}=\frac{1}{\sqrt{2}}\left(\begin{matrix}v_{+}\\\
v_{2}e^{i\xi}\end{matrix}\right)\,.$ (3)

In the following section we will use $v_{1},v_{2},v_{+}$ and $\xi$ to
parameterise the fields. However, for the moment let us think of them as just
the vacuum and expand around a neutral vacuum with $v_{+}=0$ and $\xi=0$. The
2HDM has 5 physical scalar particles: 2 neutral CP-even states, $h$ and $H$,
one CP-odd neutral state, $A$, and 2 charged states, $H^{\pm}$. The other
three scalar degrees of freedom correspond to would-be Goldstone bosons,
$G^{0}$ and $G^{\pm}$, which are absorbed into the longitudinal components of
the electroweak gauge bosons, $W^{\pm}$ and $Z^{0}$. Expressions for the
masses of the scalar Higgs particles, $h$, $H$, $A$ and $H^{\pm}$, are
obtained as eigenvalues of the Hessian matrix of (2) using the parametrization

$\Phi_{1}=\left(\begin{matrix}\varphi_{1}^{+}\\\
\frac{1}{\sqrt{2}}v_{1}+\varphi_{1}+ia_{1}\end{matrix}\right),\quad\Phi_{2}=e^{i\xi}\left(\begin{matrix}\varphi_{2}^{+}\\\
\frac{1}{\sqrt{2}}v_{2}+\varphi_{2}+ia_{2}\end{matrix}\right)\,,$ (4)

where $\varphi_{i}^{+}$ are complex scalar fields. One finds that, due to the
greatly reduced parameter space required to obtain $SO(3)_{\text{HF}}$
symmetry, one of the CP-even scalars along with the CP-odd scalar are both
massless, $M_{H}=0,M_{A}=0$ and the other CP even scalar is the SM Higgs with
$M_{h}^{2}=2\lambda_{1}v_{\text{SM}}^{2}$. The charged scalar mass is given by
$M_{H^{\pm}}^{2}=-\frac{1}{2}\lambda_{4}v_{\text{SM}}^{2}.$ Therefore, to
ensure a stable vacuum, the quartic coupling parameters must satisfy the
inequalities, $\lambda_{1}>0,\lambda_{4}<0$, such that the non-zero scalar
masses are real and positive. The fact that there are massless Higgs states
would probably make this specific manifestation of the model not
phenomenologically viable, but many of the features found here are likely to
persist in more realistic models.

##  3 Field parameterization

Figure 1: Two-dimensional slices showing spatial distributions of
$\sqrt{n^{a}n^{a}}$ (left), $R^{\mu}R_{\mu}$ (middle) and $R^{a}R^{a}$ (right)
from a (3+1) dimensional simulation of the
$\text{SO}(3)_{\text{HF}}$-symmetric 2HDM with $\epsilon=0.35$. Note that the
colour bars have been chosen so that all three quantities look very similar.
The simulation with $P=512$, $\Delta t=0.2$ and $\Delta x=0.9$ was run for a
time $t=240$ and dissipation was removed at $t=30$. Time increases downwards
with the top plot being $t=30$, next two $t=60,120$ and the bottom plot
corresponds to the end of the simulation, $t=240$ The global monopoles are
located at the point where $R^{a}R^{a}=0$ which correlates strongly with
maxima in $R^{\mu}R_{\mu}$ indicating that charge neutrality condition is
violated at the core of the monopoles. Away from the monopoles
$\sqrt{n^{a}n^{a}}=1/2$ and this quantity takes its minimum value at the
monopole cores. As time evolves the density of monopoles can been seen to
visually reduce roughly scaling like $t^{-3}$. See later for further
discussion of the scaling exponent.

It turns out, that in terms of understanding the topological defect solutions,
the doublet fields $\Phi_{1}$ and $\Phi_{2}$ are not the best way to see the
structure of the vacuum. This is because the global symmetries which are being
spontaneously broken are due to the internal symmetries between the components
of the doublets. It was explained in ref.[13] that the additional degrees of
freedom found in the 2HDM can be parameterized in terms of a vector $R^{A}$
with $A=0,1,..,5$ which is null in the sense that $R^{A}R_{A}=0$ and hence
there are 5 degrees of freedom. These can be thought of as corresponding to
the 5 Higgs particles in the 2HDM. For $A=\mu=0,1,.,3$ we have that [26, 27,
28, 13]

$R^{\mu}=\Phi^{\dagger}\left(\sigma^{\mu}\otimes
I_{2}\right)\Phi=\left(\begin{matrix}\Phi_{1}^{\dagger}\Phi_{1}+\Phi_{2}^{\dagger}\Phi_{2}\\\
\Phi_{1}^{\dagger}\Phi_{2}+\Phi_{2}^{\dagger}\Phi_{1}\\\
-i[\Phi_{1}^{\dagger}\Phi_{2}-\Phi_{2}^{\dagger}\Phi_{1}]\\\
\Phi_{1}^{\dagger}\Phi_{1}-\Phi_{2}^{\dagger}\Phi_{2}\end{matrix}\right)\,,$
(5)

where $\sigma^{\mu}$ are the Pauli matrices including the identity,
$\sigma^{0}=I_{2}$ and we have introduced the multiplet,

$\Phi\equiv\left(\begin{matrix}\Phi_{1}\cr\Phi_{2}\end{matrix}\right)\,.$ (6)

The other two components are

$R^{4}=\Phi_{1}^{T}i\sigma^{2}\Phi_{2}-\Phi_{2}^{\dagger}i\sigma^{2}\Phi_{1}^{*}\,,\quad
R^{5}=-i\left(\Phi_{1}^{T}i\sigma^{2}\Phi_{2}+\Phi_{2}^{{\dagger}}i\sigma^{2}\Phi_{1}^{*}\right)\,,$
(7)

which can be brought together in a complex field
$\tilde{R}=R^{4}+iR^{5}=2\Phi^{T}_{1}i\sigma^{2}\Phi_{2}$. In a neutral
vacuum, where the photon mass is zero, we have that $v_{+}=0$ and hence we
find that $R^{\mu}R_{\mu}=0$ and $\tilde{R}=0$. It has already been
established [17] that simulations of 2HDM domain walls predict a violation of
this neutral vacuum condition in the core of the defect as a general feature
emerging from random initial conditions and we also expect this to be the case
for the global monopoles considered here.

Now define the 3-component field,

$n^{a}=-\Phi^{\dagger}(I_{2}\otimes\sigma^{a})\Phi=-\sum_{n=1}^{2}\Phi_{n}^{\dagger}\sigma^{a}\Phi_{n}\,,$
(8)

for $a=1,2,3$. We will see that this encodes the 3 degrees of freedom which
will be “eaten” by the gauge fields as part of the Higgs mechanism.

We can parameterise the field as $\Phi=e^{i\chi/2}(I_{2}\otimes
U_{L}){\bar{\Phi}}$ where $0\leq\chi<2\pi$, $U_{L}$ is an element of
$SU(2)_{L}$ and
${\bar{\Phi}}={1\over\sqrt{2}}\left(0,v_{1},v_{+},v_{2}e^{i\xi}\right)$, that
is, a general element of the vacuum manifold. There are a total of eight
degrees of freedom $v_{1},v_{2},v_{+},\xi,\chi$ and the three degrees of
freedom contained within $U_{L}$. Immediately we see that
$R^{\mu}={\bar{\Phi}}^{{\dagger}}(\sigma^{\mu}\otimes
I_{2}){\bar{\Phi}}={\bar{R}}^{\mu}$ with

${\bar{R}}^{\mu}={1\over
2}\left(\begin{matrix}v_{1}^{2}+v_{2}^{2}+v_{+}^{2}\\\ 2v_{1}v_{2}\cos\xi\\\
2v_{1}v_{2}\sin\xi\\\ v_{1}^{2}-v_{2}^{2}-v_{+}^{2}\end{matrix}\right)\,,$ (9)

which only depends on $v_{1},v_{2},v_{+}$ and $\xi$, and is independent of
$\chi$ and $U_{L}$. We also see that
$\tilde{R}=2e^{i\chi}{\bar{\Phi}}_{1}^{T}i\sigma^{2}{\bar{\Phi}}_{2}=e^{i\chi}\tilde{\bar{R}}$
where $\tilde{\bar{R}}=-v_{+}v_{1}$. This is also independent of $U_{L}$ as a
consequence of $U_{L}^{T}i\sigma^{2}U_{L}=i\sigma^{2}$ and this illustrates
that the phase of the complex scalar field, ${\tilde{R}}$, isolates $\chi$.
Therefore, $R^{\mu}$ and $\tilde{R}$ encode the four new degrees of freedom
added in going from the SM to the 2HDM as well as the SM Higgs.

Under the same parameterization we see that

$n^{a}=-{\bar{\Phi}}^{\dagger}\left(I_{2}\otimes(U_{L}^{\dagger}\sigma^{a}U_{L})\right){\bar{\Phi}}\,,$
(10)

which is independent of $\chi$. If we write
$U_{L}=\exp\left[iG^{a}\sigma^{a}/(2v_{\rm SM}\right)]$ where the $G^{a}$ are
the three “would be” Goldstone bosons, then
$U_{L}^{\dagger}\sigma^{a}U_{L}={\cal R}^{ab}\sigma^{b}$ where

${\cal
R}^{ab}={\hat{G}}^{a}{\hat{G}}^{b}+(\delta^{ab}-{\hat{G}}^{a}{\hat{G}}^{b})\cos\gamma+\epsilon^{abc}{\hat{G}}^{c}\sin\gamma\,,$
(11)

$\gamma=|G|/v_{\rm SM}$ and ${\hat{G}}^{a}=G^{a}/|G|$, that is, a rotation by
an angle $\gamma$ around the unit vector ${\hat{G}}^{a}$. Hence, we can write
$n^{a}={\cal R}^{ab}{\bar{n}}^{b}$ where
${\bar{n}}^{a}=-{\bar{\Phi}}^{\dagger}(I_{2}\otimes\sigma^{a}){\bar{\Phi}}$
and in a neutral vacuum state we have $\bar{n}^{a}=\textstyle{1\over 2}v_{\rm
SM}^{2}(0,0,1)$. This shows that action of the SM degrees of freedom leads to
a rotation of ${\bar{n}}^{a}$, but remember that they have no impact on
$R^{\mu}$.

##  4 Formation and evolution of global monopoles

(a)

(b)

(c)

(d)

Figure 2: (top) Isosurfaces of $R^{\mu}R_{\mu}=0.1$ \- a proxy for the
position of the monopole - as a function of time for $\epsilon=0.35$. The
times used are $t=30$ (top-left), $t=60$ (top-right), $t=120$ (middle-left)
and $t=240$ (middle-right) as in Fig. 1. In the $t=240$ plot there is a small
blue box which isolates a small region used in Fig. 3. At the bottom we show
the evolution of the number of monopoles, for various values of $\epsilon$,
averaged over 10 realisations as described in the text. Error bars represent
the numerical scatter via the standard deviation.

We have performed (3+1) dimensional simulations for the global scalar field
theory of the 2HDM with $\text{SO}(3)_{\text{HF}}$ symmetry by evolving the
equations of motion on a regular grid of $P^{3}$ points for $P=512$ and
$P=1024$ with Minkowski metric and periodic boundary conditions (for details
of the simulation procedure, see [17, 18]). Temporal derivatives are
approximated to second order and spatial derivatives to fourth order. The
initial conditions are created by placing small random numbers in all of the
fields and then evolving them under dissipative dynamics for around 200 time
steps, before allowing the fields to evolve under the full dynamics up to the
light crossing time set by $P\Delta x/2$ where $\Delta x=0.9$ when $P=512$ and
$\Delta x=0.45$ when $P=1024$ is the spatial stepsize. We use a timestep of
$\Delta t=0.2$ for $P=512$ and $0.1$ for $P=1024$ in order to maintain
numerical stability.

The units can be rescaled so the energy density is measured in units of
$\mu_{1}^{4}/\lambda_{1}$ and length/time in units of $\mu_{1}^{-1}$. This
leaves a single dimensionless parameter $\epsilon=M_{H^{\pm}}/M_{h}$. We ran
simulations for three different choices of $\epsilon=0.35,0.5$ and $0.7$, with
the largest corresponding to $M_{H\pm}\approx 88\,{\rm GeV}$. This is quite
low for the charged Higgs mass, but having a low value of $\epsilon$ makes the
simulation more numerically tractable. We expect the qualitative features of
the simulations to be similar over a wide range of parameters, and in
particular for much larger values of $\epsilon$, but this will be confirmed in
future work.

In Fig. 1 we present two dimensional slices of the spatial distribution of the
scalar quantities $\sqrt{n^{a}n^{a}}$, $R^{\mu}R_{\mu}$ and $R^{a}R^{a}$ for
the simulation with $\epsilon=0.35$. The centres of monopoles are defined to
be when $R^{a}R^{a}=0$. The right hand column indicates that monopoles form
during the period of dissipation and that their density is initially high but
that it rapidly decreases with time. The quantity $R^{\mu}R_{\mu}$ is non-zero
at the cores of the monopoles and the maxima of this quantity are strongly
correlated with the positions of the monopoles. This indicates that inside the
monopoles the vacuum is charged with $v_{+}\neq 0$, which would correspond to
a non-zero photon mass in the centre of the monopole. We also find that
$\sqrt{n^{a}n^{a}}$ has its minimum value at the same place as the position of
the monopoles.

The 3D dynamics of the monopoles for $\epsilon=0.35$ are illustrated in Fig. 2
where we have plotted isosurfaces of the quantity $R^{\mu}R_{\mu}$ which,
based on the previous paragraph, we now take as a proxy for the positions of
the monopoles. Complex dynamics takes place, presumably with the monopoles
interacting with each other and annihilating with anti-monopoles. Close
examination of the simulation output appears to reveal a number of dumb-bell-
shaped configurations which could be monopoles/anti-monopoles interacting via
a bridge comprising a string-like structure, motivating detailed investigation
of the simulations that is beyond the scope of the present work. Similar
objects have been examined in the context of the Standard Model [29, 30] and
$SU(5)$ Grand Unified Theories [31].

We also present the number of monopoles, $N_{\text{mono}}$, as a function of
time obtained as an average over 10 realizations and we give the best fit
scaling exponents in Table 1. We note that in the case of $\epsilon=0.7$ it
was necessary to use a lower value of $\Delta x$ (keeping $P\Delta x$
constant) in order to resolve the length scale $\propto\epsilon^{-1}$ \-
simulations with the higher value of $\Delta x$ indicated some deviation from
th expected scaling law which was restored by higher resolution. We find that
monopoles in the $\text{SO}(3)_{\text{HF}}$-symmetric 2HDM scale as $N_{\rm
mono}\propto t^{-3}$ which is compatible with what is seen to take place in
the global $O(3)$ model [32, 33]. It appears, therefore, that whatever leads
to the violation of the neutral vacuum inside the monopoles does not lead to a
significant change in the expected evolution, at least within the dynamical
range of these simulations. It would be interesting to perform simulations
with a larger dynamic range in order to test this.

$\epsilon$ | $0.35$ | $0.5$ | $0.7$
---|---|---|---
Scaling Exponent | $-2.9\pm 0.1$ | $-2.93\pm 0.08$ | $-2.98\pm 0.08$
Table 1: Measured scaling exponents for three different values of $\epsilon$.
Those for $\epsilon=0.35$ and $0.5$ used $P=512$ and $\Delta x=0.9$, while
that for $\epsilon=0.7$ used $P=1024$ and $\Delta x=0.45$. In all three cases
we see that they are compatible with the number of monopoles scaling $\propto
t^{-3}$.

(a) $R^{a}$

(b) $R^{a}$ after a global rotation

(c) $n^{a}$

(d) $n^{a}$ after a global rotation

Figure 3:  We present $R^{a}$ (top) and $n^{a}$ (bottom) around one of the
monopoles at the end of the simulation. The left side shows the vectors that
are directly calculation from the simulation and on the right we show the same
vector fields after a global rotation has been applied. Note that it is a
different rotation for $R^{a}$ and $n^{a}$. We have applied a similar
procedure, but with different rotations, to some of the other monopoles in the
simulation and are able to get similar results.

To further investigate the nature of the field configurations that are formed
in the simulations, we have isolated a region at late times which appears to
contain a single monopole. In Fig. 3 we have attempted within constraints of
the resolution of the simulation, which only has a few grid squares for each
monopole, to plot the vectors $R^{a}$ and $n^{a}$. At first glance the
configurations indicate some spatial gradients, but are not immediately
identifiable as monopoles. However, we have found that we can apply global
rotations (not the same for both) to orient all of the vectors so that they
point in a direction that is approximately radial in the vicinity of the
monopole.

##  5 Global monopole solution

A global monopole solution was proposed [13] which can written as

$\Phi={1\over\sqrt{2}}v_{\rm SM}f(r)\left({\hat{\bf r}}\cdot\sigma\otimes
I_{2}\right)\left(\begin{matrix}0\cr 1\cr 0\cr 0\end{matrix}\right)\,,$ (12)

where $f(r)$ is a function just depending on the radial coordinate $r$ which
satisfies a second order differential equation that can be deduced from ref.
[13] with boundary conditions $f(0)=0$ and $f(\infty)=1$ and $v_{\rm
SM}=\mu_{1}/\sqrt{\lambda_{1}}$. It was constructed by setting $v_{1}=v_{\rm
SM}f(r)\cos\theta$, $v_{2}=v_{\rm SM}f(r)\sin\theta$, $v_{+}=0$ and
$\xi=\phi$.

This neutral vacuum solution has $R^{a}=\textstyle{1\over 2}v_{\rm
SM}^{2}[f(r)]^{2}{\cal Q}^{ab}{\hat{z}}^{b}$ where ${\cal
Q}^{ab}=2{\hat{r}}^{a}{\hat{r}}^{b}-\delta^{ab}$ and $n^{a}=\textstyle{1\over
2}v^{2}_{\rm SM}[f(r)]^{2}{\hat{z}}^{a}$, both of which are $\textstyle{1\over
2}v^{2}_{\rm SM}[f(r)]^{2}$ times a unit vector with $R^{a}={\cal
Q}^{ab}n^{b}$. We now believe that this is in fact not a monopole since
$R^{a}$ has zero winding number, and in fact is a sphaleron.

In light of $R^{a}$ and $n^{a}$ that we have discovered in the monopoles
formed from random initial conditions, we initially condsidered a monopole
ansatz that has $R^{a}=n^{a}\propto\hat{r}^{a}$

$\Phi=\frac{v_{\rm SM}f(r)}{2\sqrt{2}}\begin{pmatrix}-\sin\theta e^{-i\phi}\\\
\cos\theta+1\\\ \cos\theta-1\\\ \sin\theta e^{i\phi}\end{pmatrix}\,,$ (13)

which has $v_{1}=v_{\rm SM}f(r)\cos\frac{1}{2}\theta$, $v_{2}=v_{\rm
sM}f(r)\sin\frac{1}{2}\theta$, $v_{+}=0$ and $\xi=\phi$ with the
$\text{SU(2)}_{\rm L}$ rotation

$U_{L}=\begin{pmatrix}\cos\frac{1}{2}\theta&-\sin\frac{1}{2}\theta
e^{-i\phi}\\\ \sin\frac{1}{2}\theta
e^{i\phi}&\cos\frac{1}{2}\theta\end{pmatrix}\,.$ (14)

The expressions for $v_{1}$ and $v_{2}$ are half angle versions of (12) and
$U_{L}$ generates a Nambu monopole in the SM [34], but in this case there is
no need for a string to be attached to one of the poles. This is because the
divergent part of the gradient energy associated with the winding of $n^{a}$
cancels due to the extra winding in $R^{a}$. However, this ansatz cannot
describe the monopoles that we see in the simulations as $R_{\mu}R^{\mu}=0$
(because $v_{+}=0$) and $\Phi$ will be forced to zero at the core of the
monopole by the gradient energy terms.

A simple extension that allows for charge-breaking in the core of the monopole
is

$\Phi=\frac{v_{\rm SM}}{2\sqrt{2}}\begin{pmatrix}-(f-f_{+})\sin\theta
e^{-i\phi}\\\ (f-f_{+})\cos\theta+(f+f_{+})\\\
(f-f_{+})\cos\theta-(f+f_{+})\\\ (f-f_{+})\sin\theta
e^{i\phi}\end{pmatrix}\,,$ (15)

where $f=f(r)$ and $f_{+}=f_{+}(r)$ for which $R^{a}=n^{a}=\frac{v_{\rm
SM}^{2}}{2}(f^{2}-f_{+}^{2})\hat{r}^{a}$. This ansatz can be written as
$(I_{2}\otimes U_{L})\bar{\Phi}$ (as any configuration in the 2HDM can) but it
is a much simpler expression in this form. Note that fixing $f_{+}=0$ returns
the ansatz of equation (13) and the gradient energy does not force $\Phi=0$ at
the centre of the monopole - only $f=f_{+}$. Under this ansatz,
$R_{\mu}R^{\mu}=v_{\rm SM}^{4}f^{2}f_{+}^{2}$ and therefore there will be a
massive photon and neutral vacuum violation whenever both $f$ and $f_{+}$ are
simultaneously non-zero.

One can calculate the energy to be

$\hat{E}=\int
r^{2}dr\bigg{\\{}\frac{1}{2}\bigg{(}\frac{df}{dr}\bigg{)}^{2}+\frac{1}{2}\bigg{(}\frac{df_{+}}{dr}\bigg{)}^{2}+\frac{(f-f_{+})^{2}}{2r^{2}}-\frac{1}{2}(f^{2}+f_{+}^{2})+\frac{1}{4}(f^{2}+f_{+}^{2})^{2}+\epsilon^{2}f^{2}f_{+}^{2}\bigg{\\}}\,,$
(16)

where $E=4\pi v_{\rm SM}^{2}\hat{E}/\mu_{1}$, $v_{\rm
SM}=\mu_{1}/\sqrt{\lambda_{1}}$ and we have rescaled the length scale, $r\to
r/\mu_{1}$, so that the solution only depends upon the single parameter,
$\epsilon=M_{\rm H\pm}/M_{h}$. Such a solution must satisfy,

$\displaystyle\frac{1}{r^{2}}\frac{d}{dr}\bigg{(}r^{2}\frac{df}{dr}\bigg{)}-\frac{f-f_{+}}{r^{2}}-\bigg{[}f^{2}+f_{+}^{2}-1+2\epsilon^{2}f_{+}^{2}\bigg{]}f$
$\displaystyle=$ $\displaystyle 0\,,$ (17)
$\displaystyle\frac{1}{r^{2}}\frac{d}{dr}\bigg{(}r^{2}\frac{df_{+}}{dr}\bigg{)}+\frac{f-f_{+}}{r^{2}}-\bigg{[}f^{2}+f_{+}^{2}-1+2\epsilon^{2}f^{2}\bigg{]}f_{+}$
$\displaystyle=$ $\displaystyle 0\,,$ (18)

with the boundary conditions $f(0)=f_{+}(0)$,
$f^{\prime}(0)=-f_{+}^{\prime}(0)$, $f(\infty)=1$ and $f_{+}(\infty)=0$.

(a)

(b)

(c)

(d)

Figure 4: We present the numerical solutions of charge-breaking monopoles, for
a range of $\epsilon$, with $f(r)$ displayed in the top left plot, $f_{+}(r)$
in the top right plot, $\hat{R}_{\mu}\hat{R}^{\mu}$ on the bottom left and the
energy density on the bottom right. Note that
$\hat{R}_{\mu}=\lambda_{1}R_{\mu}/\mu_{1}^{2}$ and that the y-axis for the
$f(r)$ plot does not start at zero ($f(0)=f_{+}(0)$).

In Figure 4 we present the solutions to the above equations (found numerically
using a grid of $10000$ points and $\Delta x=0.01$) as well as
$\hat{R}_{\mu}\hat{R}^{\mu}$ \- where $R_{\mu}=v_{\rm SM}^{2}\hat{R}_{\mu}$ \-
and the energy density. It is clear that these solutions have
$R_{\mu}R^{\mu}\neq 0$ close to the monopole core and that its value decreases
with increasing $\epsilon$ \- a feature that is also present in the 3D
simulations. We can also see that the length scale associated with $\epsilon$
corresponds to the length scale of $f_{+}$, while $f$ changes only marginally.
Finally, although the energy of these solutions is infinite because they are
global monopoles, the energy density decreases with increasing $\epsilon$,
which likely means that the total energy of these monopoles in a more
realistic gauged theory will increase with $\epsilon$.

##  6 Discussion

The 2HDM can predict a variety of topological defects when accidental
symmetries are broken. Specifically, for
$\text{SU}(2)_{L}\times\text{U}(1)_{Y}$ preserving symmetries, there are three
domain wall solutions, two vortex solutions and one global monopole solution
[13]. Here, we have focused exclusively on the 2HDM with an
$\text{SO}(3)_{\text{HF}}$ symmetry that produces global monopoles.

Our simulations clearly show the formation of monopoles and that they evolve
as one might expect with $N_{\rm mono}\propto t^{-3}$ for the values of
$\epsilon$ that are easily investigated with our simulations and within the
dynamic range. However, the structure of the monopoles formed is very
different to the solution expected from ref. [13]. There appears to be a
violation of the neutral vacuum in the centre of the monopole with
$R^{\mu}R_{\mu}$ and $v_{+}$, and the field configuration has non trivial
structures in the vector field $n^{a}$ which we have shown encode the standard
model degrees of freedom. A similar phenomenon was found in the case of domain
walls [17, 22], albeit in a lower dimensional situation.

We have presented a new monopole ansatz for the 2HDM that is closely related
to the Nambu monopole, but it is stable and has no requirement for there to be
a string emerging from one of the poles that connects it to an anti-monopole.
It is spherically symmetric and has the property that
$R^{a}=n^{a}\propto\hat{r}^{a}$ and we have shown that this is also the case
for the monopoles formed in our simulations, after a constant rotation is
performed. The energy minimising solution under this ansatz always has non-
zero $R_{\mu}R^{\mu}$ (as $f_{+}=0$ is not a solution) in the core of the
monopole, but it decreases as $\epsilon$ is increased, presumably going to
zero in the limit where $\epsilon\to\infty$.

We note as a final point, that topological defects which violate the neutral
vacuum conditions in the core and hence become superconducting can have some
novel interaction properties with photons and other relativistic particles
[35]. Here, we have confirmed that this indeed happens in the 2HDM model with
$\text{SO(3)}_{\rm HF}$ symmetry and it is clearly worth exploring the
cosmological consequences of this. There are a number of other possibilities,
for example, mixing of photons with Z and W bosons, generation of primordial
magnetic fields and possibly even the violation of SM constraints on the
interaction of leptons and quarks, although this will depend on the structure
of the Yukawa sector of the theory which we have not discussed here. Moreover,
Nambu monopoles in the SM have been discussed in the context of primordial
magnetic fields [36, 37]. We have shown that the production of such monopoles
maybe much more natural in the context of the 2HDM and this deserves further
investigation.

## Acknowledgements

RB and DV would like to thank Apostolos Pilaftsis for his collaboration on
related work on domain walls.

## References

  * [1] A. Vilenkin and E.P.S Shellard.  Cosmic strings and other topological defects.  Cambridge monographs on mathematical physics. Cambridge University Press, Cambridge, 1994.
  * [2] Gerard ’t Hooft.  Magnetic Monopoles in Unified Gauge Theories.  Nucl. Phys., B79:276–284, 1974.  [,291(1974)].
  * [3] C. J. A. P. Martins and A. Achúcarro.  Evolution of local and global monopole networks.  Phys. Rev. D, 78:083541, Oct 2008.
  * [4] Sergey Burdin, Malcolm Fairbairn, Philippe Mermod, David Milstead, James Pinfold, Terry Sloan, and Wendy Taylor.  Non-collider searches for stable massive particles.  Phys. Rept., 582:1–52, 2015.
  * [5] Georges Aad et al.  Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC.  Phys. Lett., B716:1–29, 2012.
  * [6] Serguei Chatrchyan et al.  Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC.  Phys. Lett., B716:30–61, 2012.
  * [7] Georges Aad et al.  Combined Measurement of the Higgs Boson Mass in $pp$ Collisions at $\sqrt{s}=7$ and 8 TeV with the ATLAS and CMS Experiments.  Phys. Rev. Lett., 114:191803, 2015.
  * [8] Abdelhak Djouadi.  The Anatomy of electro-weak symmetry breaking. I: The Higgs boson in the standard model.  Phys. Rept., 457:1–216, 2008.
  * [9] Georges Aad et al.  Measurements of the Higgs boson production and decay rates and constraints on its couplings from a combined ATLAS and CMS analysis of the LHC pp collision data at $\sqrt{s}=7$ and 8 TeV.  JHEP, 08:045, 2016.
  * [10] T. D. Lee.  A Theory of Spontaneous T Violation.  Phys. Rev., D8:1226–1239, 1973.  [,516(1973)].
  * [11] Apostolos Pilaftsis and Carlos E. M. Wagner.  Higgs bosons in the minimal supersymmetric standard model with explicit CP violation.  Nucl. Phys., B553:3–42, 1999.
  * [12] G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, Marc Sher, and Joao P. Silva.  Theory and phenomenology of two-Higgs-doublet models.  Phys. Rept., 516:1–102, 2012.
  * [13] Richard A. Battye, Gary D. Brawn, and Apostolos Pilaftsis.  Vacuum Topology of the Two Higgs Doublet Model.  JHEP, 08:020, 2011.
  * [14] Minoru Eto, Masafumi Kurachi, and Muneto Nitta.  Constraints on two Higgs doublet models from domain walls.  Phys. Lett. B, 785:447–453, 2018.
  * [15] Minoru Eto, Masafumi Kurachi, and Muneto Nitta.  Non-Abelian strings and domain walls in two Higgs doublet models.  JHEP, 08:195, 2018.
  * [16] Ning Chen, Tong Li, Zhaolong Teng, and Yongcheng Wu.  Collapsing domain walls in the two-Higgs-doublet model and deep insights from the EDM.  JHEP, 10:081, 2020.
  * [17] Richard A. Battye, Apostolos Pilaftsis, and Dominic G. Viatic.  Simulations of domain walls in Two Higgs Doublet Models.  JHEP, 01:105, 2021.
  * [18] Richard A. Battye, Apostolos Pilaftsis, and Dominic G. Viatic.  Domain wall constraints on two-higgs-doublet models with ${Z}_{2}$ symmetry.  Phys. Rev. D, 102:123536, Dec 2020.
  * [19] Minoru Eto, Yu Hamada, and Muneto Nitta.  Topological structure of a Nambu monopole in two-Higgs-doublet models: Fiber bundle, Dirac’s quantization, and a dyon.  Phys. Rev. D, 102(10):105018, 2020.
  * [20] Minoru Eto, Yu Hamada, Masafumi Kurachi, and Muneto Nitta.  Dynamics of Nambu monopole in two Higgs doublet models. Cosmological Monopole Collider.  JHEP, 07:004, 2020.
  * [21] Minoru Eto, Yu Hamada, and Muneto Nitta.  Stable $Z$-strings with topological polarization in two Higgs doublet model.  11 2021.
  * [22] Kai Hong Law and Apostolos Pilaftsis.  On the Charged and CP-violating Kink Solutions in the Two Higgs Doublet Model.  10 2021.
  * [23] Theodore Garagounis and Mark Hindmarsh.  Scaling in numerical simulations of domain walls.  Phys. Rev., D68:103506, 2003.
  * [24] T.W.B. Kibble, George Lazarides, and Q. Shafi.  Walls Bounded by Strings.  Phys. Rev. D, 26:435, 1982.
  * [25] Kazunori Nakayama, Fuminobu Takahashi, and Norimi Yokozaki.  Gravitational waves from domain walls and their implications.  Phys. Lett., B770:500–506, 2017.
  * [26] M. Maniatis, A. von Manteuffel, and O. Nachtmann.  CP violation in the general two-Higgs-doublet model: A Geometric view.  Eur. Phys. J., C57:719–738, 2008.
  * [27] Igor P. Ivanov.  Minkowski space structure of the Higgs potential in 2HDM. II. Minima, symmetries, and topology.  Phys. Rev., D77:015017, 2008.
  * [28] C.C. Nishi.  CP violation conditions in N-Higgs-doublet potentials.  Phys. Rev. D, 74:036003, 2006.  [Erratum: Phys.Rev.D 76, 119901 (2007)].
  * [29] George Lazarides, Qaisar Shafi, and Tanmay Vachaspati.  Dirac plus Nambu monopoles in the Standard Model.  Phys. Rev. D, 104(3):035020, 2021.
  * [30] Teerthal Patel and Tanmay Vachaspati.  Structure of electroweak dumbbells.  2 2023.
  * [31] G. Lazarides and Q. Shafi.  Electroweak monopoles and magnetic dumbbells in grand unified theories.  Phys. Rev. D, 103:095021, May 2021.
  * [32] Masahide Yamaguchi.  Cosmological evolution of global monopoles.  Phys. Rev. D, 64:081301, 2001.
  * [33] D. P. Bennett and S. H. Rhie.  Cosmological evolution of global monopoles and the origin of large scale structure.  Phys. Rev. Lett., 65:1709–1712, 1990.
  * [34] Y. Nambu.  String-like configurations in the weinberg-salam theory.  Nuclear Physics B, 130(3):505–515, 1977.
  * [35] Richard A. Battye and Dominic G. Viatic.  Photon interactions with superconducting topological defects.  Phys. Lett. B, 823:136730, 2021.
  * [36] Tanmay Vachaspati.  Progress on cosmological magnetic fields.  Rept. Prog. Phys., 84(7):074901, 2021.
  * [37] Teerthal Patel and Tanmay Vachaspati.  Kibble mechanism for electroweak magnetic monopoles and magnetic fields.  Journal of High Energy Physics, 2022(1):59, Jan 2022.


# Control- & Task-Aware Optimal Design of Actuation System
for Legged Robots using Binary Integer Linear Programming

Youngwoo Sim, Guillermo Colin and Joao Ramos This work is supported by the
National Science Foundation via grant CMMI-2043339.The authors are with the
Department of Mechanical Science and Engineering at the University of Illinois
at Urbana - Champaign<EMAIL_ADDRESS>
###### Abstract

Athletic robots demand a whole-body actuation system design that utilizes
motors up to the boundaries of their performance. However, creating such
robots poses challenges of integrating design principles and reasoning of
practical design choices. This paper presents a design framework that guides
designers to find optimal design choices to create an actuation system that
can rapidly generate torques and velocities required to achieve a given set of
tasks, by minimizing inertia and leveraging cooperation between actuators. The
framework serves as an interactive tool for designers who are in charge of
providing design rules and candidate components such as motors, reduction
mechanism, and coupling mechanisms between actuators and joints. A binary
integer linear optimization explores design combinations to find optimal
components that can achieve a set of tasks. The framework is demonstrated with
200 optimal design studies of a biped with 5-degree-of-freedom (DoF) legs,
focusing on the effect of achieving multiple tasks (walking, lifting),
constraining the mass budget of all motors in the system and the use of
coupling mechanisms. The result provides a comprehensive view of how design
choices and rules affect reflected inertia, copper loss of motors, and force
capability of optimal actuation systems.

##  I Introduction

Athletic robots such as the MIT cheetah and Cassie have emerged to show
superior performance to their semi-static counterparts. These machines possess
the ability to execute agile movements, such as utilizing momentum to jump or
push heavy objects. Achieving such athletic motions requires a comprehensive
approach to whole-body control and systematic actuation system design [1, 2,
3], wherein multiple actuators are built to operate cooperatively across the
robot’s body. While the literature has made significant strides in whole-body
control of humanoids, the design problem pertaining to high-DoF systems still
demands further investigation.

The ability of athletic robots to generate large forces and acceleration is
harnessed by incorporating design principles at two distinct layers of the
actuation system. At the actuator level, the quasi-direct drive (QDD) paradigm
[4] suggests minimizing the inertia of the limbs. Using smaller gear ratios
minimizes the reflected inertia of individual actuators and reduces inertial
torques required to swing the limbs. In cases where gear ratios are not large
enough to amplify torque sufficiently, this can be addressed at the system
level by employing cooperative (coupled or parallel) actuation [5], where
actuators are coupled together to combine their torque output and meet the
torque requirements of specific joints.

However, as the number of actuators of a system grows, designers face multi-
faceted issues in integrating the aforementioned principles. First, the design
considerations for multi-actuator systems, such as humanoid legs, are not
trivial. Design principles often provide general intuition but leave specific
design choices to the discretion of designers. Consequently, designers need to
analyze the trade-offs between component choices and their impact on
performance metrics in a multi-dimensional sense. A popular choice of
performance metrics are, namely, inertia matrix which measures kinetic energy
of the entire system [6] and force capability which measures available force
that end-effectors can exert on the environment [7]. In legged robot design,
coupled actuation has shown potential benefits of amplifying joint torques
without increasing reflected inertia extensively [5]. However, under coupling,
modifying even a single gear ratio out of many actuators affects torque and
velocity capabilities of other actuators, which complicates design analysis.
Second, it is difficult to smoothly parameterize characteristics of available
design candidates. Designers typically source motors or reduction mechanisms
from the market, where there can be significant jumps in properties between
available options. This demands designers reflect the discontinuities in
design studies, even complicating the design process.

Figure 1: Conceptual overview of the proposed design optimization of multi-
actuator system. Video: https://youtu.be/DEXYp0TsqRI

Fortunately, a broad spectrum of methods exists to address the challenges
associated with designing multi-DoF systems. One straightforward approach
involves defining high-level tasks and identifying optimal component choices,
such as gear ratios, motors, and limb lengths, with the goal of optimizing
system performance metrics such as power consumption [8, 9, 10] or force
capability [11]. Another family of methods is called the co-design paradigm,
wherein the design and control aspects are jointly optimized through the
addition of a controller and simulation of design cases [12]. This co-design
approach emphasizes control-related aspects such as sensitivity [13] and
robustness [14, 15].

However, existing optimal design tools lack certain crucial aspects of the
design process, often yielding theoretical designs that lack practicality in
real-world applications. First, the expert knowledge of designers, which
includes non-parametric design rules, proves essential in humanoid design. For
example, the actuation system design in humanoids involves intricate space
optimization to accommodate multiple actuators within the robot’s body. These
considerations are better defined as logical constraints rather than
parametric ones. Secondly, designers typically rely on optimal design software
to gain a comprehensive view of multiple design scenarios, rather than seeking
a singular theoretically optimal solution. Moreover, designers engage in an
iterative process of obtaining optimal designs from software and making
necessary modifications. In this process, design tools assist in reasoning and
justifying these modifications. Hence, the software employed must be
relatively fast and capable of incorporating human intuitions. Regrettably,
current state-of-the-art frameworks heavily rely on computationally intensive
algorithms, such as genetic algorithms, or the number of optimization
variables is limited.

This paper proposes a design tool that assists designers in creating high-DoF
actuation systems, addressing the aforementioned issues. First, we want a
design tool that guarantees a system design that is capable of various
humanoid tasks such as carrying objects, weight-lifting, and/or running.
However, it is a chicken-and-egg problem; without a robot design, one cannot
reason about tasks and vice versa. Thus, tasks or motions are planned based on
a reduced-order model (ROM) of a robot consisting of high-level design
parameters such as overall mass, inertia, and center of mass (CoM) height with
fixed morphology and ordering of joints. These assumptions are consistent with
existing motion planners based on single rigid body dynamics (SRBD) [16] or
centroidal dynamics [1], allowing us to use similar ROM-based motion planners.
Second, in this framework, designers are encouraged to use their expert
knowledge to populate and curate libraries of design candidates such as
motors, gear ratios, and coupling mechanisms between actuators and joints.
Moreover, the high-level intuition of designers is translated into design
rules, which designers can experiment with to observe how optimal solutions
behave under these rules. Third, our objective is to create dynamic humanoids
capable of rapid limb acceleration and impact mitigation. To achieve these
goals, [4] suggests that minimizing the reflected inertia of actuation systems
is the key, which sets the objective of the proposed design optimization to
choosing components from design libraries that minimize reflected inertia
while ensuring task completion. Lastly, to solve the optimization more
efficiently, the optimal design problem is translated into a binary integer
linear program which yields optimal solutions faster than existing methods.

This framework is demonstrated by design studies of a biped with 5-DoF legs
focusing on the effects of 1) achieving a set of tasks, 2) employing actuator-
joint couplings (parallel mechanism) and 3) limiting the mass budget of all
motors in a system. The optimization study was given two tasks that involve
different joint usage patterns: walking and weight-lifting in a snatch style.
Next, design libraries of motors, gear ratio, and coupling mechanisms were
compiled from existing humanoids and unprecedented-but-practical designs for
comparison. From the results of 200 optimal design studies with varied motor
mass budgets and usage of couplings, several trends in reflected inertia,
force/torque capability, and copper loss of motors were observed.

The contribution of the proposed framework is as follows. First, this tool
allows designers to quickly navigate vast design spaces and reveals a
comprehensive overview of optimal multi-actuator system designs for legged
robots. Second, the framework serves as a reasoning tool for high-DoF system
design. Designers can quickly run a large number of optimal studies to observe
how optimal solutions respond to modifying tasks, design libraries, or design
rules. For example, providing a larger set of tasks would result in general-
purpose humanoid design whereas providing a single task would yield a more
dedicated system design. Moreover, novice engineers can update components
available in the market to create better designs, while experts can explore
combinations of novel mechanisms for comparison.

##  II Background: Representation of Design Library & Combinations

The binary variables are efficient tools for handling two essential aspects of
design: 1) determining whether to select or exclude component candidates from
libraries and 2) expressing combinations of components.

Definition 1: (Binary Combination) Let
$\bm{\mathrm{p}}\\!=\\!\\{p_{1},\cdots,p_{n}\\}$ be a library of $n$ objects,
$\bm{\mathrm{x}}\\!=\\!(x_{1},\cdots,x_{n})\\!\in\\!\mathbb{\\{}0,1\\}^{n}$ be
a vector of binary variables. An arbitrary singular choice of an object
$p(\bm{\mathrm{x}})_{B}$ out of the library $\bm{\mathrm{p}}$ is obtained by a
binary combination of $\bm{\mathrm{p}}$ and $\bm{\mathrm{x}}$ as follows,

$\displaystyle
p(\bm{\mathrm{x}})_{B}\coloneqq\textstyle\sum_{i=1}^{n}p_{i}x_{i}\>\>\>$ (1)
$\displaystyle\Longrightarrow\>\>\>$
$\displaystyle\textstyle\sum_{i=1}^{n}x_{i}=1\quad(\textrm{exclusivity}).$ (2)

For instance, a selection of $p_{j}$ is achieved by assigning 1 to the
corresponding variable $(x_{j}\\!=\\!1)$ and 0 to the rest
$(x_{k}\\!=\\!0,k\\!\neq\\!j,k\\!\in\\![n])$. The definition itself implies an
exclusivity constraint that enforces a single choice. It should be noted that
the objects $p_{i}$ can be scalars such as maximum torque or matrices like the
Jacobian of coupling mechanism that maps actuator velocities to joint
velocities.

Remark 1: (Linearization of Multilinear Polynomials) Products of binary
variables commonly arise when combining components; e.g., an actuator is a
combination of a motor and a reduction mechanism. A useful property of the
product of binary variables is that it can be perfectly linearized (as opposed
to linearly approximated). The product of binary variables $x$ and $y$,
denoted as $xy$, is linearized by introducing a new binary variable $z$ in
place of $xy$ and adding the following linear constraints:

$z\leq x,\quad z\leq y,\quad z\geq x+y-1.$ (3)

For linearization of general multilinear forms, refer to [17].

##  III Optimal Design of Transmission System for Humanoid Robots

###  III-A Overview

We introduce a design optimization tool that finds mechanical components for
the actuation system of athletic humanoids. Most importantly, this framework
is composed as an interactive tool for designers and addresses difficulties in
designing a high-DoF actuation system. Hence, the framework expects designers
to provide their knowledge and intuition and run several iterations to gain an
understanding of optimal designs of humanoids.

The workflow of the proposed tool is divided into three steps. (Fig. 2). The
first step is to prepare key ingredients for design optimization; 1) ROMs for
the generation of task trajectories, 2) a set of task trajectories generated
by motion planning software using ROMs, 3) libraries of candidate components
to create the robot, 4) relevant design rules, and 5) design goals. As to the
preparation of task trajectories, first, ROMs are created from known
kinematics and parameters of high-level design envelopes and behavior targets
set by designers (e.g., overall mass, inertia, CoM height or desired velocity)
depending on the tasks. Then, using ROM-based motion planning approaches
(e.g., SRBD-based trajectory optimization [18] or kino-dynamic planning [19]),
task-space trajectories of desired tasks are obtained. Finally, the task-space
trajectories are converted to joint-space trajectories using inverse
kinematics, inverse dynamics, or whole-body controllers. It is important to
note that the ROM assumes fixed limb lengths and joint ordering. This is
because humanoids are typically designed to have similar proportions and
kinematics to humans. Also, as there exist only a few practical joint
orderings that mimic human limb kinematics, design studies for every joint
ordering can be carried out. It is also assumed that the actuators are placed
very closely to the center of the robot (proximal actuation) and actuator
dynamics is neglected. Hence, the actuator placements and their inertial
torques do not affect the dynamics of the robot. Lastly, the framework does
not take into account the structural design of limbs because we assume that
the apparent inertia and mass could be significantly minimized compared to the
reflected inertia of the actuation system.

The second step of the framework is to integrate tasks, design library, and
design rules and translate them to linear binary integer programming. The
overarching concept is expressed as follows with $\bm{\mathrm{x}}$ as design
choices of mechanical components, —s— xReflectedInertia(x), x∈DesignLibrary
TaskTrajectory⊆OperationDomain(x) DesignRules(x) . The contents of (III-A) are
briefed as follows.

Control-Awareness (III-A): The objective is to minimize the reflected inertia
of the entire actuation system perceived at the joint-level (joint-space
reflected inertia). This objective aims to improve force control performance,
enhance the system’s ability to regulate external impacts and enable agile
limb movements [4]. The objective aligns with the QDD paradigm, although
designers have the flexibility to modify it to prioritize energy efficiency or
minimum torque expenditure, depending on their interest.

Figure 2: Overview of design optimization of high-DoF actuation system for
humanoid robots

Discrete Design Space (III-A): In this framework, designers are responsible
for providing and modifying libraries of candidate motors, reduction
mechanisms, and actuator-joint couplings. An optimization solver then searches
for the optimal design by exploring all possible combinations of provided
components. Designers are encouraged to use actuator-joint coupling mechanisms
from existing robots, utilize their expert knowledge to curate practical
component candidates, or include novel components whose utility is in
question.

Task-Awareness (III-A): Robots are expected to achieve multiple tasks. This
framework ensures the feasibility of these tasks within the capability of
motors. In other words, the given tasks are mapped to trajectories of motor
torques and velocities. Then, these trajectories are constrained to remain
within the motor limits such as voltage and thermal limit (operation domain)
as in Fig. 3(b). The safety of a controller can be tuned by parameters that
define allowable closeness between extremes of task trajectory and the
boundary of the operating domain.

Design Rules (III-A): Additional constraints can be incorporated into the
optimization process, such as limiting the sum of motor mass for a subset of
joints or limiting the copper loss on motors for a certain duration. As long
as the constraints are multilinear polynomial, they can be perfectly
linearized (see Sec. II). Although not presented in this paper, logical
constraints can be imposed for compact arrangement of large components such as
motors.

As the last step of the workflow, designers need to analyze the resulting
optimal solutions and modify the aforementioned ingredients for design
optimization. For instance, if the optimization is infeasible due to high
torque requirements, one can expand the design library by adding more torque-
dense motors, relax the design rules, or relax the task requirements (e.g.,
smaller total mass, smaller desired CoM velocity of the ROM).

###  III-B Maximizing Force Control Bandwidth

Smaller reflected inertia of actuators contributes to faster limb swing or
acceleration and improves the ability to mitigate impacts from contact by
quickly retracting from contacts (e.g., ground contact). Both capabilities are
critical requirements for dynamic humanoids.

Hereafter, the goal of minimizing reflected inertia of the actuation system is
reconnected with the earlier optimization sketch (III-A) to formulate
objective functions using binary variables. Although [4] suggests minimizing
reflected inertia in task-space, obtaining such expression using binary
variables in linearizable form is prohibitive due to rational forms arising
from matrix inversions [6, Eqn. (17)]. As a proxy to task-space inertia, a
joint-space reflected inertia matrix is used. The downside of using this proxy
is that it neglects the lever effect of limb lengths to joints which is beyond
the scope of this paper. Nevertheless, some relevant observations are stated
in Sec. V-C.

Let us build the joint-space reflected inertia matrix of an actuation system
consisting of $n_{d}$ joints (or actuators) from libraries of motors,
gearboxes, and actuator-joint couplings with associated binary variables
$\bm{\mathrm{x}}_{i}\\!\in\\!\mathbb{B}^{n_{{m}_{i}}}$,
$\bm{\mathrm{y}}_{i}\\!\in\\!\mathbb{B}^{n_{g_{i}}}$, $i\\!\in\\![n_{d}]$, and
$\bm{\mathrm{z}}\\!\in\\!\mathbb{B}^{n_{c}}$, where $n_{(\cdot)}$ is the
number of candidates in the respective libraries. The reflected inertia of an
$i$’th actuator is the product of binary combinations of rotor inertia
$I_{i}^{r}(\bm{\mathrm{x}}_{i})_{B}$ and gear ratio squared
$N^{2}_{i}(\bm{\mathrm{y}}_{i})_{B}$,

$I_{i}^{a}(\bm{\mathrm{x}}_{i},\bm{\mathrm{y}}_{i})=I^{r}_{i}(\bm{\mathrm{x}}_{i})_{B}N_{i}^{2}(\bm{\mathrm{y}}_{i})_{B}.$
(4)

Next, the reflected inertia of each actuator forms a diagonal matrix of
actuator-space reflected inertia
$\bm{\mathrm{H}}_{\textrm{r}}^{a}(\bm{\mathrm{x}},\bm{\mathrm{y}})=\operatorname{diag}({I}^{a}_{1},\dots,{I}^{a}_{n_{d}})$.
This inertia matrix is projected on to joint-space by a binary combination of
actuator-joint coupling Jacobian $\bm{\mathrm{C}}(\bm{\mathrm{z}})_{B}$ which
represents a mechanical connection between actuators and joints. The coupling
Jacobian maps actuator velocities
$\dot{\bm{\mathrm{\psi}}}\in\mathbb{R}^{n_{a}}$ to joint velocities
$\dot{\bm{\mathrm{q}}}\in\mathbb{R}^{n_{a}}$ as
$\bm{\mathrm{C}}:\dot{\bm{\mathrm{\psi}}}\rightarrow\dot{\bm{\mathrm{q}}}$.
The joint-space reflected inertia is as follows,

$\displaystyle\overline{\bm{\mathrm{H}}}^{\mathrm{j}}_{r}(\bm{\mathrm{x}},\bm{\mathrm{y}},\bm{\mathrm{z}},k)=\bm{\mathrm{C}}^{-\top}\bm{\mathrm{H}}_{\textrm{r}}^{a}\bm{\mathrm{C}}^{-1},$
(5)

where $k$ denotes the time instance at which the inertia is evaluated along
the given tasks and $\bm{\mathrm{x}}$, $\bm{\mathrm{y}}$ are lifted binary
variables of $\bm{\mathrm{x}}_{i}$, $\bm{\mathrm{y}}_{i}$.

Lastly, the size of the inertia matrix is measured using the
$\operatorname{trace}()$ operator. The trace of an inertia matrix is the sum
of its eigenvalues which relate to the lengths of semi-axes of the inertia
ellipsoid (Fig. 3(a)). Hence, the trace operator allows the optimizer to
minimize reflected inertia on every dimension until motor mass becomes too
small and, subsequently, a motor does not have sufficient torque to achieve
given tasks.

Figure 3: (a) The size of reflected inertia matrix is measured as the sum of
semi-axes lengths which is obtained by trace operation (b) Task trajectory
projected on to motor-level is constrained to remain within the operation
domain of a motor defined by vertices $\bm{\mathrm{v}}_{i}$.

###  III-C Task Feasibility

Task feasibility assures that the desired task trajectories are accommodated
within the capabilities of the motors. The torque and velocity capability or
an operation domain of a motor is approximated as a convex polygon
$\operatorname{convhull}(\bm{\mathrm{v}}_{1},\dots,\bm{\mathrm{v}}_{n})$
(colored polygon in Fig. 3(b)). This polygon comprises vertices
$\bm{\mathrm{v}}_{i}=(\alpha_{i}\omega_{\circ},\beta_{i}\tau_{\circ})$, whose
coordinates are defined by motor parameters $\omega_{\circ},\tau_{\circ}$ such
as rated velocity and torque with coefficients $\alpha_{i},\beta_{i}$. Next,
the vertices are associated with the motor library through a binary
combination,
$\bm{\mathrm{v}}_{i}(\bm{\mathrm{x}})=(\alpha_{i}\omega_{\circ}(\bm{\mathrm{x}})_{B},\beta_{i}\tau_{\circ}(\bm{\mathrm{x}})_{B})$.
Consequently, the operation domain is expressed as follows,

$P(\bm{\mathrm{x}},\bm{\mathrm{\alpha}},\bm{\mathrm{\beta}})=\operatorname{convhull}\big{(}\bm{\mathrm{v}}_{1}(\bm{\mathrm{x}})_{B},\dots,\bm{\mathrm{v}}_{k}(\bm{\mathrm{x}})_{B}\big{)},$
(6)

where $\bm{\mathrm{\alpha}}$ and $\bm{\mathrm{\beta}}$ are vectors of
coefficients $\alpha_{(\cdot)}$ and $\beta_{(\cdot)}$.

Next, the joint-space trajectories of given tasks are concatenated and
projected to motor-space. The joint-space trajectory of length $T$,
$U^{\mathrm{j}}=\\{(\bm{\mathrm{\omega}},\bm{\mathrm{\tau}})_{k}\>|\>\bm{\mathrm{\omega}},\bm{\mathrm{\tau}}\in\mathbb{R}^{m},k\in[T]\\}$,
is projected to actuator-space using actuator-joint coupling Jacobian
$\bm{\mathrm{C}}(\bm{\mathrm{z}})_{B}$, then projected to motor-space using
gear ratio Jacobian
$\bm{\mathrm{N}}(\bm{\mathrm{y}})=\operatorname{diag}(N_{1}(\bm{\mathrm{y}})_{B},\dots,N_{m}(\bm{\mathrm{y}})_{B})$.
The tasks trajectory on motor-space
$U^{\textrm{m}}(\bm{\mathrm{y}},\bm{\mathrm{z}})$ is obtained as a function of
choices of gearboxes and couplings,

$\displaystyle U^{\mathrm{m}}(\bm{\mathrm{y}},\bm{\mathrm{z}})=\\{($
$\displaystyle\bm{\mathrm{\omega}},\bm{\mathrm{\tau}})\>|\>\bm{\mathrm{\omega}}=\bm{\mathrm{N}}\bm{\mathrm{C}}^{-1}\bm{\mathrm{\omega}}_{\mathrm{j}},$
(7)
$\displaystyle\bm{\mathrm{\tau}}=\bm{\mathrm{N}}^{-\top}\bm{\mathrm{C}}^{\top}\bm{\mathrm{\tau}}_{\mathrm{j}},(\bm{\mathrm{\omega}}_{\mathrm{j}},\bm{\mathrm{\tau}}_{\mathrm{j}})\in
U^{\mathrm{j}}\\}.$

Finally, the following task feasibility constraint ensures the task trajectory
to be contained within the operation domain

$U^{m}(\bm{\mathrm{y}},\bm{\mathrm{z}})\subseteq
P(\bm{\mathrm{x}},\bm{\mathrm{\alpha}},\bm{\mathrm{\beta}}),$ (8)

which easily converts to a set of linear constraints.

###  III-D Design Rules

In this study, we bound the motor mass budget. If the motor mass budget is
left unbounded, the optimizer will always prefer heavier, larger, and more
torque-dense motors. However, there is a limit to the diameter of the motor
that fits inside the robot’s body. Moreover, the sum of motor mass may take up
to $40\%$ of the entire mass of a robot [20]. Hence, a smaller motor mass
budget may leave significant room for a larger payload. The motor mass budget
constraint is enforced as follows,

$m_{t}(\bm{\mathrm{x}})_{B}=\textstyle\sum_{i=1}^{m}m_{i}x_{i}\leq
m_{t}^{\textrm{UB}},$ (9)

where $m_{t}$ is total mass of a design case, $m_{i}$ is mass of candidate
motors and $m_{t}^{\textrm{UB}}$ is the total mass budget.

Designers can further employ various design rules such as limits on torque,
power, etc. at various levels such as motor, actuator, and joint-levels. It is
because system capability and task trajectories can be explicitly expressed
using binary combinations. Second, logical constraints are also allowed. A
statement, ’If a component $A_{1}$ is used, then $A_{2}$ also is used’, is
translated as $x_{1}-x_{2}\geq 0$ where $x_{1}$, $x_{2}$ are binary variables
that associate components $A_{1}$, $A_{2}$. For more logical constraints,
refer to [21].

##  IV Humanoid Design Study

###  IV-A Robot Model and Tasks

A humanoid robot is modeled as a single rigid body with massless limbs and
weight $20$ kg that closely resembles Tello in [5]. The joint ordering is hip
yaw, hip roll, hip pitch, and knee pitch followed by ankle pitch.

Next, based on the ROM, two task trajectories were generated (Fig. 4). The
first task corresponds to moderately fast walking, with a peak speed of 1.3
m/s, a step frequency of 0.3 s. Joint torques and velocities of single support
phase (0.15 s) are used. The walking trajectories are generated by simulating
the SRBD of Tello. A step placement feedback controller that tracks a
predefined walking speed is implemented [22]. An optimization-based balance
controller, similar to in [23], computes the ground reaction forces required
to realize the desired center-of-mass walking dynamics. The SRBD task-space
trajectories are converted to joint-space trajectories using inverse
kinematics. The resulting joint velocities, which contain high-frequency
signals, are low-pass filtered at a cutoff frequency of 30 Hz to obtain
meaningful solutions. The second task trajectory is to lift a point mass of 15
kg in a snatch style over 3.6 s. This trajectory is generated by extracting
the joint trajectory from a video of an athlete performing a snatch. The
trajectory is scaled down to match the size of the robot, and joint torques
are obtained using inverse dynamics. The lifting task only recruits the pitch
axes of the hip, knee, and ankle joints.

It is important to note that designers may choose other tasks and this study
is merely showcasing an example. Depending on the choice of task, optimal
designs may vary significantly.

###  IV-B Design Library and Rules

The design space is formulated by combinations of the following component
libraries. The motor library contains $n_{m}=10$ motors from CubeMars (AK, R,
RI series [24]). The reduction library consists of ideal transmission without
frictional loss and gear ratio ranges from 1 to 12 at whole number interval
$(n_{r}=12)$. The vertices of the operation domains of the motors are defined
by peak torque and rated velocity from datasheets and vectors of tuning
coefficients
$\bm{\mathrm{\alpha}}=(1.5,-2.5,-2.5,1.5),\,\bm{\mathrm{\beta}}=(1.2,1.2,-1.2,-1.2)$.

(a)

(b)

Figure 4: Sketches of tasks used in humanoid design study. (a) walking forward at moderately fast speed and (b) lifting an object in snatch style.  TABLE I: List of Actuator-Joint Couplings Name | Jacobian |  | Name | Jacobian
---|---|---|---|---
serial | $\operatorname{diag}({1,1,1,1,1})$ |  | par-45 | $\operatorname{diag}(1,1,1,\bm{\mathrm{C}}_{d})$
par-12 | $\operatorname{diag}(\bm{\mathrm{C}}_{d},1,1,1)$ |  | par-12-34 | $\operatorname{diag}(\bm{\mathrm{C}}_{d},\bm{\mathrm{C}}_{d},1)$
par-23 | $\operatorname{diag}(1,\bm{\mathrm{C}}_{d},1,1)$ |  | par-12-45 | $\operatorname{diag}(\bm{\mathrm{C}}_{d},1,\bm{\mathrm{C}}_{d})$
par-34 | $\operatorname{diag}(1,1,\bm{\mathrm{C}}_{d},1)$ |  | par-23-45 | $\operatorname{diag}(1,\bm{\mathrm{C}}_{d},\bm{\mathrm{C}}_{d})$

There are 8 couplings studied in this paper including 1 serial and 7 parallel
couplings. The parallel couplings are based on modules of 2-DoF differential
coupling whose Jacobian is

$\bm{\mathrm{C}}_{d}=\frac{1}{2}\begin{bmatrix}1&1\\\ -1&1\end{bmatrix}.$ (10)

The couplings are labeled as ‘serial’, ’par-($n_{1}m_{1}$)’ and
‘par-($n_{1}m_{1}$)-($n_{2}m_{2}$)’ where $n_{i},m_{i}$ denotes the index of
joints involved in $i$’th differential coupling (see Table.
LABEL:table:Jacobians). Since a mechanical coupling between faraway joints is
typically impractical, the couplings are limited to consecutive joints. The
hardware realization of Tello [5] utilizes ‘par-23-45’.

The total mass of 5 motors is subject to a constraint (motor mass budget),
with 25 values ranging from 2.2 kg to 4.0 kg at intervals of 0.075 kg. After a
few runs of modifying the values of the motor mass budget, it was found that
the optimization was infeasible below 2.2 kg and solutions did not vary beyond
4.0 kg.

###  IV-C Design Optimization

The optimization problem is detailed as follows. Inclusion of library data is
implied by the use of binary optimization variables
$\bm{\mathrm{x}},\bm{\mathrm{y}},\bm{\mathrm{z}}$. —s— x, y, z ∑k=1Ttrace(¯Hjr
(x,y, z, k)) U^m(y,z) ⊆P(x, α, β) m_t(x)_B ≤m_t^UB x∈B^n_a ×n_m, y∈B^n_a ×n_r,
z∈B^n_c exclusivity(x,y,z) To explicitly study the effect of constraining
total mass budget for all 5 motors and employing different coupling Jacobian,
the optimization was run for all combination of mass budget constraints and
coupling Jacobians (total $25\times 8=200$ runs). After optimal solutions are
obtained, system characteristics such as inertia matrix and force/torque
capability polytopes (FCP/TCP) in joint and task-space and copper loss on
actuators are retrieved. The FCP/TCP represents the maximum force/torque that
can be applied for a given configuration, based on the limits of actuators
[7].

The optimization problem (IV-C) is linearized using techniques introduced in
Sec. II. The design optimization is written in MATLAB and calls an integer
linear programming solver (intlinprog, Gurobi). The solver was simultaneously
run on CPU (20 cores, Intel i7-12700H) with default solver parameters.
Computation of total 200 runs took less than 100 s.

##  V Results and Discussion

(a)

(b)

Figure 5: (a) Optimal costs (trace of joint-space reflected inertia) vs. total
mass budget for 5 motors for every coupling case. Total copper losses per
walking and snatch lifting task (b) Optimal designs’ motor mass budget,
optimal cost and copper loss of snatch task. The bigger dots represent Pareto
front.  Figure 6: Comparison of two optimal designs under identical motor mass
budget and with different usage of coupling mechanism. (a1), (a2) FCP of
contact leg in contact along walking. (b1), (b2) Frontal views show
significant difference in size the FCP. (c1), (c2) FCPs of both designs are
large enough to contain ground reaction forces to achieve moderately fast
walking. (d1), (d2) Couplings affect how task trajectories are projected onto
motor-level and motor torque capability margins (T1 vs. T2).

###  V-A Motor Mass Budget vs. Reflected Inertia vs. Copper Loss

The first row of Fig. 5(a) gives an overview of optimal designs obtained with
varied motor mass budgets and use of different couplings. First, as the motor
mass budget decreases, availability of larger and more torque-dense motors
decreases. Consequently, the optimal solutions leans towards smaller and
weaker motors, necessitating the usage of elevated gear ratios. In the end,
the use of higher gear ratios inflates reflected inertia or contributes to
higher cost (optimal objective function value).

As to the effect of couplings, a few coupling cases were found infeasible
(par-34, par-12-34) across all motor mass budgets, whereas other couplings
become infeasible as the motor mass budget decreases. When the motor mass
budget was limited to extremely small values, only coupled mechanisms provided
sufficient torques and velocities (feasible). The serial actuation (no
coupling) demonstrated the least reflected inertia at the cost of high copper
loss for walking task across all cases of motor mass budget. In conclusion,
this result suggests that 1) the framework can filter out couplings that
cannot generate sufficient torque or speed, 2) at extremely low motor mass
budget, actuators need to cooperate so that they can amplify joint torques to
achieve desired tasks, and 3) use of coupling increases the reflected inertia
for designs with a large motor mass budget.

The second and third row of Fig. 5(a) shows the combined copper loss of all
motors of each task. As the motor mass budget decreases, the copper loss also
decreases for most couplings. This tendency is opposite to that of reflected
inertia because larger reflected inertia implies higher gear ratios, which
again implies less motor torque and motor current draw. Two optimal designs
following this trend are compared in Table. LABEL:table:DesignComparison.
However, these implications are not strict: in the walking task, a few
coupling cases did not follow this tendency (Fig. 5a, second row, serial,
par-45, par-23-45, around 3.2$\sim$3.7 kg). This is because different choices
of gear ratio and motor can emerge as long as they are optimal and satisfy
design rules.

If this study is posed as a multi-objective optimization with minimization of
reflected inertia and the copper loss, the optimal solutions can be plotted in
3D space as Fig 5(b). Pareto optimal solutions are emphasized as larger dots
over other optimal solutions. To aid visual understanding, optimal solutions
are projected onto 2D slices.

###  V-B Effect of Actuator-Joint Coupling

The influence of the usage of couplings on task-space FCP and motor-level
trajectory is visualized in Fig. 6 which presents a comparative analysis
between two design cases, specifically, ‘par-12’ and ‘par-23-45’. The latter
utilizes more a higher degree of coupling than the former while both are
designed with an identical motor mass budget.

A notable point in this comparison involves the velocity-dependent nature of
FCP, attributed to the contour of the operational domain of motors. When motor
velocities are high, the available torque (T1) near the boundary of the
operation domain is reduced due to high motor velocity (Fig. 6(d1), (d2)).

The latter case, with additional degrees of coupling, has more chances for
cooperative actuation which means that the load on a joint can be shared among
the coupled actuators. Thanks to the load sharing, the individual actuators
experience lower torque requirements. This allows employing smaller gear
ratios which renders smaller motor velocities. Hence, employing more degrees
of coupling has a benefit of keeping the motor-level trajectory (velocity) at
a safe distance from the boundary of the operational domain. The degradation
of motor torque capability in design case ‘par-12’ manifests as a
significantly smaller FCP at 129 ms (Fig. 6(a1), (b1)), in contrast to a
larger FCP in (Fig. 6(a2), (b2)).

The state-dependent nature of the FCP could be related to the margin of
admissible error of feedback controllers. For instance, if joint velocity
deviates from a (pre)planned trajectory towards the boundary of the operation
domain of motors, the torque capability diminishes. Simultaneously, a feedback
controller will attempt to bring the state back to reference. At that moment,
if the torque margin is insufficient, the torque commands from the feedback
controller could suffer from saturation which may lead to failure. Designers
can prevent this by modifying the approximated operation domain to be smaller
than the actual operation domain of motors by tuning the coefficients
$\bm{\mathrm{\alpha}}$ and $\bm{\mathrm{\beta}}$.

TABLE II: Comparison of Two Optimal Design Cases

| Joint (Axis)
---|---
|  Hip(z) |  Hip(x) |  Hip(y) |  Knee(y) |  Ankle(y)
Design Case A (lightest, coupling $=$ (par-12), $m_{t}=2.28$ kg)
Inertia [$\textrm{gm}^{2}$] | ( 0.31, 0.31 ) |  1.06 |  1.06 |  0.75
Peak Torq. [Nm] | ( 25.7, 25.7 ) |  92.4 |  92.4 |  28.9
Motor | ( RI70, RI70 ) |  R100 |  R100 |  RI70
Gear Ratio | ( 3, 5 ) |  7 |  7 |  9
Design Case B (heaviest, coupling $=$ (par-23-45), $m_{t}=3.7$ kg)
Inertia [$\textrm{gm}^{2}$] |  0.068 | ( 0.69, 0.69 ) | ( 0.73, 0.73 )
Peak Torq. [Nm] |  26.4 | ( 105.6, 105.6 ) | ( 105.6, 105.6 )
Motor |  R100 | ( R100, R100 ) | ( R100, R100 )
Gear Ratio |  2 | ( 4, 4 ) | ( 3, 5 )

  * •

Note: Two entities within a parenthesis denotes two joints involved in 2-DoF
actuator-joint coupling.

###  V-C Characterization of Inertia and Capability Polytopes

The relation between reflected inertia and FCP/TCP in both joint- and task-
space is analyzed to address two possible claims; 1) minimization of reflected
inertia should take place at task-space level, not joint-space and 2)
maximization of TCP/FCP could be also taken into account as well as
minimization of the reflected inertia. Unfortunately, the authors could not
obtain direct quantification of reflected inertia in task-space and the volume
of TCP/FCP using binary integer linear programming. Still, a few trends were
observed by plotting the relation between the size of reflected inertia versus
volume of TCP/FCP in joint-space and task-space as in Fig. 7. Optimal designs
with smaller motor mass budgets are plotted with more opacity. The TCP/FCP
were calculated without considering the degradation in torque margin due to
high velocity. Since FCP varies as configuration (joint angles) changes, its
distribution is illustrated as an 0.7-sigma ellipse. In each space, two
similar trends were observed; 1) serial actuation shows more cohesiveness in
the relation between the size of reflected inertia and FCP/TCP compared to
coupled ones, 2) a combination of coupled actuation and smaller mass budget
leads to larger size and variation of reflected inertia, and 3) larger motor
mass budget allows larger FCP/TCP.

Figure 7: Normalized optimal objective function value (trace of reflected
inertia matrix) vs. normalized volume of FCP/TCP in joint-space (left) and
task-space (right).

##  VI Conclusion

The existing design principles and optimal design tools for the design of
multi-DoF actuation systems often fall short in providing practical designs or
understanding of optimal designs. To resolve this issue, an interactive
optimal design tool for actuator system for humanoids is introduced. This tool
is designed to incorporate expert designers’ knowledge and rapidly solve large
numbers of design studies to expand understanding of optimal designs.

This paper closes with several limitations of the present study and
suggestions for future research. First, ROM employed in this study assumes
massless limbs and actuator dynamics is neglected. Hence, the inertial forces
exerted on the limbs and rotors within the actuation system are not accounted
for. Future work may incorporate inertial forces associated with the limbs’
apparent inertia and the reflected inertia of the actuators into the
optimization process. Second, the applicability of optimal design solutions
produced by the proposed tool remains unclear. An essential next step would be
to integrate this tool with simulation software for the validation of proposed
designs.

## ACKNOWLEDGMENT

To research advisors Patrick Wensing and Donghyun Kim, Johannes Englsberger,
and my colleagues Donghoon Baek and Seunghyun Bang for guidance and fruitful
discussion.

## References

  * [1] M. Chignoli, D. Kim, E. Stanger-Jones, and S. Kim, “The mit humanoid robot: Design, motion planning, and control for acrobatic behaviors,” in _2020 IEEE-RAS 20th International Conference on Humanoid Robots (Humanoids)_ , 2021, pp. 1–8.
  * [2] S. H. Jeon, S. Kim, and D. Kim, “Online optimal landing control of the mit mini cheetah,” in _2022 International Conference on Robotics and Automation (ICRA)_ , 2022, pp. 178–184.
  * [3] M. S. Ahn, _Development and Real-Time Optimization-based Control of a Full-sized Humanoid for Dynamic Walking and Running_. University of California, Los Angeles, 2023.
  * [4] P. M. Wensing, A. Wang, S. Seok, D. Otten, J. Lang, and S. Kim, “Proprioceptive actuator design in the mit cheetah: Impact mitigation and high-bandwidth physical interaction for dynamic legged robots,” _IEEE Transactions on Robotics_ , vol. 33, no. 3, pp. 509–522, 2017.
  * [5] Y. Sim and J. Ramos, “Tello leg: The study of design principles and metrics for dynamic humanoid robots,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 4, pp. 9318–9325, 2022.
  * [6] O. Khatib, “Inertial properties in robotic manipulation: An object-level framework,” _The international journal of robotics research_ , vol. 14, no. 1, pp. 19–36, 1995.
  * [7] P. Chiacchio, Y. Bouffard-Vercelli, and F. Pierrot, “Force polytope and force ellipsoid for redundant manipulators,” _Journal of Robotic Systems_ , vol. 14, no. 8, pp. 613–620, 1997.
  * [8] Y.-H. Shin, S. Hong, S. Woo, J. Choe, H. Son, G. Kim, J.-H. Kim, K. Lee, J. Hwangbo, and H.-W. Park, “Design of kaist hound, a quadruped robot platform for fast and efficient locomotion with mixed-integer nonlinear optimization of a gear train,” in _2022 International Conference on Robotics and Automation (ICRA)_ , 2022, pp. 6614–6620.
  * [9] M. Chadwick, H. Kolvenbach, F. Dubois, H. F. Lau, and M. Hutter, “Vitruvio: An open-source leg design optimization toolbox for walking robots,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 6318–6325, 2020.
  * [10] W. Pawlus, G. Hovland, M. Choux, D. Frick, and M. Morari, “Drivetrain design optimization for electrically actuated systems via mixed integer programing,” in _IECON 2015-41st Annual Conference of the IEEE Industrial Electronics Society_. IEEE, 2015, pp. 001 465–001 470.
  * [11] A. Hassan and M. Abomoharam, “Modeling and design optimization of a robot gripper mechanism,” _Robotics and Computer-Integrated Manufacturing_ , vol. 46, pp. 94–103, 2017.
  * [12] M. Garcia-Sanz, “Control co-design: An engineering game changer,” _Advanced Control for Applications_ , vol. 1, no. 1, p. e18, 2019, e18 ADC2-19-0039. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/adc2.18
  * [13] S. Ha, S. Coros, A. Alspach, J. Kim, and K. Yamane, “Joint optimization of robot design and motion parameters using the implicit function theorem.” in _Robotics: Science and systems_ , vol. 8, 2017.
  * [14] G. Bravo-Palacios and P. M. Wensing, “Large-scale admm-based co-design of legged robots,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2022, pp. 8842–8849.
  * [15] K. M. Digumarti, C. Gehring, S. Coros, J. Hwangbo, and R. Siegwart, “Concurrent optimization of mechanical design and locomotion control of a legged robot,” in _Mobile Service Robotics_. World Scientific, 2014, pp. 315–323.
  * [16] R. Batke, F. Yu, J. Dao, J. Hurst, R. L. Hatton, A. Fern, and K. Green, “Optimizing bipedal maneuvers of single rigid-body models for reinforcement learning,” in _2022 IEEE-RAS 21st International Conference on Humanoid Robots (Humanoids)_. IEEE, 2022, pp. 714–721.
  * [17] C.-T. Chang and C.-C. Chang, “A linearization method for mixed 0–1 polynomial programs,” _Computers & Operations Research_, vol. 27, no. 10, pp. 1005–1016, 2000.
  * [18] J. Ahn, S. J. Jorgensen, S. H. Bang, and L. Sentis, “Versatile locomotion planning and control for humanoid robots,” _Frontiers in Robotics and AI_ , p. 257, 2021.
  * [19] M. Chignoli, D. Kim, E. Stanger-Jones, and S. Kim, “The mit humanoid robot: Design, motion planning, and control for acrobatic behaviors,” in _2020 IEEE-RAS 20th International Conference on Humanoid Robots (Humanoids)_. IEEE, 2021, pp. 1–8.
  * [20] G. Kenneally, A. De, and D. E. Koditschek, “Design principles for a family of direct-drive legged robots,” _IEEE Robotics and Automation Letters_ , vol. 1, no. 2, pp. 900–907, 2016.
  * [21] A. Schrijver, _Theory of linear and integer programming_. John Wiley & Sons, 1998.
  * [22] Y. Gong and J. Grizzle, “One-step ahead prediction of angular momentum about the contact point for control of bipedal locomotion: Validation in a lip-inspired controller,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 2832–2838.
  * [23] M. Chignoli and P. M. Wensing, “Variational-based optimal control of underactuated balancing for dynamic quadrupeds,” _IEEE Access_ , vol. 8, pp. 49 785–49 797, 2020.
  * [24] “Cubemars,” accessed: 12-Jul.-2023. [Online]. Available: https://store.cubemars.com/category-116-Exoskeleton+Robot.html