# Inverse design of arbitrary optical helicity patterns

Romuald Kilianski<EMAIL_ADDRESS>Robert Bennett  School of
Physics & Astronomy, University of Glasgow, Glasgow G12 8QQ, United Kingdom

###### Abstract

Superposing multiple plane waves can generate helicity lattices in which the
optical helicity varies regularly in space. Here we propose an inverse design
method for constructing arbitrary helicity structures based on placing a
digital object consisting of dielectric inclusions in three-dimensional space.
We apply the method to design structures that reproduce two-dimensional
lattices embedded within a three-dimensional region using only a single plane
wave as an input. In order to demonstrate the power and flexibility of our
method, we go beyond the paradigm of a regular lattice and propose structures
which can create arbitrary images consisting of regions of varying helicity,
again using only a single plane wave as an input.

##  I Introduction

The concept of chirality appears across the sciences and details the
characteristics of geometrical arrangements with identical constituents but
inverted, non-superimposable forms. Circularly polarised light fits the
criteria of a chiral object as spatial inversion transforms left-handed
circularly polarised light into its right-handed form — the chiral structure
is the combination of its propagation direction and polarisation. This
property makes light with circular polarisation exert influence on other
chiral objects with which it interacts. For example, the dependence of a
chiral molecule’s absorption cross section on the type of circularly polarised
light that illuminates it is described by an effect known as circular
dichroism (see, for example, [1]). In essence, the rate of excitation is
sensitive to the molecule’s handedness and that of the illuminating light. The
degree to which an electromagnetic field influences such an interaction can be
related to a product of the field vectors, a parity-odd scalar density called
optical chirality [2] (originally introduced by Lipkin as the scalar component
of the ‘00-zilch’ [3], but not assigned a physical interpretation at that
early stage). In our example, we are interested in monochromatic fields. For
this case, the optical chirality is equal — up to a constant [4] — to a
quantity called _electromagnetic helicity_ [5]. Familiar from particle
physics, as well as analogous quantities in plasma physics [6], helicity
describes the projection of the spin angular momentum onto the direction of
propagation. Our specific focus is the electromagnetic helicity _density_
(see, e.g., $\mathcal{H}$ [7]), which for a monochromatic plane wave of
frequency $\omega$ reads

$\displaystyle\mathcal{H}=-\dfrac{i}{4c\omega}\left(\mathbf{E}\cdot\mathbf{H^{*}}-\mathbf{E^{*}}\\!\\!\cdot\mathbf{H}\right),$
(1)

where $\mathbf{E}$ and $\mathbf{H}$ are the complex electric and magnetic
fields, respectively, and $c$ is the speed of light.

###  I.1 Helicity lattices

We will consider circumstances where the helicity density given by Eq. (1)
varies regularly in space and forms a so-called helicity lattice [8], a
periodic structure resembling the arrangement of intensity patterns in optical
lattices. Helicity lattices are defined by having a helicity density that
varies in space but an electric field intensity that remains constant — the
latter feature is of particular importance when considering the interactions
of chiral matter with the lattice. A chiral optical force mediated by the
helicity gradients in the field can accelerate enantiomers (molecular species
of opposite handedness) in different directions [9], providing a means of
separating them. Importantly, the homogeneity of the electric field intensity
in a helicity lattice means that achiral optical forces are (in principle)
completely suppressed, allowing the much weaker chirality-dependent forces to
be isolated. Oppositely handed species often exhibit drastically different
properties, so that methods for efficiently separating and distinguishing them
are an active area of research — for a comprehensive review, see [10].
Helicity lattices are already finding their place in this field, for example a
scheme was recently proposed to utilise helicity lattices in separating cold,
chiral molecules via their chirality-sensitive quantum phase transitions [11].

With their potential in industrial applications, the study of helicity
lattices is ready to move from theoretical curiosity towards an exploration of
more practical methods for their generation. In their original formulation
[8], helicity lattices were discovered by superimposing up to six idealised
plane waves with particular well-defined polarisations. Aside from the
impossibility of generating a true plane wave in the laboratory, the need to
align and manipulate up to six differently-polarised beams is a significant
technical challenge — our goal in this work is to define a proof-of-principle
method for generating helicity lattices in a more robust way. The problem is
similar to that encountered in generating magneto-optical trapping potentials,
where again up to six beams can be required. A way of circumventing this issue
in that context has been to use reflective or diffractive optical elements to
generate the required number of beams in an integrated and compact apparatus
(see, for example, [12, 13, 14]). While it may indeed be possible to follow
similar design principles to create a helicity lattice, the additional
complexity found therein (especially regarding polarisation-sensitivity) leads
us to take a different approach, employing the techniques of _inverse design_.

###  I.2 Inverse design

Traditional ‘forward’ design of optical components relies on a human designer
using some combination of general principles, intuition and previous
experience to produce a candidate design which is then tested and refined
against some desired figure of merit. Inverse design turns this paradigm
around — only the figure of merit is specified, following which a computer
algorithmically builds up an appropriate device. The vast space of possible
designs originally necessitated some imposition of geometrical constraints in
early work in this direction [15, 16], for example the sizes and positions of
fixed circles and rectangles would be altered such that an improved design is
reached. The introduction to photonics of _adjoint methods_ from aerospace
engineering [17], coupled with increased computing power, allowed such
constraints to be lifted, kicking off the era of completely free-form photonic
inverse design [18, 19]. The combination of this and the relatively recent
ability to manipulate matter on the micro- and nano-scales has initiated an
explosion in the application of inverse design methods (for a comprehensive
review, see [20]). In contrast to the earlier constrained approaches, the pool
of available designs is unlimited — the shapes of the structures are dictated
by their impact on performance (and manufacturability) and are no longer
guided by intuition or arbitrary constraints. Some of the more recent diverse
applications of these ideas have been concerned with the conversion of solar
energy [21], integrated photonic devices [22], sub-wavelength focusing [23]
and heat transfer [24]. Most recently, ID has been applied to atomic coherence
[25], light-matter interactions [26] and nuclear quantum optics [27].

As discussed in the previous sections, it has been shown analytically in Ref.
[8] that helicity lattices can be constructed from up to six idealised plane
waves in various configurations, and that we would like to come up with a way
of relaxing this requirement. To accomplish this, we will use a freeform
inverse design approach to generate a dielectric structure that can refract a
single beam multiple times in such a way that a lattice structure emerges in a
given plane. The placement of the elements forming such a structure will turn
out to be non-obvious and not easily inferred from a wave analysis point of
view. The desired outcome will be chosen, and an algorithm will adjust the
geometry to bring a chosen merit function as close as possible to a designated
goal.

This paper is organised as follows: in Section II, we will set up the problem
of helicity lattice generation in terms of an optimisation task, develop the
mathematical structure and propose an algorithm capable of reproducing a
helicity lattice. In Section III we will then outline the specifics of the
actual inverse design process. Subsequently, in Section IV we will present the
results of two single-beam simulations reproducing lattices originally formed
by a superposition of three sources. Additionally, in that section we present
a novel approach where we simulate a helicity density that is not periodic and
can take an arbitrary (user-defined) shape. Critical analysis of the method
and conclusions follow in sections and V and VI.

##  II Description of the problem

Figure 1: The schematic of different stages of building the dielectric
structure using a single beam source. As it grows in size and complexity, the
more the underlying field approximates the helicity lattice.

We aim to produce a two-dimensional optical pattern resembling a helicity
lattice as closely as possible. Instead of superposing multiple light sources,
we want to effect it by shining a single beam onto a crystal-like (not
uniformly periodic) structure that will refract the beam in such a way to
ensure that a helicity lattice emerges. Our idea is to create such a structure
algorithmically; instead of iterating over possible geometries, we let an
algorithm grow the shape of the crystal. This idea is shown schematically in
Fig. 1, where we see the underlying field progressively changing and starting
to resemble the target lattice as more inclusions are added. At every step,
the algorithm will decide on the optimal placement of a dielectric inclusion —
we will use an inverse design algorithm based on the adjoint method.

The adjoint method relies on the fact that Maxwell’s equations in (reciprocal)
media are symmetrical with respect to the exchange of source and observation
points; this is known as Onsager reciprocity [28]. This concept is most
naturally expressed in terms of the (dyadic) Green’s tensors that solve
Maxwell’s equation via the inhomogeneous Helmholtz equation. The
correspondence between the sources in terms of Green’s tensors has been
outlined in, for example, [29], thus we are going to state the important
results which form the foundational symmetry argument used in the adjoint
method. The Cartesian component $ij$ of a Green’s tensor describing an
electric field at $\mathbf{r}_{1}$ that originates from a point-like electric
dipole source at $\mathbf{r}_{2}$ can be written as
$\mathcal{G}^{\mathrm{EP}}_{ij}(\mathbf{r}_{1},\mathbf{r}_{2})$. Analogously,
the elements of a Green’s tensor relating magnetic source with its field are
$\mathcal{G}^{\mathrm{HM}}_{ij}(\mathbf{r}_{1},\mathbf{r}_{2})$. We can
express the source-observer symmetry in reciprocal media in terms of
$\mathcal{G}^{\mathrm{EP/HM}}$ as:

${\color[rgb]{0,0,0}\mathcal{G}^{\mathrm{EP/HM}}_{ij}(\mathbf{r}_{2},\mathbf{r}_{1})=\mathcal{G}^{\mathrm{EP/HM}}_{ji}(\mathbf{r}_{1},\mathbf{r}_{2}).}$
(2)

The $i$-th component of an electric/magnetic field at $\mathbf{r}_{2}$ arising
from the $j$-th component of an electric/magnetic dipole at $\mathbf{r}_{1}$
is equal to the $j$-th component of an electric/magnetic field at
$\mathbf{r}_{1}$ from the $i$-th component of an electric/magnetic dipole at
$\mathbf{r}_{2}$. This relation shows that a purely electric or magnetic
source can be exchanged with its respective field at the observation point.
The relationship is not as straightforward for an electric/magnetic field
radiating from a magnetic/electric dipole. If we let
$\mathcal{G}_{ij}^{\mathrm{HP}}$ define a Green’s tensor encapsulating a
magnetic field from an electric dipole, and $\mathcal{G}_{ij}^{\mathrm{EM}}$
contain an electric field from a magnetic dipole, the relation reads (see,
e.g. [29]):

$\displaystyle{\color[rgb]{0,0,0}\mathcal{G}_{ij}^{\mathrm{EM}}(\mathbf{r}_{2},\mathbf{r}_{1})=-\mathcal{G}_{ji}^{\mathrm{HP}}(\mathbf{r}_{1},\mathbf{r}_{2}).}$
(3)

Equations (2) and (3) enable us to simplify the expression describing the
interplay between the fields and our goal structure. This is discussed in
detail in the next section.

###  II.1 Choosing a merit function

To create an algorithm capable of producing a helicity lattice with a desired
geometrical pattern, we need to choose a merit function that will output a
number that encodes the closeness of a given field to the goal, which we will
aim to optimize. Inspired by [30], wherein a method of controlling 3D optical
fields and their intensity is proposed, we seek to control the value of
helicity density at points in 3D space. Since we possess a goal pattern
$\mathcal{H}_{0}$ given by the helicity lattices reported in [8], a natural
choice would be to define a set $S$, of $N$ critical points, where
$\mathcal{H}_{0}$ is “close” to maximum and minimum values respectively. A
simulated helicity density $\mathcal{H}[\mathbf{E},\mathbf{H}]$, can be
compared with the values $\mathcal{H}_{0}$ takes at these coordinates; hence a
measure of similarity between the two fields can be established. Assuming that
our goal structure lies in the plane $z=z_{0}$, we define the set of points
$S$, as a collection of $N$ tuples of the form
$\mathbf{s}_{i}=\\{x_{i},y_{i},z_{0}\\}$. We then seek points $\mathbf{s}_{i}$
for which $\mathcal{H}_{0}(\mathbf{s}_{i})$ is close to its minimum and
maximum. We can then define an arbitrary bound $a\in(0,1]$ and form the set
$S$ of size $N$, such that $\mathbf{s}_{i}\in S$ if

$\displaystyle\mathcal{H}_{0}(\mathbf{s}_{i})<a\cdot\mathrm{\mathrm{min}}[\mathcal{H}_{0}]\text{
or
}\mathcal{H}_{0}(\mathbf{s}_{i})>a\cdot\mathrm{\mathrm{max}}[\mathcal{H}_{0}].$
(4)

Our simulation space is discretised, thus the $N$ is finite and will depend on
the size of the domain and the resolution, as well as the value of $a$. The
merit function $\Theta[\mathbf{E},\mathbf{H}]$ will thus take the form:

$\displaystyle\Theta[\mathbf{E},\mathbf{H}]=\sum_{k=1}^{N}{\color[rgb]{0,0,0}\theta_{k}^{2}(\mathbf{s}_{k}),}$
(5)

where,

$\displaystyle{\color[rgb]{0,0,0}\theta_{k}(\mathbf{s}_{k})\equiv\mathcal{H}[\mathbf{E},\mathbf{H}](\mathbf{s}_{k})-\mathcal{H}_{0}(\mathbf{s}_{k}),}$
(6)

is the error between the goal pattern and the simulated field at a coordinate
$\mathbf{s}_{k}$. The closer $\Theta$ is to zero, the closer
$\mathcal{H}[\mathbf{E},\mathbf{H}]$ will approximate the desired pattern
$\mathcal{H}_{0}$. The merit function $\Theta$ is a sum of $N$ functionals
$\theta_{k}^{2}$, that are evaluated at $\mathbf{s}_{k}$; by linearity, their
variation can be evaluated as

$\displaystyle\delta\Theta[\mathbf{E},\mathbf{H}]$
$\displaystyle=\delta\sum_{k=1}^{N}\theta_{k}^{2}(\mathbf{s}_{k})=2\sum_{k=1}^{N}\theta_{k}(\mathbf{s}_{k})\delta\theta_{k}(\mathbf{s}_{k})$
$\displaystyle=2\sum_{k=1}^{N}\theta_{k}(\mathbf{s}_{k})\delta\mathcal{H}[\mathbf{E},\mathbf{H}](\mathbf{s}_{k}).$
(7)

We want to calculate the $\delta\mathcal{H}[\mathbf{E},\mathbf{H}]$ at each
site in $S$. $\mathcal{H}[\mathbf{E},\mathbf{H}]$ is dependent on the complex
fields $\mathbf{E}$ and $\mathbf{H}$, suggesting that we write

$\displaystyle\delta\mathcal{H}[\mathbf{E},\mathbf{H}]=\,$
$\displaystyle\dfrac{\partial{\mathcal{H}}}{\partial{\mathbf{E}}}\cdot\delta\mathbf{E}+\dfrac{\partial{\mathcal{H}}}{\partial{\mathbf{E}^{*}}}\cdot\delta\mathbf{E}^{*}$
$\displaystyle+\dfrac{\partial{\mathcal{H}}}{\partial{\mathbf{H}}}\cdot\delta\mathbf{H}+\dfrac{\partial{\mathcal{H}}}{\partial{\mathbf{H}^{*}}}\cdot\delta\mathbf{H}^{*}.$
(8)

By plugging in the definition of helicity density from Eq. (1) into Eq. (II.1)
and simplifying, we obtain

$\displaystyle\delta\mathcal{H}[\mathbf{E},\mathbf{H}]=-\dfrac{i}{2c\omega}\mathrm{Re}\left\\{\mathbf{H^{*}}\\!\cdot\delta\mathbf{E}-\mathbf{E^{*}}\\!\cdot\delta\mathbf{H}\right\\},$
(9)

where the all fields are evaluated at $\mathbf{s}_{k}$.

The final step in the method is calculating $\delta\mathbf{E}$ and
$\delta\mathbf{H}$. We employ an approach called shape calculus, which is
established on understanding the relationship between Maxwell’s equations and
small changes in geometry [31]. In our present problem, this manifests as the
link between introducing small pieces of dielectric material and their effect
on the existing electric and magnetic fields. The principle behind evaluating
these can be illustrated with the example of the electric field $\mathbf{E}$.
The effect of modifying the geometry by adding a piece of a polarizable
dielectric will cause a change
$\delta\mathbf{E}=\mathbf{E}_{\mathrm{after}}-\mathbf{E}_{\mathrm{before}}$ in
the electric field, with a similar argument applying for the magnetic field
$\mathbf{H}$;
$\delta\mathbf{H}=\mathbf{H}_{\mathrm{after}}-\mathbf{H}_{\mathrm{before}}$.
As discussed in detail in, for example, [29], the variations in the fields
$\mathbf{E}$ and $\mathbf{H}$ are given by a product of a Green’s tensor
$\mathcal{G}$ with a source current $\mathbf{P}$,

$\displaystyle\delta\mathbf{E}(\mathbf{s}_{k})=$
$\displaystyle\int_{\phi}\mathrm{d}^{3}\mathbf{r}\
\mathcal{G}^{\mathrm{EP}}(\mathbf{s}_{k},\mathbf{r})\mathbf{P}(\mathbf{r})$
(10) $\displaystyle\delta\mathbf{H}(\mathbf{s}_{k})=$
$\displaystyle\int_{\phi}\mathrm{d}^{3}\mathbf{r}\
\mathcal{G}^{\mathrm{HP}}(\mathbf{s}_{k},\mathbf{r})\mathbf{P}(\mathbf{r}),$
(11)

where $\phi$ is a 3D region containing the inclusion. Taking into account the
reciprocity relations in Eq. (3), we can write the variation in
$\theta_{k}^{2}$ at a point $\mathbf{s}_{k}$ as:

$\displaystyle\delta\theta_{k}^{2}(\mathbf{s}_{k})=-\dfrac{i}{2c\omega}\int_{\phi}\mathrm{d}^{3}\mathbf{r}\
\mathrm{Re}\left\\{\mathbf{P}(\mathbf{r})\cdot\mathbf{\mathbf{F}}_{k}(\mathbf{r},\mathbf{s}_{k})\right\\},$
(12)

where we defined;

$\displaystyle\mathbf{\mathbf{F}}_{k}(\mathbf{r},\mathbf{r}^{\prime})=[\mathbf{H}^{*}(\mathbf{s}_{k})$
$\displaystyle\mathcal{G}^{\mathrm{EP}}(\mathbf{r},\mathbf{r}^{\prime})$
$\displaystyle+\mathbf{E^{*}}(\mathbf{r}^{\prime})\mathcal{G}^{\mathrm{EM}}(\mathbf{r},\mathbf{r}^{\prime})]\theta_{k}(\mathbf{r}^{\prime}).$
(13)

Here lies the most important step of the adjoint method; since the structure
of Eq. (II.1) is that of an electric field, that is, a Green’s tensor
multiplying a source current at the observation point $\mathbf{s}_{k}$,
$\mathbf{F}_{k}(\mathbf{r},\mathbf{s}_{k})$ is the $k$th adjoint field.
Exploiting the source-observer symmetry in reciprocal media, each source at
$\mathbf{s}_{k}$ provides us with a field profile at $\mathbf{r}$.

We can further simplify Eq. (12) by taking the dielectric inclusion’s relative
permittivity $\varepsilon_{\mathrm{I}}$ to be homogenous and isotropic. The
induced polarizability, $\mathbf{P}(\mathbf{r})$ can then be written as

$\displaystyle\mathbf{P}(\mathbf{r})=\Delta\varepsilon\mathbf{E}(\mathbf{r}),$
(14)

where $\Delta\varepsilon\equiv\varepsilon_{\mathrm{I}}-\varepsilon^{(0)}$, is
the difference in relative permittivity between the inclusion and the
background $\varepsilon^{(0)}$. Moreover, we assume that the inclusion is
sufficiently small that we can take its centre to be at $\mathbf{r^{\prime}}$,
thus approximating the fields in the integrand of Eq. (12) by their value at
$\mathbf{r^{\prime}}$. This gives

$\displaystyle\delta\theta_{k}^{2}(\mathbf{s}_{k})\approx-\dfrac{i\Delta\varepsilon
V}{2c\omega}\mathrm{Re}\left\\{\mathbf{E}(\mathbf{r^{\prime}})\cdot\mathbf{\mathbf{F}}_{k}(\mathbf{r^{\prime}},\mathbf{s}_{k})\right\\},$
(15)

where $V$ is the volume of the inclusion. We are only interested in relative
changes in the merit function between sites, thus the scalar prefactor
$\left(-i\Delta\varepsilon V/2c\omega\right)$ can be ignored. This allows us
to finally write the variation in $\theta^{2}_{k}$ at a site $\mathbf{s}_{k}$
as

$\delta\theta^{2}_{k}(\mathbf{s}_{k})\approx\mathrm{Re}\left\\{\mathbf{E}(\mathbf{r^{\prime}})\cdot\mathbf{\mathbf{F}}_{k}(\mathbf{r^{\prime}},\mathbf{s}_{k})\right\\},$
(16)

and by linearly adding $N$ adjoint fields we have for the total variation:

$\displaystyle\delta\Theta\left[\mathbf{E},\mathbf{H}\right]$
$\displaystyle\approx\mathrm{Re}\left\\{\mathbf{E}(\mathbf{r^{\prime}})\cdot\sum_{k=1}^{N}\mathbf{\mathbf{F}}_{k}(\mathbf{r^{\prime}},\mathbf{s}_{k})\right\\}$
$\displaystyle=\mathrm{Re}\left\\{\mathbf{E}(\mathbf{r^{\prime}})\cdot\mathbf{F}(\mathbf{r^{\prime}})\right\\}.$
(17)

The change (variation) in the merit function $\Theta$ is thus given by the
overlap between the forward field $\mathbf{E}$ and the adjoint field
$\mathbf{F}$.

###  II.2 Advantages of adjoint method over brute force optimization

Figure 2: A brute-force optimization in a) and adjoint method in b). In a) the
algorithm has to iterate through every coordinate in 3D space placing a test
inclusion there and running a separate simulation before it can determine the
optimal location. In b) the adjoint method finds this coordinate by doing only
two simulations. Both a) and b) continue until the merit function value can no
longer decrease.

To find a structure producing a helicity lattice, one might either use
intuition or follow a brute-force process of building it. The latter involves
trying every coordinate one by one and measuring what effects placing the
inclusion there has on the merit function. One then would narrow it down to
the single coordinate where the dielectric inclusion exerts the largest
influence and place it there. The whole procedure would have to be repeated
for every element added to the structure until the goal field is reached. A 3D
simulation would require a number of iterations equal to the number of
inclusions needed to reach the goal, multiplied by $n^{3}$, where $n$ is the
number of grid points along the side of the simulation box.

By employing the adjoint method instead, we no longer need to check each
coordinate to see how the function responds to adding an inclusion. This
information is contained in the $\delta\Theta$, an array of values
corresponding to each coordinate in the domain. By picking the highest value
of the array, i.e. placing the dielectric at the coordinate with the greatest
number, we are guaranteed to produce the highest change in the function
$\Theta$. This process saves computational resources and makes certain
optimization goals feasible by dramatically decreasing simulation time.

The comparison between the brute force and adjoint approaches is illustrated
schematically in Fig. 2. The brute force approach (left) scans each coordinate
and measures the impact on the merit function by placing a test inclusion
there. Only after finding the location with the highest influence the
algorithm commits to placing the inclusion permanently. The adjoint method
(right) scans the 3D space for the highest value of $\delta\Theta$ from Eq.
(17). It then places an inclusion there, circumventing the need for additional
simulations.

##  III Procedure

Figure 3: Schematic illustration showing the adjoint optimisation cycle (boxes
ABCD), and the depiction of the comparison between the goal structure and the
simulated field in E. The control points are depicted in yellow, both
pictorially, and symbolically, i.e., $\mathbf{s}_{k}$, inside the expressions.
The red dots, along with the red colouring of the fields’ $\mathbf{E}^{*}$ and
$\mathbf{H^{*}}$symbols, signify the sources of the adjoint field (situated at
the yellow control points). Starting the cycle at A, we obtain the values of
the fields $\mathbf{E}(\mathbf{s}_{k})$ and $\mathbf{H}(\mathbf{s}_{k})$ at
the yellow control points, arising due to the single light source. In B, we
use the field values $\mathbf{E}(\mathbf{s}_{k})$ and
$\mathbf{H}(\mathbf{s}_{k})$ from A, and use them as amplitudes of the adjoint
dipoles (in red) that we excite at every control point $\mathbf{s}_{k}$, to
produce the adjoint field $\mathbf{F}$. Then in C, we combine the fields from
A and B obtaining the derivative of the merit function, i.e.,
$\delta\Theta=\mathbf{E}\cdot\mathbf{F}$, and find the coordinate where the
value of $\delta\Theta$ is the highest. In D, we place the dielectric
inclusion at this coordinate, and re-run the simulation in A with the updated
geometry. The cycle repeats until the value of the merit function $\Theta$ is
sufficiently low, i.e., the fields in the box E look alike. The mesh in C and
D shows that a maximum value of $\delta\Theta$ is sought across the 3D space;
the algorithm is not favouring any particular plane, only a south facing wall
is being shown for clarity of the image.

We aim to create a two-dimensional pattern embedded in a three-dimensional
region; it is, therefore, sufficient to generate EM waves that agree with the
goal lattice only at the particular plane where the structure resides (say,
$z=0$). Every point in space apart from the region immediately surrounding the
observation slice is thus a potential target for the placement of a dielectric
inclusion. Since we want our structure to allow for physical access to the
helicity lattice, we leave a margin of width $l$ at each side of the $z=0$
slice. The simulation space is discretized in units of wavelength; any design
that we produce can be adapted to an arbitrary length regime as Maxwell’s
equations are scale invariant (in the absence of charges or currents, as is
the case here). We define our optimisation target as $\mathcal{H}_{0}$, which
is the normalized (between $-1$ and $1$), analytical expression for the chosen
helicity lattice. We will choose the parameter $a$, in line with Eq. (4),
effectively controlling the tolerance in any small mismatch between the
lattice produced and the goal.

The first step of the optimization process is to run the simulation in the
specified region and record the values of the steady-state fields $\mathbf{E}$
and $\mathbf{H}$,  as shown in the box A in Fig. 3. We then construct the
sources for the adjoint simulation (box B in Fig. 3); at each site
$\mathbf{s}_{k}$, we place dipoles whose amplitude depends on the steady-state
values of the fields and the pattern defined at $\mathbf{s}_{k}$, i.e.
$\mathbf{E}(\mathbf{s}_{k}),\mathbf{H}(\mathbf{s}_{k})$ and
$\mathcal{H}_{0}(\mathbf{s}_{k})$. At each $\mathbf{s}_{k}$, we position both
an electric dipole with the amplitude
$\theta_{k}(\mathbf{s}_{k})\mathbf{H^{*}}(\mathbf{s}_{k})$, and a magnetic
dipole with the amplitude
$\theta_{k}(\mathbf{s}_{k})\mathbf{E^{*}}(\mathbf{s}_{k})$, where
$\theta(\mathbf{s}_{k})$ is defined in Eq. (6) and depicted schematically in
box E of Fig. 3. Running the backward (adjoint) simulation concludes the first
step. Upon recording the values of the steady-state fields $\mathbf{F}$, we
combine them with the values of $\mathbf{E}$ via Eq. (17). The result
encapsulates the relationship between the placement of a dielectric inclusion
and a corresponding change to the merit function $\Theta$. In other words,
picking a coordinate $\mathbf{r}_{\mathrm{max}}$ such that
$\mathrm{Re}\left\\{\mathbf{E}(\mathbf{r}_{\mathrm{max}})\cdot\mathbf{F}(\mathbf{r}_{\mathrm{max}})\right\\}$
is maximum, ensures that placing an inclusion at $\mathbf{r}_{\mathrm{max}}$
will result in the largest change in $\Theta$ (boxes C and D in Fig. 3). We
then update the geometry and repeat the previous step until the desired
pattern has been achieved.

The overall algorithm can be summarized in the following way;

  1. 1.

First stage:

    1. (i)

Define a 3D simulation space; choose a plane where the 2D helicity pattern is
to be created and situate its normalized analytical expression there.

    2. (ii)

In the chosen plane, pick a set $S$ of critical points $\mathbf{s}_{k}$, where
the analytical function is “close” to its maximum absolute values (box E in
Fig. 3).

    3. (iii)

By running the forward simulation, obtain steady-state values of the fields
$\mathbf{E}$ and $\mathbf{H}$ at every point in the 3D domain, including the
points of interest $\mathbf{s}_{k}$ (box A in Fig. 3).

    4. (iv)

Place electric and magnetic dipoles at every $\mathbf{s}_{p}$ as functions of
the $\mathbf{E}(\mathbf{s}_{p})$ and $\mathbf{H}(\mathbf{s}_{p})$ at the first
step. Obtain the steady state values of the fields by running the adjoint
simulation (box B in Fig. 3).

    5. (v)

Combine the forward and adjoint fields, as in Eq. (17). Pick the coordinate
where the $\delta\Theta$ is the highest and place the dielectric inclusion
there (boxes C and D in Fig. 3).

  2. 2.

Second stage (iterative):

    1. (i)

Repeat the steps (iii)-(v) from the first stage.

    2. (ii)

Repeat until the forward field approximates the desired pattern at the chosen
plane.

##  IV Results

Using the single beam technique supported by the adjoint method, we have
created the patterns approximating the helicity density originally arising as
a superposition of three waves in two different configurations. While we have
chosen just two, many other non-interfering superpositions are realizable, and
the list of their explicit constituent light sources can be found in [8].

The parameters defining the simulations are as follows. All distances in all
computations are expressed in units of wavelength, and each has a chosen
resolution parameter $R$ such that there are $R$ pixels per wavelength. The
dielectric inclusions are chosen to be cubic with side length equal to $1/R$
wavelengths (i.e., one pixel) and are referred to hereafter as ‘blocks’. In
all simulations the relative permittivity of the blocks was $\varepsilon=1.3$,
the coefficient $a$ in Eq. (4) was set to $0.3$, the gap $l$ was chosen to be
three blocks (i.e. three pixels, or $3/R$ wavelengths), and the incident plane
wave is $y$-polarised and propagating in the $x$ direction. All simulations
were performed using the FDTD library MEEP [32], within which we used a built-
in implementation of perfectly matched layers [33] as boundary conditions at
the edge of the simulation box. The underlying code, along with detailed
documentation can be found online at Ref. [34].

###  IV.1 Rectangular three wave superposition

We first use the above-described method to produce the simplest two-
dimensional helicity lattice, namely a rectangular three-wave superposition
that results in a checkerboard pattern (see Table. 2 and Fig. 3 in [8]), but
here from a single plane wave input. The goal and optimised patterns are shown
in the upper two panels in Fig. 4, with the lower two panels here (and in the
two subsequent figures) included for later discussion in Section V.

Figure 4: Comparison between the goal pattern (upper left panel) and that
obtained by using an inverse-designed structure and a single source (upper
right panel) for an originally three-source helicity lattice in a checkerboard
configuration. The lower panels show the obtained field’s normalised
inhomogeneity in the square of the electric (left) and magnetic (right)
fields, respectively. The simulation was run for 4000 iterations with
resolution of 10 pixels per wavelength $\lambda$.

We observe that since the source plane wave oscillates perpendicular to the
lattice’s vertical lines of symmetry, the resulting simulation approximates
the helicity lattice well (showing clear periodicity by agreeing with the goal
pattern at the $\mathcal{H}=0$ points) even in lower resolution settings (i.e.
less than six pixels per wavelength). The reason for this is suggested by
noting that the goal structure $\mathcal{H}_{0}$, has the form:

$\displaystyle\mathcal{H}_{0}=\mathcal{N}\cos\big{[}(2-\sqrt{2})\pi
x\big{]}\sin(\sqrt{2}\pi y),$ (18)

where the factor $\mathcal{N}$ carries the units of
$\varepsilon_{0}|E|^{2}/2\omega$. This lattice has two orthogonal lines of
symmetry: $x=(2-\sqrt{2})n/2$, and $y=n/\sqrt{2}$, $n\in\mathbb{N}$, where
$\mathcal{H}_{0}$ is zero. The algorithm can thus independently influence the
fields along those degrees of freedom more “naturally”, resulting in the
outline of the correct structure emerging at a very early point of the
simulation ($<100$ iterations, well before the 4000th iteration shown in Fig.
4).

###  IV.2 Triangular (diamond) three-wave superposition

The same number of sources as in the previous example, albeit angled
differently [8] produces a new lattice. The analytical expression in this case
reads:

$\displaystyle\mathcal{H}_{0}=\mathcal{N}\left\\{-\sin(2\pi
x)+\sin[2\pi(x-y)]+\sin(2\pi y)\right\\},$ (19)

The results of the optimisation are plotted in Fig. 5,

Figure 5: Comparison between the goal pattern (upper left panel) and that
obtained by using an inverse-designed structure (see Fig. 7a) and a single
source (upper right panel) for an originally three-wave helicity lattice in a
triangular (diamond) configuration. The lower panels show the obtained field’s
normalised inhomogeneity in the square of the electric (left) and magnetic
(right) fields, respectively. The simulation was run for 2000 iterations with
a resolution of 14 pixels per wavelength $\lambda$.

and for this calculation we also show in Fig. 7a the geometry of the resulting
3D dielectric crystal which, as anticipated, shows no recognisable structures
that could have been reached by intuition.

Figure 6: Comparison between an arbitrary user-defined helicity image (upper
left panel) and that obtained by using an inverse-designed structure (see Fig.
7b) and a single source (upper right panel) producing a helicity density
pattern in the shape of a letter “G” . The lower panels show the obtained
field’s normalised inhomogeneity in the square of the electric (left) and
magnetic (right) fields, respectively. The simulation was run for 16000
iterations with a resolution of 10 pixels per wavelength $\lambda$. Figure 7:
The three-dimensional models of the crystal-like structure. In a), the
arrangement produces the three-wave triangular (diamond) lattice as shown in
Fig. 5. In b), the structure produces an arbitrary pattern (letter “G”) as
shown in Fig. 6. In both cases the gap at the level $z=0$ is three blocks
high.

In contrast to that described by Eq. (18), this lattice possesses more than
two lines of symmetry. This can be inferred from inspecting the arguments of
the $\sin$ function in (19); for $n\in\mathbb{N}$, lines $x=n/2$,$y=n/2$ and
$y=x-n/2$ are the zeros of $\mathcal{H}_{0}$. As a result, it takes more
iterations for the pattern to become apparent as the initial blocks force the
simulated lattice to have the zeros at $y=x-n/2$, exploiting the wave nature
of the source and modulating the helicity density in the diagonal direction.
The other lines of symmetry emerge later as more blocks are added to the
structure, approximating the original lattice to a precision controlled by the
overall resolution (number of pixels per wavelength).

###  IV.3 Arbitrary helicity pattern

Although our technique has shown itself to be effective in constructing
helicity lattices, there is nothing in principle tying it to a regular,
repeating pattern. The algorithm understands any helicity pattern as a series
of pixels, which can be viewed as an input to the algorithm and therefore
chosen arbitrarily. To demonstrate this, we constructed a greyscale image of
the letter “G”, whose grey values were then normalised to lie between $-1$ and
$1$ in the same way as the lattices previously investigated. The algorithm was
tasked with reproducing this image; the results are shown in Fig. 6 and the
corresponding structure is shown in Fig. 7b. We note that the algorithm
correctly distributes the regions of positive and negative helicity density
(the inside of the G is positive), whereas the outline of the letter has
helicity density close to zero. The shape itself is recognisable; however,
some distortion occurs and artifacts are present even though the simulation
was left to run for much longer than the lattices discussed in the previous
two sections (16000 iterations, rather than 4000 and 2000, respectively). The
process constructs the base and the top of the G first as these align with
oscillations of the source field; the “ripples” of positive helicity spread to
the bounds of the shape, and further details are added successively later. The
size of the simulation volume constrains the quality of the resulting lattice;
such images without periodicity and increased level of detail require a large
computational volume. Additionally, if one was mostly concerned with the
fidelity of the reproduced image, another input beam could be incorporated
into the process to improve the quality and speed up the process.

##  V Analysis

In the case of analytical solutions for helicity lattices (four sources or
fewer), the electric field intensity can always be made homogeneous [8]. This
is not always the case for simulated structures, as can be observed in Figs.
4-6. This important detail influences their potential experimental use. As
mentioned in the introduction, forces on atoms and molecules in trapping
potentials is constrained by optical forces induced by gradients in the
electric intensity. Thus, to create an efficient helicity lattice, one would
require the absence of intensity gradients in the region spanned by the
structure to allow for interrupted chiral interactions between the enantiomers
and the lattice. The fields generated by our simulation possess
inhomogeneities due to the finite width of the plane waves generating the
helicity. For the lattices shown in Figs. 4 and 5, the inhomogeneity is
dependent on the size of a computational vs observational region; the larger
the gap between the boundaries of a computational volume and an observation
box, the closer the source resembles a plane wave. This, in turn, produces an
electric field intensity with smaller variations; an infinitely sized
computational box would allow one to obtain a perfectly homogenous intensity.
Consequently, a trade-off exists between a reasonable computational time and
fidelity of the lattice compared to the goal structure. Routes towards
reducing the larger inhomogeneity in intensity exhibited by arbitrary patterns
such as that shown in Fig. 6 are not clear to us. Since they are artificially
introduced as a result of interacting fields, the relationship between their
shapes and the resulting intensity gradients is unclear. Future work will
address this issue and focus on developing more sophisticated methods for
producing an arbitrary pattern. An additional merit function could be defined,
specifically controlling the intensity gradients. Additionally, a weight
function between the two merit functions would be introduced. A link could
potentially be established between certain shapes and the homogeneity of the
intensity; this could be revealed by a machine-learning algorithm.

As can be seen in Fig. 7, some of the elements of the crystal are “floating”.
This is a particularly extreme example of the structural integrity issues that
arise in general in inverse design (see, e.g. [35]). Such problems can be
rectified by running a longer simulation, effectively allowing the algorithm
to fill the larger gaps by chance and then manually removing any remaining
“islands” of the dielectric material. Alternatively, a weight function can be
implemented to penalise the algorithm for placing isolated blocks. A similar
issue occurs as the structure comprises two halves, allowing access to the
lattice; the crystal can be grown within the bounds of an existing structure
providing support, or additional dielectric elements can be incorporated to
join the two halves near the edges. When considering manufacturing challenges
for structures such as this, it is worth recalling that Maxwell’s equations
(without charges or currents) are scale-invariant, and therefore so are the
helicity lattices and the structures proposed here. Thus, at this proof-of-
principle stage, the choice of the size regime will dictate the difficulty of
potential manufacturing challenges.

##  VI Summary and conclusions

In this work, a gradient-based formalism was used to develop a proof-of-
principle scheme capable of reproducing a helicity lattice pattern, only with
a reduced number of input sources. The problem was cast as a minimum-seeking
algorithm, iteratively growing a crystal-like structure. This object refracts
a single beam in such a way as to reproduce the effect of superposing multiple
beams, resulting in the formation of a helicity lattice. We have presented the
results of reproducing two originally three-wave helicity lattices using a
single source input. Additionally, we introduced a modified scheme where the
pattern to be reproduced was not a result of superposing light sources but
rather an arbitrary image provided by the user. This approach opens up a
possibility for designing bespoke helicity patterns suited to the needs of a
particular experiment; it is worth mentioning that using the same principles,
a different merit function can be provided, thus it appears to be possible to
use versions of this method to achieve structures arising from various mutual
arrangements of the electric and magnetic fields, not just those where a
helicity structure arises.

The simulations reproducing known helicity patterns serve as proof-of-
principle concepts and demonstrate the feasibility of using a single light
source and a refracting object to mimic the behaviour of a field composed of
multiple sources. We have described how increasing the computational volume
approximates the plane wave more closely and, as a result, produces an
increasingly homogeneous intensity in the case of known helicity patterns. We
outlined possible strategies to resolve the issue of intensity gradients in
arbitrary patterns bringing practical helicity lattices a step closer to
reality.

###### Acknowledgements.

The authors acknowledge financial support from the UK Research and Innovation
Council through a doctoral training programme grant (EPSRC/DTP
2020/21/EP/T517896/1), New Horizons grant EP/V048449/1 and New Investigator
grant EP/W016486/1. The authors also acknowledge fruitful discussions with
Stephen M. Barnett, Jörg B. Götte, Ben W. Butler, Thomas M. Jones and Zhujun
Ye. We also want to thank Claire Cisowski for testing the early version of the
code, Owen Miller for enlightening correspondence and Neel MacKinnon for
critical reading of the manuscript.

## References

  * Barron [2009] L. D. Barron,  _Molecular light scattering and optical activity_ (Cambridge University Press, 2009).
  * Tang and Cohen [2010] Y. Tang and A. E. Cohen, Optical chirality and its interaction with matter, Phys. Rev. Lett. 104, 163901 (2010).
  * Lipkin [1964] D. M. Lipkin, Existence of a new conservation law in electromagnetic theory, Journal of Mathematical Physics (New York) (U.S.) Vol: 5, 10.1063/1.1704165 (1964).
  * Mackinnon [2019] N. Mackinnon, On the differences between helicity and chirality, Journal of Optics 21, 125402 (2019).
  * Trueba and Rañada [1996] J. L. Trueba and A. F. Rañada, The electromagnetic helicity, European Journal of Physics 17, 141 (1996).
  * Woltjer [1958] L. Woltjer, A theorem on force-free magnetic fields, Proceedings of the National Academy of Sciences 44, 489 (1958), https://www.pnas.org/doi/pdf/10.1073/pnas.44.6.489 .
  * Cameron  _et al._ [2012] R. P. Cameron, S. M. Barnett, and A. M. Yao, Optical helicity, optical spin and related quantities in electromagnetic theory, New Journal of Physics 14, 053050 (2012).
  * Kruining  _et al._ [2018] K. Kruining, R. Cameron, and J. Götte, Superpositions of up to six plane waves without electric-field interference, Optica 5, 1091 (2018).
  * Cameron  _et al._ [2017] R. P. Cameron, J. B. Götte, S. M. Barnett, and A. M. Yao, Chirality and the angular momentum of light, Philos Trans A Math Phys Eng Sci 375, 20150433 (2017).
  * Cheng  _et al._ [2022] Q. Cheng, Q. Ma, H. Pei, and Z. Mo, Chiral membranes for enantiomer separation: A comprehensive review, Separation and Purification Technology 292, 121034 (2022).
  * Isaule  _et al._ [2022] F. Isaule, R. Bennett, and J. B. Götte, Quantum phases of bosonic chiral molecules in helicity lattices, Physical Review A 106, 10.1103/physreva.106.013321 (2022).
  * Nshii  _et al._ [2013] C. C. Nshii, M. Vangeleyn, J. P. Cotter, P. F. Griffin, E. A. Hinds, C. N. Ironside, P. See, A. G. Sinclair, E. Riis, and A. S. Arnold, A surface-patterned chip as a strong source of ultracold atoms for quantum technologies, Nature Nanotechnology 8, 321 (2013), 23563845 .
  * Barker  _et al._ [2019] D. Barker, E. Norrgard, N. Klimov, J. Fedchak, J. Scherschligt, and S. Eckel, Single-Beam Zeeman Slower and Magneto-Optical Trap Using a Nanofabricated Grating, Physical Review Applied 11, 064023 (2019).
  * Zhu  _et al._ [2020] L. Zhu, X. Liu, B. Sain, M. Wang, C. Schlickriede, Y. Tang, J. Deng, K. Li, J. Yang, M. Holynski, S. Zhang, T. Zentgraf, K. Bongs, Y.-H. Lien, and G. Li, A dielectric metasurface optical chip for the generation of cold atoms, Science Advances 6, eabb6667 (2020).
  * Dobson and Cox [1999] D. C. Dobson and S. J. Cox, Maximizing band gaps in two-dimensional photonic crystals, SIAM Journal on Applied Mathematics 59, 2108 (1999).
  * Felici and Engl [2001] T. Felici and H. W. Engl, On shape optimization of optical waveguides using inverse problem techniques, Inverse problems 17, 1141 (2001).
  * Volpe and Melnik [1986] G. Volpe and R. E. Melnik, The design of transonic aerofoils by a well‐posed inverse method, International Journal for Numerical Methods in Engineering 22, 341 (1986).
  * Borel  _et al._ [2004] P. I. Borel, A. Harpøth, L. H. Frandsen, M. Kristensen, P. Shi, J. S. Jensen, and O. Sigmund, Topology optimization and fabrication of photonic crystal structures, Optics express 12, 1996 (2004).
  * Jensen and Sigmund [2004] J. S. Jensen and O. Sigmund, Systematic design of photonic crystal structures using topology optimization: Low-loss waveguide bends, Applied physics letters 84, 2022 (2004).
  * Molesky  _et al._ [2018] S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, Inverse design in nanophotonics, Nature Photonics 12, 659 (2018).
  * Xiao  _et al._ [2016] T. P. Xiao, O. S. Cifci, S. Bhargava, H. Chen, T. Gissibl, W. Zhou, H. Giessen, K. C. Toussaint Jr, E. Yablonovitch, and P. V. Braun, Diffractive spectral-splitting optical element designed by adjoint-based electromagnetic optimization and fabricated by femtosecond 3d direct laser writing, ACS Photonics 3, 886 (2016).
  * Piggott  _et al._ [2015] A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vučković, Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer, Nature Photonics 9, 374 (2015).
  * Dobson and Simeonova [2009] D. C. Dobson and L. B. Simeonova, Optimization of periodic composite structures for sub-wavelength focusing, Applied Mathematics and Optimization 60, 133 (2009).
  * Jin  _et al._ [2017] W. Jin, R. Messina, and A. W. Rodriguez, Overcoming limits to near-field radiative heat transfer in uniform planar media through multilayer optimization, Optics Express 25, 14746 (2017).
  * Bennett [2021] R. Bennett, Inverse design of environment-induced coherence, Phys. Rev. A 103, 013706 (2021).
  * Matuszak  _et al._ [2022] J. Matuszak, S. Y. Buhmann, and R. Bennett, Shape optimizations for body-assisted light-matter interactions, Phys. Rev. A 106, 013515 (2022).
  * Diekmann  _et al._ [2022] O. Diekmann, D. Lentrodt, and J. Evers, Inverse design approach to x-ray quantum optics with mössbauer nuclei in thin-film cavities, Phys. Rev. A 105, 013715 (2022).
  * Onsager [1931] L. Onsager, Reciprocal Relations in Irreversible Processes. I., Physical Review 37, 405 (1931).
  * Miller [2012] O. Miller,  _Photonic Design: From Fundamental Solar Cell Physics to Computational Inverse Design_ , Ph.D. thesis, EECS Department, University of California, Berkeley (2012).
  * Zhan  _et al._ [2019] A. Zhan, R. Gibson, J. Whitehead, E. Smith, J. R. Hendrickson, and A. Majumdar, Controlling three-dimensional optical fields via inverse mie scattering, Science Advances 5, eaax4769 (2019), https://www.science.org/doi/pdf/10.1126/sciadv.aax4769 .
  * Costabel and Louër [2011] M. Costabel and F. L. Louër, Shape derivatives of boundary integral operators in electromagnetic scattering. part ii: Application to scattering by a homogeneous dielectric obstacle, Integral Equations and Operator Theory 73, 17 (2011).
  * Oskooi  _et al._ [2010] A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. Joannopoulos, and S. G. Johnson, Meep: A flexible free-software package for electromagnetic simulations by the fdtd method, Computer Physics Communications 181, 687 (2010).
  * Berenger [1994] J. P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, Journal of Computational Physics 114, 185 (1994).
  * Kilianski [2022] R. Kilianski, Inverse design of arbitrary optical helicity structures, https://researchdata.gla.ac.uk/1374/ (2022).
  * Augenstein and Rockstuhl [2020] Y. Augenstein and C. Rockstuhl, Inverse design of nanophotonic devices with structural integrity, ACS Photonics 7, 2190 (2020), https://doi.org/10.1021/acsphotonics.0c00699 .


# The Effect of Strange Quark on the Chiral Symmetry Breaking in Magnetic
Background in the framework of logarithmic quark-sigma model

M. Abu-Shady  Department of Applied Mathematics, Faculty of Science, Menoufia
University, Egypt

###### Abstract

The chiral symmetry breaking in the presence of an external magnetic field is
studied in the framework of logarithmic quark-sigma model with three flavors
$\left(u,d,s\right)$. The effective logarithmic mesonic potential is extended
to strange quark and is numerically solved in the mean-field approximation.
The present results show that the chiral symmetry breaking is sensitive for
the strange quark flavour in the presence an external magnetic field. In
addition, the effect of the free parameters is studied. A comparison with
original sigma model is discussed.

Keywords: Chiral Lagrangian density, Magnetic catalysis, Mean-field
approximation

Chiral Lagrangian density, Magnetic catalysis, Chiral symmetry breaking.

###### pacs:

PACS number

††preprint:

LABEL:FirstPage1 LABEL:LastPage#1

##  I Introduction

The quantum chromodyanmic (QCD) theory is an acceptable theory for strong
interactions. A direct technique to calculate the hadron properties through
using the lattice QCD but this is not easy task, in particular, at finite
density. So the investigations are preformed by using effective models such as
the Nambu-Jona-Lasinio model and quark sigma model that share the QCD theory
in same properties. The linear sigma model is introduced by Gell-Man and Levy
$\left[1\right]$ to describe pion-nucleon interactions. This model is extended
to quark level by Birse and Banerjee $\left[2\right]$ to calculate hadron
properties at low energy. Birse and Banerjee model $\left[2\right]$ has some
difficulties that conflict with experimental data. So, the model is extended
to avoid these difficulties as in Refs. $\left[3-5\right]$. In addition, the
model is extended to include finite temperature as Refs. $\left[6-8\right]$.

The study of the influence of external magnetic fields on the fundamental
properties of quantum chromodynamic (QCD) theory, confinement and dynamical
chiral symmetry breaking is still a matter of great interest theoretical and
experimental activities. In Ref. $\left[9\right]$, the fact that an external
magnetic field enhances the generation of a fermion mass in 3+1 dimensions was
first established in the framework of the Nambu-Jona-Lasinio (NJL) model made
in Refs. $\left[10,11\right]$. Also, in the context of the (2+1)-dimensional
Gross-Neveu model, which was expected to give a simple effective description
for certain condensed matter planar systems, the authors of Refs.
$\left[12,13\right]$ showed that the dynamical generation of nonzero fermion
mass takes place in a magnetic field as soon a there an attractive interaction
between fermions and antifermions. The symmetric quark matter and quark matter
in $\beta$ equilibrium are investigated in the NJL model $\left[14\right]$ in
the presence external magnetic field at zero temperature $\left[15,16\right]$.
The effective quark model that takes into account chiral symmetry
$\left[17,18\right]$. It was shown that non-negligible effects on the equation
of state and single-particle quark properties. The linear sigma model with two
flavors has also applied to determine the critical point temperature in the
presence of external magnetic field $\left[19,20\right]$. The effect of the
higher-order mesonic interactions on the chiral symmetry breaking is
investigated in the presence of an external magnetic field $\left[21\right]$.

In Ref. $\left[22\right],$ the authors have modified the linear sigma model by
including the logarithmic mesonic potential and study its effect on the phase
transition at finite temperature. In addition, the comparison with other
models is done. On the same hand, the logarithmic sigma model successfully to
describe the nucleon properties at finite temperature and density
$\left[23\right]$. In addition, the logarithmic sigma model give good
descriptions of magnetic catalysis for two flavors $\left[24\right]$.

To continue the investigation that started in Ref. $\left[22\right]$, the
effect of strange quark field on the chiral symmetry breaking in the framework
of logarithmic sigma model is investigated in the presence of an external
magnetic field. So for no attempts have done to investigate the logarithmic
sigma model for three flavors $\left(u,d,s\right)$ in the presence of an
external magnetic field. In addition, a new parametrization set for sigma mass
and coupling constant on the behavior of chiral symmetry breaking is studied.

This paper is organized as follows: The background of previous papers is given
in Sec. 1. The original sigma model is briefly presented in Sec. 2. Next, the
effective logarithmic mesonic potential in the presence of external magnetic
field is presented in Sec. 3. The results are discussed and are compared with
other models in Secs. 4 and 5, respectively. Finally, the summary and
conclusion are presented in Sec. 6.

##  II The chiral sigma model with original effective potential

The interactions of quarks via the exchange of $\sigma-$ and $\mathbf{\pi}$
-meson fields are given by the Lagrangian density $[2]$ as follows:

${\small
L}\left(r\right){\small=i}\overline{\Psi}{\small\partial}_{\mu}{\small\gamma}^{\mu}{\small\Psi+}\frac{1}{2}\left(\partial_{\mu}\sigma\partial^{\mu}\sigma+\partial_{\mu}\mathbf{\pi}.\partial^{\mu}\mathbf{\pi}\right){\small+g}\overline{\Psi}\left(\sigma+i\gamma_{5}\mathbf{\tau}.\mathbf{\pi}\right){\small\Psi-U}_{1}{\small(\sigma,\pi),}$
(1)

with

${\small
U}_{1}{\small(\sigma,\pi)=}\frac{\lambda^{2}}{4}\left(\sigma^{2}+\mathbf{\pi}^{2}-\nu^{2}\right)^{2}+{\small
m}_{\pi}^{2}{\small f}_{\pi}{\small\sigma.}$ (2)

$U_{1}\left(\sigma,\mathbf{\pi}\right)$ is the meson-meson interaction
potential where $\Psi,\sigma$ and $\mathbf{\pi}$ are the quark, sigma, and
pion fields, respectively. In the mean-field approximation, the meson fields
are treated as time-independent classical fields. This means that we replace
the power and products of the meson fields by corresponding powers and the
products of their expectation values. The meson-meson interactions in Eq. (2)
lead to hidden chiral $SU(2)\times SU(2)$ symmetry with $\sigma\left(r\right)$
taking on a vacuum expectation value

$\ \ \ \ {\small\ \
}\left\langle\sigma\right\rangle{\small=-f}_{\pi}{\small,}$ (3)

where $f_{\pi}=93$ MeV is the pion decay constant. In Eq. (2), the final term
is included to break the chiral symmetry explicitly. It leads to the partial
conservation of axial-vector isospin current (PCAC). The parameters
$\lambda^{2}$and $\nu^{2}$ can be expressed in terms of$\ f_{\pi}$, sigma and
pion masses as,

${\small\lambda}^{2}{\small=}\frac{m_{\sigma}^{2}-m_{\pi}^{2}}{2f_{\pi}^{2}}{\small,}$
(4)
${\small\nu}^{2}{\small=f}_{\pi}^{2}{\small-}\frac{m_{\pi}^{2}}{\lambda^{2}}{\small.}$
(5)

##  III The effective logarithmic potential in the presence of magnetic field

In this section, the logarithmic mesonic potential $U_{2}(\sigma,\pi)$ is
employed. In Eq. (6), The effective logarithmic potential is extended to
include the external magnetic field at zero temperature and chemical potential
as follows,

$U_{eff}(\sigma,\pi)=U_{2}(\sigma,\pi)+U_{Vaccum}+U_{Matter},$ (6)

where

$\ \ \ \ \ \
U_{2}\left(\sigma,\mathbf{\pi}\right)=-\lambda_{1}^{2}\left(\sigma^{2}+\mathbf{\pi}^{2}\right)+\lambda_{2}^{2}\left(\sigma^{2}+\mathbf{\pi}^{2}\right)^{2}\log\left(\frac{\sigma^{2}+\mathbf{\pi}^{2}}{f_{\pi}^{2}}\right)+m_{\pi}^{2}f_{\pi}\sigma,$
(7)

In Eq. 7, the logarithmic potential satisfies the chiral symmetry when
$m_{\pi}\longrightarrow 0$ as well as in the standard potential in Eq. 2.
Spontaneous chiral symmetry breaking gives a nonzero vacuum expectation for
$\sigma$ and the explicit chiral symmetry breaking term in Eq. 6 gives the
pion its mass.

$\left\langle\sigma\right\rangle=-f_{\pi}.$ (8)

Where

$\ \ \ \ \ \lambda_{1}^{2}=\frac{m_{\sigma}^{2}-7m_{\pi}^{2}}{12},$ (9) $\ \ \
\ \ \lambda_{2}^{2}=\frac{m_{\sigma}^{2}-m_{\pi}^{2}}{12f_{\pi}^{2}}.$ (10)

For details, see Refs. $\left[22,23\right]$. To include the external magnetic
field in the present model, we follow Ref. $\left[20\right]$ by including the
pure fermionic vacuum contribution in the free potential energy. Since this
model is renormalizable the usual procedure is to regularize divergent
integrals using dimensional regularization and to subtract the ultra violet
divergences. This procedure gives the following result

${\small
U}_{Vaccum}=\frac{N_{c}N_{f}~{}g^{4}}{(2\pi)^{2}}(\sigma^{2}+\pi^{2})^{2}(\frac{3}{2}-\ln(\frac{g^{2}(\sigma^{2}+\pi^{2})}{\Lambda^{2}})),$
(11)

where $N_{c}=3$ and $N_{f}=3$ are color and flavor degrees of freedom.
respectively, and $\Lambda$ is mass scale,

$U_{Matter}=\frac{N_{c}}{2\pi^{2}}\sum_{f~{}=u,d,s}(\left|q_{f}\right|B)^{2}[\zeta^{(1,0)}(-1,x_{f})-\frac{1}{2}(x_{f}^{2}-x_{f})\ln
x_{f}+\frac{x_{f}^{2}}{4}]$ (12)

In Eq. 12, we have used
$x_{f}=\frac{g^{2}(\sigma^{2}+\pi^{2})}{(2\left|q_{f}\right|B)}$ and
$\zeta^{(1,0)}(-1,x_{f})=\frac{d~{}\zeta(z,x_{f})}{dz}\mid_{z~{}=-1}$ that
represents the Riemann-Hurwitz function, and also $\left|q_{f}\right|$ is the
absolute value of quark electric charge in external magnetic field with
intense $B$ and the potential in Eq. 12 is extended to include strange quark
based on Ref. $\left[25\right]$. In Eq. 6, the of effect of the finite
temperature and chemical potential is not included in the present model, in
which the present model focus on the study of magnetic catalysis at low energy
due to the enhancement in the chiral symmetry breaking in the bound state.

##  IV  Discussion of results

In this section, we study the effective potential of logarithmic sigma model.
For this purpose, we numerically calculate the effective potential in Eq. (6).
The parameters of the present model are the coupling constant $\left(g\right)$
and the sigma mass $\left(m_{\sigma}\right)$. The choice of free parameters of
$\left(g\right)$ and $\left(m_{\sigma}\right)$ based on Ref.
$\left[20\right].$ The parameters are usually chosen so that the chiral
symmetry is spontaneously broken in the vacuum and the expectation values of
the meson fields. In this work, we consider two different sets of parameters
in order to get high and a low value for sigma mass. The first set is given by
$\Lambda=16.48$ MeV which yields $m_{\pi}=138$ MeV and $m_{\sigma}=600$ MeV.
The second set as the first, yielding $m_{\sigma}=400$ MeV.

In Fig. 1, the effective potential is plotted as a function of sigma field for
different values of magnetic field $(B)$. By ignoring the $B$ independent one
loop term in Eq. (6). At zero magnetic field, we note that qualitative
agreement between the effect of light quark field and strange quark field. The
difference appears that the effective potential shifts to lower values at
lower values of sigma field. Therefore, we deduce that the spontaneous
symmetry breaking remains unchanged under at zero magnetic field by including
strange quark field. By increasing the external magnetic field as in Fig. (2),
we note that the spontaneous symmetry breaking is clearly appeared and the
potential has largest two minima values in comparison with their values in the
case of light quark. Therefore, the strange quark increases the generated
fermionic mass which leads to increase magnetic catalysis in the present
model.

Fig. 1: The effective potential is plotted as a function of sigma field for
$m_{\sigma}=600$ MeV and $g=4.5$ at zero magnetic field in the chiral limit

Fig. 2: The effective potential is plotted as a function of sigma field for
$m_{\sigma}=600$ MeV and $g=4.5$ strong magnetic field (eB = 0.214
GeV${}^{2})$ in the chiral limit

It is important to discuss the effect of free parameters of the model on the
chiral symmetry breaking. So, we select two sets of parameters. First set, the
change of coupling constant with fixed sigma mass as in Fig. 3. The effective
potential is plotted as a function of sigma field for different values of $g$
at dense of magnetic field $\left(eB=0.214GeV^{2}\right)$. The effect of
coupling constant $g$ strongly clarify when sigma field increases. We note
that effective potential shifts to higher values by decreasing $g$ which means
that the energy of the potential increases with decreasing coupling constant
$\left(g\right)$. In addition, increasing coupling constant $\left(g\right)$
enhances the chiral symmetry breaking and the potential has two largest minima
values in comparison with their values at $g=4.5$. Therefore, the generate
fermionic mass increases. In the second set of parameters, we fixed the
coupling constant $g$ with two values of sigma mass as in Fig. 4.

Fig. 3: The effective potential is plotted as a function of sigma field at
m${}_{\sigma}=600$ MeV and $eB=0.214$ GeV2 for two values of coupling constant
g.

We note that qualitative features of effective potential are remains
unchanged, in which the effective potential is not sensitive up to $\simeq\pm
175$ MeV and then the potential shifts to higher values by increasing sigma
mass. Also, we note that two minima values of the potential increases with
decreasing sigma mass. Therefore, the decrease in the sigma mass enhances the
chiral symmetry breaking at dense of magnetic field.

Fig. 4: The effective potential is plotted as a function of sigma field at
coupling constant g = 4.5 and $eB=0.214$ GeV2 for two values of sigma mass.

Fig. 5: A comparison between the original sigma model and the logarithmic
sigma model at fixed parameters at $4.5$ and ${\small e}B=0.214$ GeV2 and
$m_{\sigma}=600$ MeV .

In Fig. 5, we compare between the original sigma model and the logarithmic
sigma model. We fixed all parameters in the two models to show realistic
effect on the chiral symmetry breaking. Qualitative features are similar for
the two models. In addition, we note that the two minima values of logarithmic
potential decreases in comparison with their values in the original sigma
model.

##  V Summary and conclusion

In this work, the effect strange quark flavor on the spontaneous chiral
symmetry breaking in the presence of magnetic field is studied. In addition,
the effect of free parameters of the model is studied in the presence of
strong magnetic field when strange quark flavor is included. A comparison with
original sigma model is presented, showing qualitative agreement with original
sigma model. So, novelty in this work that the effect of strange quark flavor
on the spontaneous chiral symmetry breaking is studied in the framework of the
logarithmic sigma model. Most previous studies have been concentrated on the
original sigma model and NJL model at zero or finite temperature.

We hope to extend the present model to finite temperature and chemical
potential which play an important role for studying properties of the universe
and neutron star.

##  VI References

  1. 1

M. Gell-Mann and M. Levy, Nuono Cinmento 16, 705 (1960).

  2. 2

M. Birse and M. Banerjee, Phys. Rev. D 31, 118 (1985).

  3. 3

M. Abu-Shady, Int. J. Mod. Phys. A 26, 235 (2011).

  4. 4

M. Abu-Shady, Int. J. Theor. Phys. 48, 1110 (2009).

  5. 5

T. S. T. Aly, M. Rashdan, and M. Abu-Shady, Inter. J. Theor. Phys. 45, 1654
(2006).

  6. 6

M. Abu-Shady, Inter. J. Mod. Phys. E 21, 1250061 (2012).

  7. 7

M. Abu-Shady and M. Soleiman, Phys. of Part and Nucl. Lett 10, 683 (2013)

  8. 8

M. Abu-Shady, Inter. J. of Theor. Phys. 48, 115 (2009)

  9. 9

V. A. Miransky and I. A. Shovkovy, hep-ph$\backslash$150300732 (2015).

  10. 10

S. P. Klevansky, R. H. Lemmer, Phys. Rev. D 39 3478 (1989).

  11. 11

H. Suganuma, T. Tatsumi, Annals Phys. 208, 470 (1991)

  12. 12

K. G. Klimenko, Z. Phys. C 54, 323 (1992).

  13. 13

K. G. Klimenko, Theor. Math. Phys. 90, 1(1992).

  14. 14

Y. Nambu and G. Jona-Lasinio, Phys. Rev. 122, 345, (1961); 124, 246 (1961)

  15. 15

D. P. Menezes, M. B. Pinto, S. S. Avancini, A. P. Martinez, and C.
Providencia, Phys. Rev. C 79, 035807 (2009).

  16. 16

D. P. Menezes, M. B. Pinto, S. S. Avancini, and C. Providencia, Phys. Rev. C
80, 065805 (2009).

  17. 17

D. Ebert and K. G. Klimenko, Nucl. Phys. A 728, 203 (2003).

  18. 18

B. Hiller, A. A. Osipov, A. H. Blin, and J. da Providencia, SIGMA 4, 024
(2008).

  19. 19

J. O. Andersen and R. Khan, Phys. Rev. D 85, 065026 (2012).

  20. 20

G. N. Ferrari, A. F. Garcia, and M. B. Pinto, Phys. Rev. D 86, 096005 (2012).

  21. 21

M. Abu-Shady, J. Part. Phys. 1 (1), 58 (2017) hep-ph$\backslash$1612.02151.

  22. 22

M. Abu-Shady, Mod. Phys. Lett. A 29, 1450176 (2014).

  23. 23

M. Abu-shady and A. K. Abu-Nab, Euro. Phys. J. Plus 130, 248 (2015).

  24. 24

M. Abu-Shady, Appl. Math. Inf. Sci. Lett. 4, 5 (2016).

  25. 25

M. B. Pinto, D. P. Menezes, and C. Providencia, Phys. Rev. C 83, 065805
(2011).


[affiliation=1]BohanLi [affiliation=1]FeiyuShen [affiliation=1]YiweiGuo
[affiliation=2]ShuaiWang [affiliation=1]XieChen [affiliation=1, ^∗]KaiYu

# On the Effectiveness of Acoustic BPE in Decoder-Only TTS

###### Abstract

Discretizing speech into tokens and generating them by a decoder-only model
have been a promising direction for text-to-speech (TTS) and spoken language
modeling (SLM). To shorten the sequence length of speech tokens, acoustic
byte-pair encoding (BPE) has emerged in SLM that treats speech tokens from
self-supervised semantic representations as characters to further compress the
token sequence. But the gain in TTS has not been fully investigated, and the
proper choice of acoustic BPE remains unclear. In this work, we conduct a
comprehensive study on various settings of acoustic BPE to explore its
effectiveness in decoder-only TTS models with semantic speech tokens.
Experiments on LibriTTS verify that acoustic BPE uniformly increases the
intelligibility and diversity of synthesized speech, while showing different
features across BPE settings. Hence, acoustic BPE is a favorable tool for
decoder-only TTS.

###### keywords:

discrete speech token, acoustic byte-pair encoding, decoder-only text-to-
speech

##  1 Introduction

††footnotetext: $\ast$: The corresponding author.

The language modeling paradigm has been revolutionizing text-to-speech (TTS)
by its strong generation ability since the birth of large decoder-only
language models (LMs). There have been LM-based TTS models that exhibits
versatile, expressive and even emergent abilities [1, 2, 3, 4, 5, 6].

Unlike text, speech is intrinsically a continuous signal with a much lower
information density. To adapt speech into LMs and generate it in an
autoregressive manner, researchers have adopted different methods to map
speech into discrete tokens [2, 4]. By the purpose of discretization, discrete
speech tokens can be divided into acoustic tokens [7, 8, 9] that aim to
reconstruct the signal perfectly, and semantic tokens [10, 11, 12, 13], from
self-supervised models that provide a compact abstraction of speech semantics.
After discretization, decoder-only LMs like VALL-E [3] can be seamlessly
applied to TTS, where the discrete speech tokens are treated as the targets
given a text input.

However, the low rate of information behind speech still leads to excessively
lengthy discrete speech sequences compared to text transcriptions [14]. For
example, more than 500 HuBERT [11] tokens may be required for a single
sentence of around 30 words just to convey an idea or two. It is even worse
for the acoustic tokens, since most of neural speech codec models require a
short frame shift and the residual vector quantization [8] technique to
reconstruct the signal decently. Regardless of the type of tokens, this
feature of speech poses a great challenge for long-context modeling of LM-
based TTS systems.

To address this issue, one possible way is to further compress the discrete
speech sequence. A promising approach is the acoustic byte-pair encoding (BPE)
technique, which is proposed in [15]. It is a similar method to the
traditional BPE algorithm [16] in natural language processing. It treats the
discrete indexes of speech as literal characters and iteratively compresses
consecutive tokens based on the frequency in the training corpus. Such
compression will coherently reduce sequence length with the increase of
vocabulary size. For speech discrete tokens, usually a group of multiple
tokens occur together to represent a specific phoneme or syllable, and
organizing them to be a unique modeling unit would provide a higher
abstraction of semantics and morphological information. Therefore, it is
intuitively reasonable to apply acoustic BPE on discrete speech tokens in both
reducing sequence length and improving representability. In previous
researches, such acoustic BPE has been adopted encode pseudo-target label in
HuBERT pretraining [15] and automatic speech recognition [17].

Nevertheless, the effectiveness of applying acoustic BPE in the TTS task still
remains unclear to the literature. Although BASE-TTS [6] and VoxtLM [18]
mentions the use of acoustic BPE in generation, the design space of acoustic
BPE is still not fully investigated, and the gain of such technique in TTS
needs to be further explored. Despite of the possible improvements, acoustic
BPE could bring more difficulties in choosing the correct unit for generation,
since the vocabulary of the LM in TTS is greatly enlarged. Too much
abstraction of acoustic units could also make language modeling harder.
Moreover, acoustic BPE might have different behaviors and performances on
different types of discrete speech tokens.

Therefore, in this work, we conduct a comprehensive study on the effectiveness
of introducing acoustic BPE to TTS in the decoder-only LM paradigm. We
implement various settings of acoustic BPE on the semantic tokens extracted by
speech self-supervised models, and then train a VALL-E autoregressive decoder-
only transformer as the acoustic model to generate acoustic BPEs from text
inputs. Afterwards, the original semantic tokens are unfolded and fed to a
unit-based vocoder [19, 20] for waveform synthesis. We consider HuBERT [11]
and WavLM [13] as the source of semantic tokens, adjust the number of clusters
in extracting the semantic tokens from 2048 to 8192, and increase the
vocabulary size of acoustic BPE up to 20k in order to observe the effects made
by different acoustic BPE settings. We perform most of the experiments on
LibriTTS [21] and evaluate the effectiveness of such settings via speech
intelligibility, sample diversity and speech quality in objective and
subjective measurements.

Our findings suggest that models employing acoustic BPE method are capable of
generating audio with high intelligibility and quality, while also achieving
much faster training and inference processes. This indicates that acoustic BPE
can be of great value for the TTS task when properly configured. This could
lay a solid foundation for future decoder-only LM-based TTS models, as the
pros of acoustic BPE weighs much more than its con.

##  2 Decoder-Only TTS with Acoustic BPE

###  2.1 Acoustic BPE

BPE is a widely-used algorithm for data compression and text encoding. In
natural language processing, it is initialized with a vocabulary that contains
all the characters with their emergence frequency in the text corpus, and
iteratively merges the most frequent character pairs until a specified number
of merges or a target vocabulary size is reached [22]. For most languages like
English, there exists blanks as the obvious boundaries between words in the
text. However, the textual sequences of audio discrete tokens lack obvious
boundaries, akin to linguistic features found in Chinese. Following [14], we
first obtain speech discrete tokens from k-means clusters of features derived
from speech self-supervised learning (SSL) models (HuBERT, WavLM). Then we
establish a bijective mapping for conversion between the speech discrete
tokens and the Chinese-character Unicodes, simply using an integer offset.
Chinese characters mapping to speech discrete tokens are used to train the BPE
model. With these operations, the acoustic BPE encoding is composed of
conversion to Chinese characters and the encoding of pretrained BPE model. The
acoustic BPE decoding is exactly the inverse process. This research utilizes
the acoustic BPE method as a part of the tokenizer, which plays a role in
preprocessing inputs for the decoder-only language model mentioned in the next
subsection(2.2). Similar to natural language processing, we can consider the
upcoming modeling process as a spoken language modeling process.

Based on the method description above, it is evident that the method involves
multiple dimensions of configuration, such as the SSL model, k-means cluster
number, BPE vocabulary size, and so on. Configuring combinations across
different dimensions may lead to entirely new characteristics, which can have
significant implications in spoken language modeling. A comprehensive and
systematic exploration of these configurations' impact on modeling is crucial
for better constructing spoken language models.

###  2.2  Decoder-Only TTS Language Modeling

The language model can estimate the probability of the input sequence, which
is usually expressed as:

$p(\mathbf{x}|\theta)=\prod_{i=1}^{t}p(x_{i}|x_{<i},\theta)$ (1)

where $\mathbf{x}=\\{x_{1},...,x_{t}\\}$ is the input sequence with length $t$
and the language model is parameterized by $\theta$.

However, the spoken language model in TTS pays the other attention on text
conditions. Based on the decoder-only Transformer architecture, the training
and inference process of the TTS model follows the autoregressive (AR) stage
in VALL-E.

####  2.2.1 Decoder-only TTS model

During training process, we put the phoneme sequence $\mathbf{x}$ and
tokenized audio sequence $\mathbf{s}$ into the decoder-only language model
(with parameters denoted as $\theta$) , estimating the probability
autoregressively, formulated as:

$p(\mathbf{s}|\mathbf{x};\theta)=\prod_{i=1}^{t}p(\mathbf{s}_{i}|\mathbf{s}_{:i-1},\mathbf{x};\theta)$
(2)

The model optimizes its parameters by performing the task of predicting the
next acoustic feature based on the phonemized text and historical speech
features. This training process can be viewed as maximizing the sequence
probability of the tokenized audio sequence under the condition of the text,
which can be described as:

$\hat{\theta}=\arg\max\limits_{\theta}p(\mathbf{s}|\mathbf{x};\theta)=\arg\max\limits_{\theta}\prod_{i=1}^{t}p(\mathbf{s}_{i}|\mathbf{s}_{:i-1},\mathbf{x};\theta)$
(3)

where $\hat{\theta}$ is our optimization target.

####  2.2.2 Inference with prompts

During inference, the model generates speech with similar speakers and prosody
as the given audio prompt. Concatenating tokenized audio prompt
$\mathbf{s}_{\text{prompt}}$ , phonemized source text $\mathbf{x}$ , and the
phonemized text prompt $\mathbf{x}_{\text{prompt}}$ corresponding to the audio
prompt, we build the input of the spoken language model in format of
$\\{\mathbf{x}_{\text{prompt}},\mathbf{x},\mathbf{s}_{\text{prompt}}\\}$ . The
model then autoregressively generates the remaining sequence after this
sequence. The decoding process consists of acoustic BPE decoding (if used) and
the vocoding process of a discrete-unit-based vocoder [20] for speech waveform
synthesis. Additionally, the Mel-spectrogram extracted from the speech prompt
is also input into the vocoder to achieve better speaker control.

##  3 Experiments and Results

In this section, we will introduce the exploration involving the configuration
of various dimensions in the acoustic BPE method and analyze the results
regarding Decoder-only TTS performance with these configurations.

###  3.1 Experimental Setup

####  3.1.1 Datasets

The experiments were conducted on the LibriSpeech [23] and LibriTTS [21]
datasets. LibriTTS, designed for text-to-speech tasks, consists of
approximately 585 hours of English speech from multiple speakers. Its
train-960 subset was used for model training. The LibriSpeech dataset served
primarily as the test set for speech synthesis. In this experiment, sentences
shorter than 4 seconds or longer than 10 seconds were filtered out from the
test-clean subset of LibriSpeech. The remaining 1145 sentences, totaling
approximately 2.02 hours of speech, were used as test cases. All speech audio
data was downsampled to 16kHz before use.

####  3.1.2 Settings of acoustic BPE

The encoding process of acoustic BPE consists of two stages: speech
discretization and acoustic BPE model training. In the speech discretization
stage, we utilized the HuBERT-
large††https://github.com/facebookresearch/fairseq/tree/main/examples/hubert
model pretrained with masked prediction on the 60k hours of LibriLight [24]
dataset and the WavLM-large††https://huggingface.co/microsoft/wavlm-large
model pretrained on extra dataset consisting of 10k hours of Gigaspeech [25]
and 24k hours English data subset of VoxPopuli [26], as the self-supervised
speech representation model. Continuous audio features were extracted from the
final layer's output activations of their Transformer encoder, to train
k-means models with 256, 2048, 4096 and 8192 centroids, which encoded the
semantic features extracted from the LibriTTS datasets into speech discrete
token sequences. Due to memory constraints, we randomly sampled 100 hours of
speech from the LibriTTS train-960 set during k-means model training.

In the next stage, we trained the acoustic BPE model on the speech discrete
token sequences extracted from the LibriTTS train-960 set. All discrete token
sequences were first converted into Unicode strings, and then a BPE model was
training on these strings with the
SentencePiece††https://github.com/google/sentencepiece toolkit. The pretrained
BPE model encoded the speech discrete token sequences from LibriTTS into
acoustic BPE sequences. Four different acoustic BPE encodings were
experimented with in this study: no acoustic BPE encoding, and acoustic BPE
encoding with vocabulary sizes of 5,000, 10,000, and 20,000 subwords,
respectively.

####  3.1.3 Decoder-only TTS architecture and training

The TTS model is based on the AR model of VALL-E [3], which is a decoder-only
Transformer architecture consisting of 12 Transformer layers, each with 16
attention heads and a hidden feature dimension of 1024. The absolute
positional encoding based on trigonometric functions is employed separately to
the text and speech feature sequences. As there are no official
implementations, we resorted to an unofficial
repository††https://github.com/lifeiteng/vall-e.

During training, the model was optimized using the ScaledAdam [27] optimizer
and the Eden [27] learning rate scheduler with an initial learning rate of
0.05. Each batch of training data contained approximately 100 seconds of audio
data, including only speech samples with duration between 1 second and 15
seconds. To achieve a larger effective batch size, gradients were accumulated
four times before each update. The model was trained for 20 epochs on an
NVIDIA A10 GPU.

####  3.1.4 Vocoder setup

In each setting of the SSL extractor and k-means clustering, we trained a
discrete-unit-based vocoder to convert semantic tokens into
waveform††https://github.com/X-LANCE/UniCATS-CTX-vec2wav. This vocoder
contains two conformer blocks each with 2 layers and 184 attention dimensions.
After the conformer blocks, a HifiGAN [28] is cascaded and the same learning
criterions apply. We trained all the replicas for 1M steps on LibriTTS, and
shared them within different acoustic BPE settings.

####  3.1.5 Evaluation metrics

To conduct a thorough investigation on the effectiveness of the acoustic BPE,
we established the evaluation process based on a range of metrics,
encompassing both subjective and objective dimensions. Contrasting with
scenarios where this method is absent, these metrics serve to provide an
intuitive way of gauging the impact of acoustic BPE. Below is an explanation
of the metrics utilized.

Speech Intelligibility: We used a publicly available automatic speech
recognition (ASR) model based on the Conformer-Transducer
architecture††https://huggingface.co/nvidia/stt_en_conformer_transducer_xlarge
to transcribe synthesized speech into text, and then computed the word error
rate (WER) between the transcribed text and the ground truth text.

Speech Quality and Naturalness: We employed naturalness Mean Option Score
(MOS) as the subjective metric of speech quality and naturalness. In the
listening test, each participant was assigned with multiple test cases, each
containing high-intelligibility speech synthesized by different settings from
the same sentence. Participants were asked to give a score of 1-5 to each
speech based on its naturalness. Mel cepstral distortion (MCD) [29] with
dynamic time warping (DTW), measuring the MFCC distance between the
synthesized and reference mel-spectrum features with DTW helping to find the
optimal time-sequence alignment, was utilized as the objective metrics.

Inference Speed: This research focuses on the inference speed of the decoder-
only language model, which is the main architecture of the TTS system. Real-
time factor (RTF) was measured to evaluate the inference speed.

Sample Diversity: Number of statistically-different bins (NDB) [30] and
Jensen-Shannon (JS) divergence were employed to objectively evaluate the
sample diversity. Detailed calculating process will be introduced in Section
3.3.1.

Table 1: Results of the decoder-only TTS model performance (WER, MOS) and inference speed (RTF) metrics under different BPE vocabulary sizes of HuBERT+kmeans (kms) 2048 centroids. Here ``resyn." means resynthesis from vocoder, while ``aBPE" refers to the BPE vocabulary size. Tokenizer | Task | aBPE | WER$\downarrow$ | MOS$\uparrow$ | RTF$\downarrow$
---|---|---|---|---|---
| resyn. | - | 1.9 | 4.39 $\pm$ 0.13 | -
|  | - | 9.1 | 4.13 $\pm$ 0.19 | 0.129
|  | 5000 | 5.7 |  4.20 $\pm$ 0.17 | 0.069
|  | 10000 | 5.5 | 4.19 $\pm$ 0.13 | 0.052
HuBERT kms 2048 | TTS | 20000 | 5.2 | 4.11 $\pm$ 0.19 | 0.045
Table 2: The objective metrics (WER and MCD) of synthetic speech under different SSL, k-means and acoustic BPE settings. Results are presented in the format: [w/o. BPE] / [w. BPE 10000].  SSL | HuBERT | WavLM
---|---|---
K-means | 256 | 2048 | 4096 | 8192 | 256 | 2048 | 4096 | 8192
WER $\downarrow$ | 6.0 / 6.3 | 9.1 / 5.5 | 7.7 / 4.5 | 7.8 / 8.0 | 9.9 / 6.0 | 8.3 / 7.8 | 11.4 / 9.0 | 9.0 / 9.5
MCD $\downarrow$ | 11.9 / 12.6 | 12.1 / 12.0 | 11.9 / 11.8 | 11.9 / 11.7 | 11.9 / 12.5 | 12.2 / 12.4 | 12.3 / 12.3 | 11.9 / 12.0

###  3.2 Acoustic BPE Enhances TTS Model Performance

After our preliminary experiments, it was found that in general situation,
incorporating the acoustic BPE method indeed leads to significant improvements
in the TTS model's ability to synthesize speech. We primarily focus on
intelligibility and quality of synthesized speech to reflect this ability. The
results of selected outperforming settings are presented in Table 1.

####  3.2.1 Improvement of speech integibility and quality

We conducted integibility experiments on the entire test set, which consists
of 1145 sentences. The WER of the synthesized audios ASR results are
calculated. As shown in Table 1, the intelligibility performance of audio
generated by TTS models using acoustic BPE method has a significant advantage
over those without using it. Under the pre-configuration (tokenizer composed
of HuBERT-large and a 2048-centroid kmeans model), the reduction is up to 3.9%
WER. The results of the MOS metric indicate that models using the acoustic BPE
method can generate competitive audio quality. In the worst-case scenario, it
is slightly lower by 0.02 compared to not using it, and it is accompanied by
lower confidence. However, in the best-case scenario, it can be higher by
0.07, accompanied by better confidence.

####  3.2.2 Acceleration of inference and training process

The calculation of RTF in this research is equivalently defined as the ratio
of the docoder-only LM inference time to the length of its input. As the
presented result, it significantly decreases as the expansion of its
vocabulary increase. This is due to the merging operations of BPE, which
shorten the sequence length of the input language model. The resulting
acceleration effect of the model far exceeds the time spent on the increased
computational cost of token embeddings due to the expansion of the vocabulary.
From an intuitive perspective, it results in up to 2.9 times inference
speedup. Additionally, we record the time cost of the whole training process,
the settings with acoustic BPE method also perform approximately 1.5 $\sim$ 2
times acceleration decided to configurations.

###  3.3 Enrichment of sample diversity

We are interested in investigating whether the acoustic BPE method can
influence the generation of more diverse intonations or speech rates when
reading the same sentence. We opted to utilize NDB and JS divergence as
metrics for assessing the diversity of samples in this TTS model. Experimental
results from various configurations are presented in Table 3.

####  3.3.1 The NDB and JS divergence metrics

The NDB metric is proposed in [30]. In our experiment, we extracted prosodic
features††Namely Kaldi-style [31] pitch, probability of voice and energy. from
selected low-WER synthesized utterances in the test set. Additionally, with
speaker information already present in prompts, the influence of semantic and
speaker information on sample diversity can be almost negligible. Samples,
considered as the prosodic features of frames, were divided into a training
and evaluation set. We first set a fixed number $k$ of cluster bins by
training a k-means model with samples in the training set, maintaining $k$
centroids. Let $p$ and $q$ denote the statistical distributions of samples in
the training and evaluation set on these bins, and let $m$ and $n$ denote the
sample numbers of the training and evaluation set, respectively. The
calculation of NDB metric was conducted with a two-proportion z-test on $p$,
$q$ and $m$, $n$, counting the number of the bins whose final p-values are
less than a manually set value of significant level, and dividing it by $k$.

The definition of Jensen-Shannon (JS) divergence is

$JS=\frac{1}{2}\left[KLD\left(p\mid\mid\frac{p+q}{2}\right)+KLD\left(q\mid\mid\frac{p+q}{2}\right)\right]$
(4)

where $KLD(\cdot||\cdot)$ is the Kullback–Leibler (KL) divergence function. To
reduce the impact of uncertain selection of two sample sets, we calculated two
metrics using a public tool††https://github.com/eitanrich/gans-n-
gmms/blob/master/utils/ndb.py for 10 times and used their average values as
the final results. Greater values of NDB and JS indicate a larger difference
between the statistical distributions of the samples, suggesting better
performance in terms of diversity.

Table 3: Sample diversity results. Here ``kms" means kmeans centroids, and ``aBPE" refers to the BPE vocabulary size. Tokenizer | aBPE | NDB$\uparrow$ | JS$\uparrow$
---|---|---|---
HuBERT kms 2048 | - | 0.649 | 0.00313
5000 | 0.676 | 0.00439
10000 | 0.668 | 0.00396
20000 | 0.678 | 0.00388
WavLM kms 2048 | - | 0.665 | 0.00364
5000 | 0.680 | 0.00415
10000 | 0.687 | 0.00408
20000 | 0.684 | 0.00453

(a) NDB

(b) JS

Figure 1: NDB and JS divergence results under varied semantic token and
acoustic BPE settings.

####  3.3.2 Increment of TTS sample diversity with acoustic BPE

According to the presented results, TTS models employing the acoustic BPE
method showed significant advantages in synthesized audio diversity compared
to those without it. Referring to Figure 1, this phenomenon is essentially
prevalent. We will discuss extreme boundary cases in detail in next section
3.4.

###  3.4 Discussion on boundary cases and limitations

Although the benefits of the acoustic BPE in improving model performance are
substantial, every method has its applicable range and limitations. This is
particularly observed in our study, where there are configurable parameters
that can be flexibly adjusted. Exploring the boundary conditions of this
method under our research conditions can provide intuitive insights for
constructing acoustic BPE methods in other decoder-only LM-based TTS systems
with complex configurations. According to the results in Table 2, both too
small and too large number of k-means centroids can lead to a plateau or even
a decline in model performance. In practical experiments, when there is a
significant gap between the number of k-means centroids and vocabulary size,
instability may arise, e.g. leading to the model repeatedly outputting the
same token. Moreover, WavLM also exhibits instability in the TTS architecture
used in this experiment. At this time, the use of acoustic BPE method may
exacerbate this instability, resulting in worse TTS performance.

##  4 Conclusion

In conclusion, the application of the acoustic BPE method in TTS tasks brings
significant performance benefits, including improvements in the
intelligibility, quality, and diversity of generated audio, as well as notable
enhancements in training and inference speed. This undoubtedly demonstrates
the potential of the method in discrete speech language modeling and speech-
text language modeling. However, the presence of the acoustic BPE vocabulary
increases the number of discrete tokens in speech. This imposes certain
limitations on its configuration selection. In future work, we will conduct
experiments with this method under scaled up datasets and models.
Additionally, we will explore other effective methods for tokenizing audio.

##  5 Acknowledgements

This work was supported by the China NSFC Project (No. 92370206), Shanghai
Municipal Science and Technology Major Project (2021SHZDZX0102) and
Development Program of Jiangsu Province, China (No.BE2022059).

## References

  * [1] K. Lakhotia, E. Kharitonov, W.-N. Hsu, Y. Adi, A. Polyak, B. Bolte, T.-A. Nguyen, J. Copet, A. Baevski, A. Mohamed, and E. Dupoux, ``On Generative Spoken Language Modeling from Raw Audio,'' _TACL_ , vol. 9, pp. 1336–1354, 12 2021.
  * [2] Z. Borsos, R. Marinier, D. Vincent, E. Kharitonov, O. Pietquin, M. Sharifi, D. Roblek, O. Teboul, D. Grangier, M. Tagliasacchi, and N. Zeghidour, ``Audiolm: A language modeling approach to audio generation,'' _IEEE/ACM TASLP_ , vol. 31, pp. 2523–2533, 2023.
  * [3] C. Wang, S. Chen, Y. Wu, Z. Zhang, L. Zhou, S. Liu, Z. Chen, Y. Liu, H. Wang, J. Li, L. He, S. Zhao, and F. Wei, ``Neural codec language models are zero-shot text to speech synthesizers,'' 2023.
  * [4] E. Kharitonov, D. Vincent, Z. Borsos, R. Marinier, S. Girgin, O. Pietquin, M. Sharifi, M. Tagliasacchi, and N. Zeghidour, ``Speak, Read and Prompt: High-Fidelity Text-to-Speech with Minimal Supervision,'' _TACL_ , vol. 11, pp. 1703–1718, 12 2023.
  * [5] X. Zhu, Y. Lv, Y. Lei, T. Li, L. Xie, W. HE, H. Zhou, and H. Lu, ``Vec-tok speech: Speech vectorization and tokenization for neural speech generation,'' 2024. [Online]. Available: https://openreview.net/forum?id=C53xlgEqVh
  * [6] M. Łajszczak, G. Cámbara, Y. Li, F. Beyhan, A. van Korlaar, F. Yang, A. Joly, Álvaro Martín-Cortinas, A. Abbas, A. Michalski, A. Moinet, S. Karlapati, E. Muszyńska, H. Guo, B. Putrycz, S. L. Gambino, K. Yoo, E. Sokolova, and T. Drugman, ``BASE TTS: Lessons from building a billion-parameter text-to-speech model on 100k hours of data,'' 2024.
  * [7] N. Zeghidour, A. Luebs, A. Omran, J. Skoglund, and M. Tagliasacchi, ``Soundstream: An end-to-end neural audio codec,'' _IEEE/ACM TASLP_ , vol. 30, pp. 495–507, 2022.
  * [8] A. Défossez, J. Copet, G. Synnaeve, and Y. Adi, ``High fidelity neural audio compression,'' _Transactions on Machine Learning Research_ , 2023, featured Certification, Reproducibility Certification.
  * [9] R. Kumar, P. Seetharaman, A. Luebs, I. Kumar, and K. Kumar, ``High-fidelity audio compression with improved rvqgan,'' in _Advances in Neural Information Processing Systems_ , A. Oh, T. Neumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, Eds., vol. 36. Curran Associates, Inc., 2023, pp. 27 980–27 993.
  * [10] A. Baevski, S. Schneider, and M. Auli, ``vq-wav2vec: Self-supervised learning of discrete speech representations,'' in _ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_ , 2020.
  * [11] W.-N. Hsu, B. Bolte, Y.-H. H. Tsai, K. Lakhotia, R. Salakhutdinov, and A. Mohamed, ``HuBERT: Self-supervised speech representation learning by masked prediction of hidden units,'' _IEEE/ACM TASLP_ , vol. 29, pp. 3451–3460, 2021.
  * [12] Y.-A. Chung, Y. Zhang, W. Han, C.-C. Chiu, J. Qin, R. Pang, and Y. Wu, ``w2v-bert: Combining contrastive learning and masked language modeling for self-supervised speech pre-training,'' in _2021 IEEE ASRU_ , 2021, pp. 244–250.
  * [13] S. Chen, C. Wang, Z. Chen, Y. Wu, S. Liu, Z. Chen, J. Li, N. Kanda, T. Yoshioka, X. Xiao, J. Wu, L. Zhou, S. Ren, Y. Qian, Y. Qian, J. Wu, M. Zeng, X. Yu, and F. Wei, ``WavLM: Large-scale self-supervised pre-training for full stack speech processing,'' _IEEE JSTSP_ , vol. 16, no. 6, pp. 1505–1518, 2022.
  * [14] F. Shen, Y. Guo, C. Du, X. Chen, and K. Yu, ``Acoustic BPE for speech generation with discrete tokens,'' 2024.
  * [15] S. Ren, S. Liu, Y. Wu, L. Zhou, and F. Wei, ``Speech pre-training with acoustic piece,'' in _Interspeech 2022, 23rd Annual Conference of the International Speech Communication Association, Incheon, Korea, 18-22 September 2022_ , H. Ko and J. H. L. Hansen, Eds. ISCA, 2022, pp. 2648–2652.
  * [16] J. Devlin, M. Chang, K. Lee, and K. Toutanova, ``BERT: pre-training of deep bidirectional transformers for language understanding,'' in _Proc. NAACL-HLT 2019, Minneapolis, MN, USA, June 2-7, 2019_ , J. Burstein, C. Doran, and T. Solorio, Eds. Association for Computational Linguistics, 2019, pp. 4171–4186.
  * [17] X. Chang, B. Yan, Y. Fujita, T. Maekaku, and S. Watanabe, ``Exploration of Efficient End-to-End ASR using Discretized Input from Self-Supervised Learning,'' in _Proc. INTERSPEECH 2023_ , 2023, pp. 1399–1403.
  * [18] S. Maiti, Y. Peng, S. Choi, J.-W. Jung, X. Chang, and S. Watanabe, ``Voxtlm: Unified decoder-only models for consolidating speech recognition, synthesis and speech, text continuation tasks,'' in _ICASSP 2024_ , 2024, pp. 13 326–13 330.
  * [19] C. Du, Y. Guo, X. Chen, and K. Yu, ``Vqtts: High-fidelity text-to-speech synthesis with self-supervised vq acoustic feature,'' _arXiv preprint arXiv:2204.00768_ , 2022.
  * [20] C. Du, Y. Guo, F. Shen, Z. Liu, Z. Liang, X. Chen, S. Wang, H. Zhang, and K. Yu, ``UniCATS: A unified context-aware text-to-speech framework with contextual VQ-diffusion and vocoding,'' _CoRR_ , vol. abs/2306.07547, 2023.
  * [21] H. Zen, V. Dang, R. Clark, Y. Zhang, R. J. Weiss, Y. Jia, Z. Chen, and Y. Wu, ``Libritts: A corpus derived from librispeech for text-to-speech,'' in _Interspeech 2019_ , G. Kubin and Z. Kacic, Eds. ISCA, 2019, pp. 1526–1530.
  * [22] R. Sennrich, B. Haddow, and A. Birch, ``Neural machine translation of rare words with subword units,'' in _ACL 2016, August 7-12, 2016_.
  * [23] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur, ``Librispeech: An ASR corpus based on public domain audio books,'' in _Proc. ICASSP 2015_. IEEE, 2015, pp. 5206–5210.
  * [24] J. Kahn, M. Rivière, W. Zheng, E. Kharitonov, Q. Xu, P. Mazaré, J. Karadayi, V. Liptchinsky, R. Collobert, C. Fuegen, T. Likhomanenko, G. Synnaeve, A. Joulin, A. Mohamed, and E. Dupoux, ``Libri-light: A benchmark for asr with limited or no supervision,'' in _Proc. ICASSP 2020_ , 2020, pp. 7669–7673.
  * [25] G. Chen, S. Chai, G. Wang, J. Du, W. Zhang, C. Weng, D. Su, D. Povey, J. Trmal, J. Zhang, M. Jin, S. Khudanpur, S. Watanabe, S. Zhao, W. Zou, X. Li, X. Yao, Y. Wang, Z. You, and Z. Yan, ``GigaSpeech: An evolving, multi-domain ASR corpus with 10, 000 hours of transcribed audio,'' in _Interspeech 2021, Brno, Czechia 2021_. ISCA, 2021, pp. 3670–3674.
  * [26] C. Wang, M. Rivière, A. Lee, A. Wu, C. Talnikar, D. Haziza, M. Williamson, J. M. Pino, and E. Dupoux, ``Voxpopuli: A large-scale multilingual speech corpus for representation learning, semi-supervised learning and interpretation,'' in _Proc. ACL/IJCNLP 2021, August 1-6, 2021_ , pp. 993–1003.
  * [27] Z. Yao, L. Guo, X. Yang, W. Kang, F. Kuang, Y. Yang, Z. Jin, L. Lin, and D. Povey, ``Zipformer: A faster and better encoder for automatic speech recognition,'' in _12th ICLR_ , 2024.
  * [28] J. Kong, J. Kim, and J. Bae, ``HiFi-GAN: Generative adversarial networks for efficient and high fidelity speech synthesis,'' in _NeurIPS 2020, December 6-12_ , 2020.
  * [29] R. Kubichek, ``Mel-cepstral distance measure for objective speech quality assessment,'' in _Proceedings of IEEE Pacific Rim Conference on Communications Computers and Signal Processing_ , vol. 1, 1993, pp. 125–128 vol.1.
  * [30] E. Richardson and Y. Weiss, ``On GANs and GMMs,'' in _Advances in Neural Information Processing Systems_ , S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, Eds., vol. 31. Curran Associates, Inc., 2018.
  * [31] D. Povey, A. Ghoshal, G. Boulianne, L. Burget, O. Glembek, N. Goel, M. Hannemann, P. Motlicek, Y. Qian, P. Schwarz, J. Silovsky, G. Stemmer, and K. Vesely, ``The Kaldi speech recognition toolkit.'' IEEE Signal Processing Society, 2011, IEEE Catalog No.: CFP11SRW-USB.